2
$\begingroup$

I use R with the MuMIn package for Multimodel inference. My global Model is

M.global <- glmer(y ~ c1 + I(c1^2) + c2 + c3 + c4 + (1|subject),data, family =poisson(link="log))

After that I use

M.dredge <- dredge(M.global, m.max=3, subset = dc(cpl,I(cpl^2)))

To calculate all possible models from the global model and rank them by AICc. m.max limits the number of predictor to 3 and dc is a dependency chain, so the function takes the squared c1 only if unsquared c1 is included. m.max is used since the total sample size is 30, and according to Harrels 10:1 rule of thumb there is amximum of 3 predictors to avoid overfitting.

The I use

M.avg <- model.avg(M.dredge, subset= delta < 4, fit=TRUE))

to get the model average across a subset of candidate models with a maximum difference of 4 AICc points to the "best" model having the lowest AICc. fit = TRUE refits every model in the subset.

Now I can use M.avg to get the averaged (full and conditional) coefficients for all predictors included in the subset as well as for the confdidence intervals:

Coef <- coef(M.avg, full=TRUE)
CI <- confint(M.avg, full=TRUE)

Now my question: Using poisson GLMMs and model averaging, how can I calculate model predictions and according confidence intervals for these predictions, to be used in a graph showing the regression line and confidence bands for that line?

Here they fill in the coefficients in the Model formula by hand and calculate point estimates for a given value (0, and 0.25) for the predictor f. Afterwards they repeat the procedure with the lower and the upper value of the confidence interval for the coefficients. However, I repeated that and I got something like this:

enter image description here

$\endgroup$
  • $\begingroup$ how about predict(se.fit = TRUE)? $\endgroup$ – Kamil Bartoń Apr 8 '16 at 8:16
  • $\begingroup$ This causes the following error: Error in predict.averaging(ls.avgModels[[1]], se.fit = TRUE) : 'predict' for models '21', '23', '29', '19', '17', '53', '31', '51', '49', '55', '61' and '63' caused errors In addition: There were 12 warnings (use warnings() to see them) $\endgroup$ – Pharcyde Apr 8 '16 at 12:08
  • $\begingroup$ And this is the Warning: 1: In predict.merMod(object = <S4 object of class structure("lmerMod", package = "lme4")>, ... : cannot calculate predictions with both standard errors and random effects $\endgroup$ – Pharcyde Apr 8 '16 at 12:09
  • 1
    $\begingroup$ ?predict.merMod says: re.form argument - "if NA or ~0, include no random effects." $\endgroup$ – Kamil Bartoń Apr 8 '16 at 12:50
  • $\begingroup$ Aaaaahhh Oh my god, I'm sorry I did not looked good enough in the help file. Thank you very much. $\endgroup$ – Pharcyde Apr 8 '16 at 14:47
2
$\begingroup$

Hi thanks to lots of hints by Kamil Bartón, i ended up with the following procedure:

M.global <- glmer(y ~ c1 + I(c1^2) + c2 + c3 + c4 + (1|subject),data, family =poisson(link="log))
M.dredge <- dredge(M.global, m.max=3, subset = dc(cpl,I(cpl^2)))
M.avg <- model.avg(M.dredge, subset= delta < 4, fit=TRUE))

newdata <- expand.grid(c1= with(data, seq(min(c1), max(c1), length.out=50)),
                       c2 = mean(data$c2),
                       c3 = mean(data$c3),
                       c4 = mean(data$c4))
# No random effects variable is needed!

# Calculate Predicitons with stnadard errors, and confidence intervals subsequently

pred.se <- predict(M.avg, type="response", se.fit=TRUE, re.form=NA, full=T, newdata)
  newdata$fit <- pred.se$fit
  newdata$SE <- pred.se$se
  newdata$upr=newdata$fit+1.96*newdata$SE
  newdata$lwr=newdata$fit-1.96*newdata$SE

# full=T defines the use of  full model-averaged coefficients see [here][1] 
# (full=F would not work either, due to mixed models); type= "response" 
# calculates predictions on the response scale (i.e. uses the inverse link
# function, exp() in poisson); re.form = NA avoids the use of random factors
# (inclusion is not possible with se.fit=TRUE).

Note: when using glmmadmb(family="nbinom2") the use of type="response" is not possible, but type="link" is. You could take the exp() of fit, SE, lwr, and upr, afterwards in the newdata frame.

The plots can be produced by:

p1 <- ggplot(newdata, aes(x=c1, y=fit)) +
    geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .25) +
    geom_line(size = 2) +
    labs(x = "c1 continous variable (restandardized)", y = "y response variable on original scale")
  print(p1)

Don't foget to re-standadize z transformed predictor variables in the final newdata frame. i.e.:

x <- attr(data, 'scaled:scale')
scale <- data.frame(t(x))
x <- attr(data, 'scaled:center')
center <- data.frame(t(x))

newdata$c1 <- newdata$c1 * scale$c1 + center$c1
newdata$c2 <- newdata$c2 * scale$c2 + center$c2
newdata$c3 <- newdata$c3 * scale$c3 + center$c3
newdata$c4 <- newdata$c4 * scale$c4 + center$c4

Now I am just wondering, why calculation of confidence intervals from single general linear mixed models (GLMMs) are so intensively discussed. Often you hear that you need complicated bootstrap methods or MarkovChainMonteCarlo (MCMC) algorithms to produce confidence intervals. Very, very complex for non-mathematicians. And even the function confint() takes much longer to calculate profile confidence intervals from a single mixed models than does from type=averaging model objects. Probably this is a Question of Accuracy and the se.fit = TRUE method is the least accurate? Furthermore i don't know if a correction method for Confidence intervals is included in the procedure presented here.

$\endgroup$

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.