library(phylolm)

data = read.table("larinae.txt", header = T)
row.names(data) = data$Species
laritree = read.nexus("lari.nex")
data$ALat = abs(data$CLat)
data$stdWL = (data$WingLoading-mean(data$WingLoading, na.rm = T))/sd(data$WingLoading, na.rm = T)
laritree = drop.tip(laritree,laritree$tip.label[-match(data$Species, laritree$tip.label)])
#plot(laritree, cex = 0.75)

# no obvious relationship between the two predictors
cor.test(data$stdWL,data$ALat) # R = 0.11
plot(data$stdWL,data$ALat)

# using mid wingspan versus min or max is justified
cor.test(data$WSmin,data$Wingspan) # R = 0.98
cor.test(data$WSmax,data$Wingspan) # R = 0.98


## Models for comparison

modelN = lm(MantleDarkness~stdWL+ALat, data=data) # non-phylo
summary(modelN)
phmn = logLik(modelN); phmn

modelB = phylolm(MantleDarkness~stdWL+ALat, data=data, phy=laritree) # brownian
summary(modelB)
phmb = logLik(modelB); phmb

modelOU = phylolm(MantleDarkness~stdWL+ALat, model = "OUfixedRoot", data=data, phy=laritree) # ornstein-uhlenbeck (fixed root and random root give equivalent estimates)
summary(modelOU)
phmo = logLik(modelOU); phmo
halflife = log(2)/27.12717; halflife

modelL = phylolm(MantleDarkness~stdWL+ALat, model = "lambda", data=data, phy=laritree) # pagel's lambda
summary(modelL)
phml = logLik(modelL); phml


## Likelihood ratio tests

test1 = -2*(-128.21--125.1501) # no phylo vs ou
pchisq(test1, df = 1, lower.tail = F) # no phylo model sig. diff. from ou

test2 = -2*(-133.5707--125.1501) # brownian vs ou
pchisq(test2, df = 1, lower.tail = F) # ou model sig. diff. from brownian

test3 = -2*(-133.5707--127.0854) # brownian model sig. diff. from lambda
pchisq(test3, df = 1, lower.tail = F)

test4 = -2*(-128.21--127.0854) # no phylo vs lambda
pchisq(test4, df = 1, lower.tail = F) # no phylo NOT sig. diff. from lambda


# Relationship between wing loading and aspect ratio
model_ARWL = lm(AspectRatio~stdWL, data=data)
summary(model_ARWL)


#### Checking for bias in the Larinae phylogenetic tree

## With thanks to Kevin Arbuckle for providing the code

# simulate traits evolving along the Larinae tree under Brownian Motion and Ornstein-Uhlenbeck processes
# obtain difference in AIC for each pair
AreOuOK<-function(tree,nsim=500){
  
  traits<-as.data.frame(rTrait(n=nsim,tree,model="BM"))
  
  bias<-c()
  
  for (i in 1:nsim){
    
    bm<-phylolm(traits[,i]~1,data=traits,tree,model="BM")
    
    ou<-phylolm(traits[,i]~1,data=traits,tree,model="OUfixedRoot")
    
    bias[i]<-bm$aic-ou$aic
    
  }
  
  return(bias)
  
}

obj = AreOuOK(laritree)

# visually inspect results - values should be below 0
hist(obj,main="",xlab="AIC(BM)-AIC(OU)",col="light grey",xlim=c(min(obj)-1,max(obj)+1))
# mainly correctly identifies BM process, but some simulations create large AIC values for BM relative to OU

quantile(obj,probs=c(0.025,0.975))

# compare with observed AIC values
b = AIC(modelB)
o = AIC(modelOU)
x = b - o
table(obj>=x)["TRUE"]/table(obj>=x)["FALSE"]
# p = NA (i.e. very close to 0)
# This indicates support for the OU model over the BM model