##Orthodont data from Pinheiro and Bates (2000) revisited
require(nlme)
m1 <- gls(distance ~ age, correlation = corCompSymm(value = 0.5, form = ~ 1 | Subject),
data = Orthodont, method= "ML")
##compare against lme fit
logLik(m1)
logLik(lme(distance ~ age, random = ~1 | Subject, data = Orthodont,
method= "ML"))
##both are identical
##compute predictions and SE's for different ages
predictSE(m1, newdata = data.frame(age = c(8, 10, 12, 14)))
##contagious bovine pleuropneumonia example modified from lme4
require(lme4)
data(cbpp)
##create proportion of incidence
cbpp$prop <- cbpp$incidence/cbpp$size
gm1 <- glmer(prop ~ period + (1 | herd), family = binomial,
weights = size, data = cbpp)
##create a data set to make predictions
newherd<- data.frame(period = as.factor(c("1", "2", "3", "4")))
##predictions on logit link scale
predictSE(mod = gm1, newdata = newherd, se.fit = TRUE,
type = "link", level = 0, print.matrix = FALSE)
##predictions on scale of original response variable
predictSE(mod = gm1, newdata = newherd, se.fit = TRUE,
print.matrix = TRUE, level = 0, type = "response")
if(require(unmarked)) {
##example with mallard data set from unmarked package
data(mallard)
mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
obsCovs = mallard.obs)
##run model with zero-inflated Poisson abundance
fm.mall.one <- pcount(~ ivel + date ~ length + forest, mallardUMF, K=30,
mixture = "ZIP")
##make prediction
predictSE(fm.mall.one, type = "response", parm.type = "lambda",
newdata = data.frame(length = 0, forest = 0, elev = 0))
##compare against predict
predict(fm.mall.one, type = "state", backTransform = TRUE,
newdata = data.frame(length = 0, forest = 0, elev = 0))
##add offset in model to scale abundance per transect length
fm.mall.off <- pcount(~ ivel + date ~ forest + offset(length), mallardUMF, K=30,
mixture = "ZIP")
##make prediction
predictSE(fm.mall.off, type = "response", parm.type = "lambda",
newdata = data.frame(length = 10, forest = 0, elev = 0))
##compare against predict
predict(fm.mall.off, type = "state", backTransform = TRUE,
newdata = data.frame(length = 10, forest = 0, elev = 0))
detach(package:unmarked)
}
Run the code above in your browser using DataLab