reporting model estimates values

 (please keep r-sig-mixed-models in cc:)

  Without investigating too closely I'm not sure what drives the
difference.  effects could be using CIs based on t rather than Normal,
but with 641 residual df (df.residual(m1)) that wouldn't make that big a
difference.

  Ben Bolker


On 17-10-24 02:26 PM, Juan Pablo Edwards Molina wrote:

Thanks for your quick reply prof. Bolker and sorry for the misunderstanding.

 Yes, I meant the predicted values with standard errors for each factor level.

I tried both methods: predicted values and SE they present slight
differences in the CI values. I suppose it?s not a big deal...

# Extended method

pred

       [,1]
1 2.1363277
2 0.2219451

pse

[1] 0.2233638 0.1701403

exp(qnorm(c(0.025,0.975), mean=pred0[1],sd=pse[1]))

[1] 1.378924 3.309752

exp(qnorm(c(0.025,0.975), mean=pred0[2],sd=pse[2]))

[1] 0.1590091 0.3097914

# Experimental version of effect package:

as.data.frame(allEffects(m1)[[1]])

    mined       fit             se     lower     upper
1   yes 0.2219451 0.2230325 0.1433508 0.3436301
2    no 2.1363277 0.1705744 1.5292366 2.9844276

Thanks!

Juan


On Tue, Oct 24, 2017 at 1:12 PM, Ben Bolker <bbolker at gmail.com> wrote:

 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:

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

      [[alternative HTML version deleted]]