Solution to challenge problem 2: GLS vs. PIC regression

We start by loading our packages:

library(ape)
library(nlme)
library(geiger)

Then our data from file:

datos<-read.csv("Barbetdata.csv",header=TRUE,row.names=1)
arbol<-read.nexus("BarbetTree.nex")

Run name check to verify that we have the same data in the tree and data set:

obj<-name.check(arbol,datos)
arbol.cortado<-drop.tip(arbol, obj$tree_not_data)
name.check(arbol.cortado, datos)

## [1] "OK"

Now let's make a correlation structure for our PGLS model and fit it:

bm<-corBrownian(1,arbol.cortado)
pgls.model<-gls(wing~Lnalt,data=datos,correlation=bm)
summary(pgls.model)

## Generalized least squares fit by REML
##   Model: wing ~ Lnalt 
##   Data: datos 
##         AIC       BIC   logLik
##   -41.52099 -37.21902 23.76049
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## 
## Coefficients:
##                 Value  Std.Error   t-value p-value
## (Intercept)  4.808792 0.23025431 20.884700  0.0000
## Lnalt       -0.053410 0.03332457 -1.602729  0.1191
## 
##  Correlation: 
##       (Intr)
## Lnalt -0.879
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -0.8708169  0.3897004  0.6685024  0.9361152  2.7445358 
## 
## Residual standard error: 0.1986708 
## Degrees of freedom: 33 total; 31 residual

Next, we can fit our PIC regression:

wing<-setNames(datos$wing,rownames(datos))
Lnalt<-setNames(datos$Lnalt,rownames(datos))
pic.wing<-pic(wing,arbol.cortado)
pic.Lnalt<-pic(Lnalt,arbol.cortado)
pic.model<-lm(pic.wing~pic.Lnalt+0)
summary(pic.model)

## 
## Call:
## lm(formula = pic.wing ~ pic.Lnalt + 0)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.061160 -0.018176 -0.001805  0.025447  0.086057 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)
## pic.Lnalt -0.05341    0.03332  -1.603    0.119
## 
## Residual standard error: 0.03931 on 31 degrees of freedom
## Multiple R-squared:  0.07652,    Adjusted R-squared:  0.04673 
## F-statistic: 2.569 on 1 and 31 DF,  p-value: 0.1191

What have we found?

Written by Liam J. Revell. Last updated 28 Jun. 2016.