I'm using a negative binomial glmm for analysis.
The y values of my samples are summed over 12 individual counts of insect mines on leaves of a tree. You can get the data (and session info) from here!.
Treat has 2 levels; Year has 2 levels; Habitat has 3 levels; leaves is the same number of 12 leaves for each sample; location has 4 levels and is almost similar to location. Two of the locations fall in the same "habitat" category. I use sublocation as random effect, since I want to get results for habitat not for all sublocations separately while being aware of the more similar measures at a single location.
This is the model:
library(MASS) library(lme4) library(optimx) glmm1 <- glmer.nb(response ~ year*(Habitat+Treat) + offset(log(leaves)) + (1|location), data = df1, control = glmerControl(optimizer = "optimx", calc.derivs = FALSE, optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
Now I use the following code to create a prediction plot:
# Prediction plot predframe <- expand.grid(year= levels(df1$year), Habitat= levels(df1$Habitat), Treat= levels(KMM2$Treat), leaves = 10) # Set the offset!!
Here I am able to alter the offset to get results for the more commonly used reference of 10 leaves !
mm <- model.matrix(~ year*(Habitat + Treat) + offset(log(leaves)), data = predframe) predframe$fit <- predict(glmm1,newdata=predframe, re.form = NA) # Don't use type="response"!!! # Alternatively use: predframe$fit <- mm %*% fixef(glmm1) pvar1 <- diag(mm %*% tcrossprod(vcov(glmm1),mm)) # alternatively use: pvar1 <- diag(mm %*%vcov(glmm1) %*% t(mm)) tvar1 <- pvar1+VarCorr(glmm1)$Standort # must be adapted for more complex models predframe <- data.frame(predframe , p_lwr = predframe$fit-1.96*sqrt(pvar1) # for confidence interval , p_upr = predframe$fit+1.96*sqrt(pvar1) , t_lwr = predframe$fit-1.96*sqrt(tvar1) # for prediciton interval , t_upr = predframe$fit+1.96*sqrt(tvar1) ) predframe2 <- aggregate(cbind(fit, p_lwr, p_upr) ~ year + Treat,predframe, mean) ggplot(aes(y=interaction(year, Treat), x=exp(fit)), data=predframe2) + geom_point() + geom_errorbarh(aes(xmax=exp(p_upr), xmin=exp(p_lwr)))
Here I use the lsmeans package to do the same:
library(lsmeans) lsm.options(save.ref.grid = TRUE) lsm <- lsmeans(glmm1, ~ Netz*Jahr, offset=10) .Last.ref.grid
Offset was not adjusted to 10!
Besides slightly different confidence intervals (why?) these results are the same as if I used leaves=12 for the
How can I adjust the offset in lsmeans?
Thanks in advance