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]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models