two versions of OU in phylolm

[richfitz]

Issue by mwpennell from Tuesday Dec 10, 2013 at 19:48 GMT
Originally opened as traitecoevo/modeladequacy#63

---

this is a very low priority at the moment but there are two different options for OU models in phylolm -- depending on whether one want to assume stationarity or not. i am just going to ignore the difference for now and rescale them in the same way. just wanted to put the issue on record

[richfitz]

Comment by richfitz from Wednesday Dec 18, 2013 at 02:31 GMT

---

This should be simlar to the with.optimum=TRUE / with.optimum=FALSE models made by diversitree::make.ou. It's possible that the with optimum scaling could be tricky to get right.

[richfitz]

Comment by richfitz from Wednesday Dec 18, 2013 at 02:39 GMT

---

Actually, looking further into this, phylom.R contains:

        if (model=="OUrandomRoot") {
            alpha = value
            edge.length = numeric(N)
            distFromRoot <-  exp(-2*alpha*times)
            externalEdge <- des<=n
            for (i in 1:N) {
                d1 = distFromRoot[which(names(times)==anc[i])]
                if (externalEdge[i]) {d2 = exp(-2*alpha*D[des[i]])} else d2 = distFromRoot[which(names(times)==des[i])]
                edge.length[i] = d2 - d1
                }
            root.edge = min(distFromRoot)
            }

        ### OUfixedRoot model   
        if (model=="OUfixedRoot") {
            alpha = value
            edge.length = numeric(N)
            distFromRoot <-  exp(-2*alpha*times)*(1 - exp(-2*alpha*(Tmax-times)))
            externalEdge <- des<=n
            for (i in 1:N) {
                d1 = distFromRoot[which(names(times)  == anc[i])]
                if (externalEdge[i]) {d2 = exp(-2*alpha*D[des[i]])*(1-exp(-2*alpha*(Tmax-D[des[i]])))}
                    else d2 = distFromRoot[which(names(times)==des[i])]
                edge.length[i] = d2 - d1
                }
            root.edge = min(distFromRoot)
            }

which look to be equivalent to different root treatments only. The relevant lines are

# random, fixed:
distFromRoot <-  exp(-2*alpha*times)
distFromRoot <-  exp(-2*alpha*times)*(1 - exp(-2*alpha*(Tmax-times)))

# external edge: random, fixed:
exp(-2*alpha*D[des[i]])
exp(-2*alpha*D[des[i]])*(1-exp(-2*alpha*(Tmax-D[des[i]])))

# internal edge: random, fixed
distFromRoot[which(names(times)==des[i])]
distFromRoot[which(names(times)==des[i])]

which are not so different. I'm not sure I see how these relate to the root treatment actually.

I agree that this is not a priority. We could probably use the transf.branch.lengths function to do the hard work creating a unit tree anyway, without having to pick apart the models.