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[1] # 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!
plot(cld(lsm, type="response"))
Besides slightly different confidence intervals (why?) these results are the same as if I used leaves=12 for the predframe
above
How can I adjust the offset in lsmeans?
Thanks in advance
ksmeans
function does not have anoffset
argument. That might explain why it has no effect. Meanwhile, the offset in the model formula is included in the calculations. To see what is used, look atlsm@grid
$\endgroup$ – Russ Lenth Dec 23 '16 at 17:47lsmeans
where it saysksmeans
. $\endgroup$ – Russ Lenth Dec 23 '16 at 22:36at
to specify differentleaves
value(s). Note thatleaves
is one of the variables in the reference grid, and it is used to calculate the offset. $\endgroup$ – Russ Lenth Dec 27 '16 at 14:41