[R-sig-ME] reporting model estimates values

[Ben Bolker]

 (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