Today I got the following interesting user request:

*"I have been playing around with the function pgls.Ives. Is there anyway to accomodate the corPagel correlation structures into this?"*

This is not possible, actually. corStruct objects (such as an object of class "corPagel") do not permit (so far as I know) error in both x & y which is a component of the pgls.Ives model. However we can, with a little pain, optimize Pagel's λ model using pgls.Ives - just in a slightly roundabout way. We do this by building a simple wrapper around pgls.Ives that spits out the log-likelihood conditioned on a specific input λ, then we optimize for λ. Here's a demo:

require(phytools)

require(geiger)

# simulate tree & data under the regression model

tree<-pbtree(n=50)

xhat<-fastBM(tree)

yhat<-0.75*xhat+0.2*fastBM(transform(tree,lambda=0.5))

# simulate individual data for x & y

x<-sampleFrom(xhat,randn=c(1,10))

n<-summary(as.factor(names(x)))

y<-sampleFrom(yhat,n=n[names(yhat)])

lambda<-seq(0.05,1,0.05)

logL<-sapply(lambda,function(l,tree,x,y)

pgls.Ives(transform(tree,lambda=l),x,y)$logL,tree=tree,

x=x,y=y)

ii<-which(logL==max(logL))

plot(lambda,logL,type="b")

lines(c(lambda[ii],lambda[ii]),c(min(logL),max(logL)),

lty="dashed")

ml.lambda<-lambda[ii]

ml.lambda

## optimize using base function optimize

foo<-function(l,tree,x,y)

pgls.Ives(transform(tree,lambda=l),x,y)$logL

optimize(foo,c(0,1),tree=tree,x=x,y=y,maximum=TRUE)

Let's try it:

> tree<-pbtree(n=50)

> xhat<-fastBM(tree)

> yhat<-0.75*xhat+0.2*fastBM(transform(tree,lambda=0.5))

>

> # simulate individual data for x & y

> x<-sampleFrom(xhat,randn=c(1,10))

> n<-summary(as.factor(names(x)))

> y<-sampleFrom(yhat,n=n[names(yhat)])

>

> ## brute force grid search optimization

> lambda<-seq(0.05,1,0.05)

> logL<-sapply(lambda,function(l,tree,x,y)

pgls.Ives(transform(tree,lambda=l),x,y)$logL,tree=tree,

x=x,y=y)

There were 50 or more warnings (use warnings() to see the first 50)

> ii<-which(logL==max(logL))

> ml.lambda<-lambda[ii]

> ml.lambda

[1] 0.7

> plot(lambda,logL,type="b")

> lines(c(lambda[ii],lambda[ii]),c(min(logL),max(logL)),

lty="dashed")

> foo<-function(l,tree,x,y) pgls.Ives(transform(tree,lambda=l),x,y)$logL

> optimize(foo,c(0,1),tree=tree,x=x,y=y,maximum=TRUE)

$maximum

[1] 0.6893112

$objective

[1] -134.8

There were 42 warnings (use warnings() to see them)

Note that although this works fine in principle, I have received some reports of the basic function failing to optimize for larger trees - and I believe these reports, so please use pgls.Ives with considerable caution until I work out these issues!

## No comments:

## Post a Comment