[R-sig-ME] Incorrect output from nested model with mapped pars
Ben Bolker
bbo|ker @end|ng |rom gm@||@com
Fri Jul 3 21:41:28 CEST 2020
Your first question looks like it could possibly be a bug, so please
post it on the glmmTMB github issues list (ideally with a reproducible
example!)
For your second question: using a t distribution implies knowing the
appropriate residual df, which is very difficult (see the GLMM FAQ or
any of the many documents floating around on the web for why the
sampling distributions of parameters from complex models are only
approximately t-distributed anyway, and why it is so hard to find good
approximations ...)
On 6/30/20 9:22 PM, Christopher Nottingham wrote:
> One more question. Why is the summary function giving the z value and associated p-value. A gaussian error structure is assumed, so shouldn't t values be used to obtain p values.
>
> Thanks,
> Chris
>
> From: Christopher Nottingham
> Sent: Tuesday, 30 June 2020 5:06 PM
> To: 'r-sig-mixed-models using r-project.org' <r-sig-mixed-models using r-project.org>
> Subject: Incorrect output from nested model with mapped pars
>
> I have a dataset with a variable labelled n_comm that is not relevant to some factor level combinations. I am fitting a nested model to this data and fixing betas representing the irrelevant factor combinations to 0 using map. As shown following, the model output from the summary table does not match what should be produced.
>
>> map_names = list(beta = factor(c(1:6, NA, 8)))
>> fit = glmmTMB(log(Err) ~ model + n_surv + species + n_comm:geostatistical + intensity,
> + data = Bhat_all.df,
> + start = list(beta = ifelse(is.na(map_names$beta), 0, 1)),
> + map = map_names)
>> summary(fit)
> Family: gaussian ( identity )
> Formula: log(Err) ~ model + n_surv + species + n_comm:geostatistical + intensity
> Data: Bhat_all.df
>
> AIC BIC logLik deviance df.resid
> 18340.8 18394.7 -9162.4 18324.8 6212
>
>
> Dispersion estimate for gaussian family (sigma^2): 1.11
>
> Conditional model:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) 9.329e+00 5.461e-02 170.82 <2e-16 ***
> modelbiomass-dynamics -4.082e+00 5.089e-02 -80.22 <2e-16 ***
> modelsize-structured -4.092e+00 5.056e-02 -80.93 <2e-16 ***
> n_surv -8.895e-04 8.146e-05 -10.92 <2e-16 ***
> species$\\mathit{S. aequilatera}$ 5.596e-01 2.677e-02 20.90 <2e-16 ***
> intensityFishing intensity: high 3.589e-01 2.681e-02 13.39 <2e-16 ***
> n_comm:geostatisticalFALSE 0.000e+00 7.450e-05 0.00 1.000
> n_comm:geostatisticalTRUE -3.245e-04 5.461e-02 -0.01 0.995
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> Warning message:
> In cbind(Estimate = coefs, `Std. Error` = sqrt(diag(vcov))) :
> number of rows of result is not a multiple of vector length (arg 2)
>> rbind(sqrt(diag(solve(fit$obj$he()))))
> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
> [1,] 0.05458617 0.05086506 0.05053495 8.142355e-05 0.02675588 0.02679859 7.446289e-05 0.01792267
>> fit$sdr
> sdreport(.) result
> Estimate Std. Error
> beta 9.3291879616 5.461347e-02
> beta -4.0822725635 5.089050e-02
> beta -4.0920631052 5.056022e-02
> beta -0.0008895195 8.146427e-05
> beta 0.5595933469 2.676926e-02
> beta 0.3589326271 2.681199e-02
> beta -0.0003244617 7.450013e-05
> betad 0.1082286532 1.793163e-02
> Maximum gradient component: 0.001913387
>
> The output below is wrong (there should be no std err, on a mapped value and the other values are incorrect.),
> n_comm:geostatisticalTRUE -3.245e-04 5.461e-02 -0.01 0.995
>
> The dataset is attached as a rds for reproducibility.
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using 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