0
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Well, the ksmeans function does not have an offset 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 at lsm@grid $\endgroup$ – rvl Dec 23 '16 at 17:47
  • $\begingroup$ Oops, typo in above. Of course I meant lsmeans where it says ksmeans. $\endgroup$ – rvl Dec 23 '16 at 22:36
  • $\begingroup$ Hello, thanks for your answer! Yes, you are right, the reference grid uses the offset variable (.Last.ref.grid). I knew this. But, if I understand you right, the offset value can't be changed? My case is simple, I've got the same offset value for all samples, but how is it when sb uses different values? the use of an offset might be useless then? I used offset to adjust predictions when samples are from different volumes, plot size, etc. $\endgroup$ – Pharcyde Dec 27 '16 at 9:13
  • $\begingroup$ You can use at to specify different leaves value(s). Note that leaves is one of the variables in the reference grid, and it is used to calculate the offset. $\endgroup$ – rvl Dec 27 '16 at 14:41
1
$\begingroup$

This is how it works, thanks to rvl!!

lsm.options(save.ref.grid = TRUE)
lsm <- lsmeans(glmm1, ~ Treat*year, at=list(leaves=10))
.Last.ref.grid
plot(cld(lsm, type="response", sort=FALSE))

Now the reference grid is similar to the prediction frame above.

However: same same but different! Here are the two plots from both different procedures. I PREFER THE ONE FROM lsmeans!!!!!!

enter image description here

enter image description here

$\endgroup$
  • $\begingroup$ plot(lsm, type = "response") will do the same thing more efficiently, because the compact-letter-display calculations don't play a role in the resulting plot. $\endgroup$ – rvl Dec 27 '16 at 19:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.