Add r.squared and adj.r.squared to the phylolm result

[jmaspons]

Code based in caper::pgls and caper::summary.pgls

I can update the manual if you like the idea

[jmaspons]

It produce an error with the example in the manual.

Error in solve.default(comp$XX) : Lapack routine dgesv: system is exactly singular: U[2,2] = 0

Any idea on how to solve it?

[lamho86]

Hi @jmaspons,

I am out of town for a conference right now. I will have a look when I return. Thanks for the contribution.

Best, Lam

[jmaspons]

Fixed! It's working now.

We can decide to save r.squared and adj.r.squared in the phylolm object or to save NMS, NSSQ, RMS and RSSQ and calculate r.squared in the summary.phylolm (caper does the latter)

• [x] update man/phylolm.Rd
• [x] update summary.phylolm, print.summary.phylolm and print.phylolm

[jmaspons]

From my part it's ready to merge

[lamho86]

@jmaspons sigma2 in phylolm is estimated using MLE. So, your code is not quite correct. Also, the sigma2 of OU models is a little different than other models. I made some corrections. Could you please double-check for me? Thanks.

[jmaspons]

I compared the R² from phylolm and caper::pgls. Your changes seems to reduce the differences among functions for models optimizing lambda. So your changes seems correct. I can't talk about OU models.

datL<- replicate(100, {
  tre = rcoal(60)
  taxa = sort(tre$tip.label)
  b0=0
  b1=1
  x <- rTrait(n=1, phy=tre,model="BM",
              parameters=list(ancestral.state=0,sigma2=10))
  y <- b0 + b1*x + 
    rTrait(n=1,phy=tre,model="lambda",parameters=list(
      ancestral.state=0,sigma2=1,lambda=0.5))
  dat = data.frame(trait=y[taxa],pred=x[taxa])
}, simplify=FALSE)

## phylolm
fits.phylolm<- lapply(datL, function(x) phylolm::phylolm(trait ~ pred, data=x, phy=tre, model="lambda"))

## pgls
fits.pgls<- lapply(datL, function(x){
  x$species<- rownames(x)
  compDat<- caper::comparative.data(tre, x, names.col="species")
  fit.pgls<- caper::pgls(trait ~ pred, data=compDat, lambda="ML")
})

r.sqr.phylolm<- t(sapply(fits.phylolm, function(x) c(x$r.squared, x$adj.r.squared)))
r.sqr.pgls<- t(sapply(fits.pgls, function(x) c(summary(x)$r.squared, summary(x)$adj.r.squared)))

colnames(r.sqr.phylolm)<- colnames(r.sqr.pgls)<- c("r.sqr", "adj.r.sqr")

summary(r.sqr.phylolm - r.sqr.pgls)
## Implementation as in this pull request
#     r.sqr             adj.r.sqr        
# Min.   :0.0002515   Min.   :0.0002559  
# 1st Qu.:0.0012854   1st Qu.:0.0013076  
# Median :0.0022890   Median :0.0023285  
# Mean   :0.0028946   Mean   :0.0029445  
# 3rd Qu.:0.0037209   3rd Qu.:0.0037851  
# Max.   :0.0114158   Max.   :0.0116127 

## With commit 941e46a
#     r.sqr              adj.r.sqr         
# Min.   :-2.970e-07   Min.   :-3.020e-07  
# 1st Qu.:-1.600e-08   1st Qu.:-1.600e-08  
# Median : 0.000e+00   Median : 0.000e+00  
# Mean   : 4.030e-05   Mean   : 4.100e-05  
# 3rd Qu.: 2.500e-08   3rd Qu.: 2.500e-08  
# Max.   : 4.029e-03   Max.   : 4.099e-03 

[lamho86]

Thanks.

[lamho86]

@jmaspons I wonder if you would like to be added as a contributor to the package. If yes, please let me know your name. Best, Lam

[jmaspons]

Sure! person("Joan", "Maspons", email="j.maspons@creaf.uab.cat", role=c("ctb"))

Thanks