You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
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 estimatess2 - c(s2.ml, s2.reml)# This is the relationship with the REML estimates2 * 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)
The text was updated successfully, but these errors were encountered:
Comment by mwpennell from Tuesday Dec 10, 2013 at 00:26 GMT
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.
Issue by richfitz from Wednesday Dec 04, 2013 at 05:39 GMT
Originally opened as traitecoevo/modeladequacy#56
I've taken a shot at this (see 79c2650), but I'm not sure that things are correct.
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)
The text was updated successfully, but these errors were encountered: