New version of pgls.Ives

I just posted a new version of the phytools function pgls.Ives, which is a function to fit the regression model of Ives et al. (2007) incorporating sampling error in the estimation of species means (e.g., here).

The main update in this version is that it can now accept individual data, rather than just species means, variances, and covariances. The form for the individual data is the same as fitBayes (e.g., here) - in other words, just long vectors with names equal to the species names. Obviously, some names will repeat - which is allowed in vectors.

There is nothing magic about the way this function works. It just takes the individual data, and then computes the within-species variances and covariances (which are themselves used to estimate within-species sampling variances and covariances). This is done using, for example, R base functions such as aggregate. For those in need of a quick primer, the following demo shows how aggregate can be used to compute within-species variances and covariances for a vector of values x containing individual values with names(x) being the species names, some of which repeat:

> # simulate tree & data
> tree<-pbtree(n=30)
> x<-fastBM(tree) # true means
> # individual values
> xe<-sampleFrom(xhat,randn=c(1,5))
> xe
       t13         t13         t13         t13         t25
-1.83130612 -2.08011272 -1.99096486 -3.26950066 -1.17626558
       t25         t26         t26         t20         t20
0.57125210 -1.90520230 -1.96856340         ...         ...

> # get species means
> xbar<-aggregate(xe,by=list(names(xe)),FUN=mean)
> xbar<-setNames(xbar[,2],xbar[,1])
> xbar
        t1         t10         t11         t12         t13
-1.71666562 -3.92387486 -3.78403128 -1.82452852 -2.29297109
       t14         t15         t16         t17         t18
-0.83735741 -0.66027704 -1.84528088        ...         ...

> # get species variances
> xvar<-aggregate(xe,by=list(names(xe)),FUN=var)
> xvar<-setNames(xvar[,2],xvar[,1])
> xvar
        t1         t10         t11         t12         t13
0.319766112          NA 1.330431973 0.420867485 0.434420333
       t14         t15         t16         t17         t18
0.910668871 0.985273462 0.560029491        ...         ...

Something that's apparent from the little simulation above, is that if there is not at least two samples for a particular species, then the species variance cannot be computed. In that case, the program merely assumes that those species have within-species variance equal to the mean within-species variance - but it will also spit a warning to that effect.

Computing the sampling covariance from individual data requires that the same individuals be sampled in both vectors. Thus, if the argument Cxy is not supplied, pgls.Ives will try to do it. If different individuals were used for x & y, then sampling covariance is probably negligible. In this case - if individual data is being supplied - then one should also give the function a named vector of zeroes for the argument Cxy to avoid any errors of spurious results.

Finally, here's a quick demo showing how to use individual data with pgls.Ives for both zero and non-zero error covariance:

> # load source code
> source("pgls.Ives.R")
>  
> # simulate tree & data under the regression model
> tree<-pbtree(n=30)
> xhat<-fastBM(tree)
> yhat<-0.75*xhat+0.2*fastBM(tree)
>  
> # simulate individual data for x & y
> # assuming Cxy!=0
> x<-sampleFrom(xhat,randn=c(1,10))
> n<-summary(as.factor(names(x)))
> y<-sampleFrom(yhat,n=n[names(yhat)])
>  
> # fit model
> pgls.Ives(tree,x,y)
$beta
[1] -1.193712  1.124268
$sig2x
[1] 0.5043893
$sig2y
[1] 0.219358
$a
[1] 3.111079 2.303972
$logL
[1] -74.45021
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Warning messages:
1: Some species contain only one sample. Substituting mean variance.
2: Some species contain only one sample. Substituting mean variance.
3: Some species contain only one sample. Substituting mean variance.

Download the latest version of this function on the phytools page - or wait a bit and it will be in the next CRAN release.