Test of ML vs. REML Lambda estimation

ml_reml_thunderdome

#Headers
require(geiger)
require(phytools)
require(nlme)

#Wrappers
get.lam <- function(x) attr(x$apVar, "Pars")[1]
get.vars <- function(n, lam){
    tree <- rcoal(n=n)
    x <- fastBM(geiger::rescale(tree, 'lambda', lambda = lam),sig2=0.5) # Brownian motion
    
    ml <- tryCatch(get.lam(gls(x ~ 1, correlation = corPagel(1, tree, fixed=F), method = "ML")), error=function(z) NA)
    reml <- tryCatch(get.lam(gls(x ~ 1, correlation = corPagel(1, tree, fixed=F), method = "REML")), error=function(z) NA)
    return(c(ml,reml))
}

#Simulate
set.seed(123456)
params <- expand.grid(c(33, 66, 100), c(0.25, 0.5, 0.75))
reps <- 30
x <- 1
results <- matrix(NA, nrow=nrow(params)*reps, ncol=2)
for(i in seq_len(nrow(params))){
    for(j in seq_len(reps)){
        results[x,] <- get.vars(params[i,1], params[i,2])
        x <- x+1
    }
}

#t-test of difference between the two
abs.diff <- abs(results - rep(params[,2], each=reps))
t.test(abs.diff[,1], abs.diff[,2], paired=TRUE)
diff(apply(abs.diff, 2, mean, na.rm=TRUE))
#...marginally significant, but yes, REML is better by 0.03

#Boxplot of how much the two vary across tree size
ml.diff <- abs.diff[,1]-abs.diff[,2]
boxplot(ml.diff ~ rep(params[,1], each=reps))
boxplot(ml.diff ~ paste(rep(params[,1], each=reps), rep(params[,2], each=reps)))
#...no consistent trend - increase in variance (obviously) in smaller trees
# so, potentially, REML just got lucky this time on the average