[R-sig-ME] reporting model estimates values
Fox, John
jfox at mcmaster.ca
Tue Oct 24 22:10:01 CEST 2017
Dear Ben and Juan,
Sorry to chime in late -- I was traveling when this thread started.
I'm not sure what Effect() method (the function in the effects package that does the computations) is getting invoked here -- probably the default method -- but in any event for a Poisson model it will use qnorm() not qt() to get confidence limits. I noticed that the effect estimates are identical (to the precision shown) but the standard errors of the effects are slightly different. I believe that's because Ben's computation of the SEs has removed the fixed-effect intercept, which is also subject to sampling variability and thus should figure in the SE of the "predicted" values.
Best,
John
-----------------------------
John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
Web: socserv.mcmaster.ca/jfox
> -----Original Message-----
> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org]
> On Behalf Of Ben Bolker
> Sent: Tuesday, October 24, 2017 2:52 PM
> To: Juan Pablo Edwards Molina <edwardsmolina at gmail.com>; r-sig-mixed-
> models at r-project.org
> Subject: Re: [R-sig-ME] 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="glmmTM
> >> B"))
> >> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list