[R-sig-ME] Incorrect output from nested model with mapped pars

[Ben Bolker]

    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