Incorrect output from nested model with mapped pars

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 at r-project.org' <r-sig-mixed-models at 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]]

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 at r-project.org' <r-sig-mixed-models at 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]]