[R-sig-ME] reporting model estimates values

[Juan Pablo Edwards Molina]

Thanks for your help prof. Bolker.
Best wishes,
Juan Edwards
Juan


On Tue, Oct 24, 2017 at 4:51 PM, Ben Bolker <bbolker at gmail.com> wrote:
>
>
>  (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