Extraction of ML and REML s2 estimates from gls models

[richfitz]

I've taken a shot at this (see 79c2650), but I'm not sure that things are correct.

library(arbutus)
library(diversitree)
library(nlme)
set.seed(1)
phy <- tree.bd(pars=c(1,0), max.taxa=100)
tx <- sim.character(phy, 1)
ty <- sim.character(phy, 1) + 3
data <- data.frame(x=tx, y=ty, row.names=names(tx))

fit.gls.bm.ml   <- gls(y ~ x, data, corBrownian(1, phy), method="ML")
fit.gls.bm.reml <- gls(y ~ x, data, corBrownian(1, phy), method="REML")

# Likelihood function for the same model
lik.pgls.bm <- make.pgls(phy, y ~ x, data)

# Extract s2 from the fitted object and 
s2.ml <- arbutus:::estimate.sigma2.gls(fit.gls.bm.ml)
s2.reml <- arbutus:::estimate.sigma2.gls(fit.gls.bm.reml)
p.ml <- c(coef(fit.gls.bm.ml), s2=s2.ml)

# Run a ML fit with diversitree
fit.pgls.bm <- find.mle(lik.pgls.bm, p.ml)

n <- fit.gls.bm.ml$dims$N # 100, or length(phy$tip.label)
k <- fit.gls.bm.ml$dims$p # 2, or length(coef(fit.gls.bm.ml))
s2 <- coef(fit.pgls.bm)[[3]] # diversitree's s2 estimate

# Not the same as that from arbutus, and half way between REML and ML estimates
s2 - c(s2.ml, s2.reml)

# This is the relationship with the REML estimate
s2 * n / (n - (k - 1)) - s2.reml

which all suggests that in arbutus:::estimate.sigma2.gls() we should use
(n - (k - 1)) / n as the scaling factor. This doesn't quite match my view of how these work, but I might be being dense.

Ideas?

(this depends on very recent versions of both arbutus and diversitree)

[mwpennell]

sorry i missed this one -- shouldn't have opened a new issue. i am not quite sure about why this is but the (n-(k-1))/n seems to be correct. i will change this.

[richfitz]

Yeah, it seems empirically to be correct. Is there any reference for this? I had a quick look through Felsenstein's book but didn't see much there.

[mwpennell]

i will go look this up. suspect if it is anywhere, it is in Blomberg's 2012 paper.

[richfitz]

Migrated to mwpennell/arbutus#9