## Solution to challenge problem 2: GLS vs. PIC regression

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

Then our data from file:

``````datos<-read.csv("Barbetdata.csv",header=TRUE,row.names=1)
``````

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

``````obj<-name.check(arbol,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.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.