Skip to content

reporting model estimates values

2 messages · Juan Pablo Edwards Molina, Ben Bolker

Original
Juan Pablo Edwards Molina
# 
Dear glmmTMB authors & list members ,

For glmer models I use "lsmeans" or "effects" package to report the
coefficient estimated values, but any of them works with glmmTMB
objetcs:

lsmeans::lsmeans(mod, ~ factor, type = "response"); or

effects::allEffects(mod)

How would you suggest me to do that with an glmmTMB model?

Here is a short example (with only two factor levels):

library(glmmTMB)
m1 <- glmmTMB(count~ mined + (1|site),
               zi=~0,
               family=poisson, data=Salamanders)
summary(m1)

m2 <- glmer(count ~ mined + (1|site), family=poisson, data=Salamanders)
summary(m2)

library(lsmeans)
lsmeans(m2, ~ mined, type = "response")

library(effects)
eff.m2 <- allEffects(m2); as.data.frame(eff.m2[[1]])

lsmeans(m1, ~ mined, type = "response")
eff.m2 <- allEffects(m1); as.data.frame(eff.m1[[1]])

Thanks in advance,

Juan Edwards
Ben Bolker
# 
Depends what you mean by estimated coefficient values.  You already
have the summaries; I assume you mean you want predicted values with
standard errors for each factor level?

  You could do it like this: (this brute-force method will work for any
model where you can extract the formula, fixed-effect parameters,
variance-covariance matrix of fixed-effect parameters)

## set up prediction frame
pframe <- with(Salamanders,
               data.frame(mined=levels(mined)))
## get formula and keep only fixed effects
fform <- lme4::nobars(formula(m1))
## drop LHS from formula
fform_noresp <- fform[-2]
## model matrix for new predictions
X <- model.matrix(fform_noresp,data=pframe)
## fixed-effect parameters
beta <- fixef(m1)[["cond"]]
## predictions (linear predictor scale)
pred0 <- X %*% beta
## predictions (back-transformed)
pred <- exp(pred0)
## standard errors of predictions (log scale)
pse <- sqrt(diag(X %*% crossprod(X,vcov(m1)[["cond"]])))
## back-transformed confidence intervals
exp(qnorm(c(0.025,0.975),mean=pred0[1],sd=pse[1]))
exp(qnorm(c(0.025,0.975),mean=pred0[2],sd=pse[2]))

Alternatively, you can use an experimental version that has the glue
needed for the effects package.

devtools::install_github("glmmTMB/glmmTMB/glmmTMB",ref="effects")
library(glmmTMB)
source(system.file("other_methods","effectsglmmTMB.R",package="glmmTMB"))
allEffects(m1)


  Caveats:

   - this will only work at present for families known to base-R
(poisson, gaussian, binomial, Gamma, etc.) - not for nbinom2 etc.
   - it makes lots of assumptions. In particular, for zero-inflated
models it ignores the zero-inflation part completely and gives
predictions etc etc only for the conditional model
On 17-10-24 09:16 AM, Juan Pablo Edwards Molina wrote: