`glmmTMB`: Standard Errors and `REML=TRUE`

[strengejacke]

I noticed this when writing a vignette for ggeffects. The standard errors for glmmTMB predictions of random effects are smaller than what predict() returns. Is this intended or a bug?

Reprex requires datawizard dev version...

library(ggeffects)   # predictions and significance testing
library(insight)     # extracting random effects variances
library(datawizard)  # data wrangling and preparation
library(glmmTMB)     # multilevel modelling

# sample data set
data(efc, package = "ggeffects")

efc <- efc |>
  # numeric to factors, set labels as levels
  to_factor(select = c("c161sex", "c172code", "c175empl")) |>
  # recode age into three groups
  recode_values(
    select = "c160age",
    recode = list(`1` = "min:40", `2` = 41:64, `3` = "65:max")
  ) |>
  # rename variables
  data_rename(
    pattern = c("c161sex", "c160age", "quol_5", "c175empl"),
    replacement = c("gender", "age", "qol", "employed")
  ) |>
  # age into factor, set levels, and change labels for education
  data_modify(age = factor(age, labels = c("-40", "41-64", "65+")))

m_null <- glmmTMB(qol ~ 1 + (1 | gender:employed:age), data = efc)

do.call(cbind, predict(
  m_null,
  newdata = data_grid(m_null, c("gender", "employed", "age")),
  se.fit = TRUE
))
#>                  fit    se.fit
#> eta_predict 15.24091 0.8758925
#> eta_predict 15.17225 0.6872397
#> eta_predict 16.12010 0.7850133
#> eta_predict 16.03567 0.6105341
#> eta_predict 14.93677 0.6589485
#> eta_predict 13.75407 0.3305532
#> eta_predict 15.46836 0.5570139
#> eta_predict 14.09512 0.3461792
#> eta_predict 14.41308 0.6287481
#> eta_predict 13.53274 0.4361047
#> eta_predict 15.34190 1.0851461
#> eta_predict 14.85284 1.0382235

ggpredict(
  m_null,
  c("gender", "employed", "age"),
  type = "random",
  interval = "confidence"
) |> as.data.frame()
#>         x predicted std.error conf.low conf.high group facet
#> 1    Male  15.24091 0.8758925 13.52419  16.95763    no   -40
#> 2    Male  14.93677 0.6589485 13.64525  16.22829    no 41-64
#> 3    Male  14.41308 0.6287481 13.18075  15.64540    no   65+
#> 4    Male  16.12010 0.7850133 14.58150  17.65870   yes   -40
#> 5    Male  15.46836 0.5570139 14.37663  16.56008   yes 41-64
#> 6    Male  15.34190 1.0851461 13.21505  17.46875   yes   65+
#> 7  Female  15.17225 0.6872397 13.82528  16.51921    no   -40
#> 8  Female  13.75407 0.3305532 13.10620  14.40195    no 41-64
#> 9  Female  13.53274 0.4361047 12.67799  14.38749    no   65+
#> 10 Female  16.03567 0.6105341 14.83904  17.23229   yes   -40
#> 11 Female  14.09512 0.3461792 13.41662  14.77362   yes 41-64
#> 12 Female  14.85284 1.0382235 12.81795  16.88772   yes   65+

marginaleffects::predictions(
  m_null,
  newdata = data_grid(m_null, c("gender", "employed", "age"))
)
#> 
#>  Estimate Std. Error     z Pr(>|z|) 2.5 % 97.5 % gender employed   age
#>      15.2     0.2519  60.5   <0.001  14.7   15.7 Male        no  -40  
#>      15.2     0.1633  92.9   <0.001  14.9   15.5 Female      no  -40  
#>      16.1     0.1737  92.8   <0.001  15.8   16.5 Male        yes -40  
#>      16.0     0.1167 137.5   <0.001  15.8   16.3 Female      yes -40  
#>      14.9     0.1533  97.4   <0.001  14.6   15.2 Male        no  41-64
#>      13.8     0.0391 351.4   <0.001  13.7   13.8 Female      no  41-64
#>      15.5     0.1072 144.2   <0.001  15.3   15.7 Male        yes 41-64
#>      14.1     0.0435 324.3   <0.001  14.0   14.2 Female      yes 41-64
#>      14.4     0.1388 103.8   <0.001  14.1   14.7 Male        no  65+  
#>      13.5     0.0646 209.4   <0.001  13.4   13.7 Female      no  65+  
#>      15.3     0.3578  42.9   <0.001  14.6   16.0 Male        yes 65+  
#>      14.9     0.3568  41.6   <0.001  14.2   15.6 Female      yes 65+  
#> 
#> Columns: rowid, estimate, std.error, statistic, p.value, conf.low, conf.high, gender, employed, age, qol

^{Created on 2023-06-08 with reprex v2.0.2}

[vincentarelbundock]

Thanks for the report. I’m travelling right now so I can’t investigate, but one potential culprit is that get_parameters() and get_varcov() are returning objects of different dimensions:

insight::get_parameters(m_null, components = "all")
#     Parameter  Estimate   Component
# 1 (Intercept) 14.913652 conditional
# 2 (Intercept)  3.308191  dispersion

insight::get_varcov(m_null, components = "all")
#             (Intercept)
# (Intercept)   0.1589324

[strengejacke]

There's a typo, the argument is named component. But still, get_varcov() also returns the parameter for the random effects. so the mismatch between dimensions remains:

m_null <- glmmTMB(qol ~ 1 + (1 | gender:employed:age), data = efc)

insight::get_parameters(m_null, component = "all")
#>     Parameter  Estimate   Component
#> 1 (Intercept) 14.913652 conditional
#> 2 (Intercept)  3.308191  dispersion
insight::get_varcov(m_null, component = "all")
#>                                 (Intercept) d~(Intercept)
#> (Intercept)                    0.1589323969 -0.0004477001
#> d~(Intercept)                 -0.0004477001  0.0022586882
#> theta_1|gender:employed:age.1  0.0280608157 -0.0003900293
#>                               theta_1|gender:employed:age.1
#> (Intercept)                                    0.0280608157
#> d~(Intercept)                                 -0.0003900293
#> theta_1|gender:employed:age.1                  0.0950342589

^{Created on 2023-06-08 with reprex v2.0.2}

I guess d~(Intercept) is the vcov for the dispersion parameter. In insight, we should ensure that each value for component refers to the same parameter values for the different functions (though there already has been some effort to harmonize this), and sometimes, "all" is not desired, but rather something like "everything but random effects" or similar.

[vincentarelbundock]

@strengejacke, this is really bothering me and I can’t put my finger on the problem. There are many weird things going on. For example, the glmmTMB and lme4 SEs don’t seem compatible at all:

library(glmmTMB)
library(lme4)
library(insight)
library(marginaleffects)

m1 <- glmmTMB(Sepal.Length ~ Sepal.Width + (1 | Species), data = iris)
m2 <- lmer(Sepal.Length ~ Sepal.Width + (1 | Species), data = iris)

# glmmTMB
insight::get_predicted(m1, ci = .95) |> data.frame() |> head(2)
#   Predicted         SE   CI_low  CI_high
# 1  5.069506 0.06228249 4.947434 5.191577
# 2  4.672506 0.07681069 4.521960 4.823053
predict(m1, se.fit = .95) |> data.frame() |> head(2)
#        fit     se.fit
# 1 5.069506 0.06228249
# 2 4.672506 0.07681069

# lme4
insight::get_predicted(m2, ci = .95) |> data.frame() |> head(2)
#   Predicted        SE   CI_low  CI_high
# 1  5.067641 0.5860123 3.909478 6.225803
# 2  4.669063 0.5841551 3.514571 5.823556
marginaleffects::predictions(m2) |> head(2)
# 
#  Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
#      5.07      0.586 8.65   <0.001 57.4  3.92   6.22
#      4.67      0.584 7.99   <0.001 49.4  3.52   5.81
# 
# Columns: rowid, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, Sepal.Length, Sepal.Width, Species

[vincentarelbundock]

I opened a related issue here: https://github.com/glmmTMB/glmmTMB/issues/915

[saudiwin]

Totally not a problem @vincentarelbundock . In my next version, I'm likely going to stress using stan's optimizer vs. glmmTMB as it turns out ordbeta can be estimated fairly easily with MLE (as far as I can tell). Removing a package dependency is probably good anyways as CRAN faults on every vignette build. As long as you keep supporting ordbetareg/brms I'm happy! :)

[strengejacke]

Removing SE computation completely for glmmTMB was probably a too radical step ;-) It seems to work for non-Gaussian models - maybe only Gaussian is affected? See my comment here: https://github.com/glmmTMB/glmmTMB/issues/915#issuecomment-1600252192

[vincentarelbundock]

Yes, you're probably right. I've recently been freaking out about the possibility that people publish incorrect results because of me, so I prefer to take drastic conservative steps when I see bad stuff is possible.

Also, even for Gaussian, some models are very similar between glmmTMB and lme4. So I don't find single case comparisons sufficiently reassuring...

[saudiwin]

@vincentarelbundock I might consider adding a standard text in the marginaleffects function that says something to the effect of "note: calculation of standard errors assumes model fit produces asymptotically correct standard errors" or something to that effect.

[vincentarelbundock]

@vincentarelbundock I might consider adding a standard text in the marginaleffects function that says something to the effect of "note: calculation of standard errors assumes model fit produces asymptotically correct standard errors" or something to that effect.

Isn't that implicit in "we use the delta method"?

[huftis]

Any chance you could whitelist the SE calculations (i.e., allow vcov = TRUE) for the models where there aren’t any problems? This seems to be true for at least OLS models and GLM models (i.e., models without any random effects). Example:

library(glmmTMB)
library(insight)

diff_vcov = function(m1, m2) {
  v1 = get_varcov(m1, components = "all")
  v2 = get_varcov(m2, components = "all")
  max(abs(v1 - v2))
}

# Linear model
m1 = lm(Sepal.Width ~ poly(Sepal.Length, 3) +
  Species * Petal.Length, data = iris)
m2 = glmmTMB(formula(m1), data = iris, REML = TRUE)
diff_vcov(m1, m2)
#> [1] 7.36973e-10

# Logistic model (GLM)
dat = tidyr::uncount(as.data.frame(Titanic), Freq)
m1 = glm(Survived ~ Class * Sex + Age,
  data = dat, family = "binomial")
m2 = glmmTMB(formula(m1), data = dat, family = "binomial")
diff_vcov(m1, m2)
#> [1] 2.199957e-06

# Poisson model (GLM)
m1 = glm(breaks ~ wool * tension,
  data = warpbreaks, family = "poisson")
m2 = glmmTMB(formula(m1),
  data = warpbreaks, family = "poisson")
diff_vcov(m1, m2)
#> [1] 2.991652e-08

[vincentarelbundock]

I appreciate where you're coming from, but I don't think I'll do that because:

1. I have to admit that I still don't understand exactly where the problem is coming from, and I prefer not to do guesswork about which models are problematic or not.
2. I prefer not to write a bunch of package specific code to cover weird exceptions.
3. It seems trivial for users to add vcov = vcov(model) to the call.

So there are risks for users and work for me involved, and I don't think the benefits are that great.

Hopefully I'll get a chance to dive in to really understand what's happening and eventually come up with a better solution.

[huftis]

I don’t fully understand where the SE problem is coming from either. But based https://github.com/glmmTMB/glmmTMB/issues/915, it seems to be related to the random effects, and on whether we take into account the uncertainty in the random effects (including the correlations of the random effect estimators). Perhaps both SEs are technically correct, they’re just taking into account different things?

I certainly understand your point no. 2. 😃

Regarding point 3, the problem I’m having is that I’m often fitting various similar models (with either different predictors, different outcomes or on different subsets), like this:

# Original model
mod1 = glmmTMB(something ~ bunch_of_predictors, data = dat)

mod2 = glmmTMB(something_else ~ ., data = dat)
# or
mod3 = update(mod1, . ~ . + extra_predictors)
# or
mod4 = update(mod1, data = dat_subset)

And then I’m doing something like:

predictions(mod1, newdata = ..., other_options, vcov = vcov(mod1))
# or
plot_predictions(mod1, newdata = ..., other_options, vcov = vcov(mod1))

And then I’m partly reusing this code for the other models, so I end up with:

predictions(mod2, newdata = ..., other_options, vcov = vcov(mod1))

And I get completely wrong SEs / confidence intervals, since I have forgotten to also update the model object in the vcov argument (note: I don’t get an error message, so this is easy to miss). Of course, this is all my own fault. But it’s been happening more times than I’m willing to admit. 😉

So it would be nice to at least have an option like vcov = "yes-i-really-do-trust-the-vcov-from-my-model", so we don’t have to remember to put in the correct model in the vcov argument all the. (A more natural name might be vcov = TRUE, with the default instead being changed to vcov = NULL, taking on the current default behaviour.)

[vincentarelbundock]

I see what you mean and understand that my "no" is a bit inconvenient.

Perhaps you should consider using fitting and updating functions, instead of copying the code, that way, the difference between mod1 and mod2 would be irrelevant.

[saudiwin]

So my take on this is that @vincentarelbundock is right, we have to be careful of standard errors and it's difficult to alert the end user to these issues. What needs to happen in my opinion is to implement an MLE version of ordbetareg in the package. Right now I can't do that because brms doesn't support optimization, but I'm planning on adding that support before the end of the year.

[bwiernik]

It looks like it might be the lme4 SEs that are problematic rather than glmmTMB?

https://github.com/glmmTMB/glmmTMB/issues/915

[strengejacke]

It looks like it might be the lme4 SEs that are problematic rather than glmmTMB?

glmmTMB/glmmTMB#915

I agree, and maybe calculating SEs for glmmTMB by default should be added back, while a warning should be given for merMod?

[vincentarelbundock]

@bbolker, I would really appreciate your advice on what to do here. In the other thread, you write:

Yes. I’m a little bit puzzled that the solution to this problem was to disable calculations based on glmmTMB. Arguably glmmTMB’s internal SE calculation is the most accurate method we have:

The reason I disabled calculations based on glmmTMB in marginaleffects is that marginaleffects does not use the internal TMB estimates from predict(se.fit=TRUE). We are not interested only in SEs for predictions, but also for a bunch of other transformations of the parameters: differences, odds, slopes, etc. To get them, we just use a simple delta method approach:

• Construct an NxP matrix called J where:
  • N: number of estimates, ex: 1 prediction or slope per observation
  • P: number of parameters in coef(model)
  • Entries in column 2 are partial derivatives of the quantity (prediction, slope, etc.) w.r.t. to the 2nd parameter.
• Use this jacobian and vcov(model) by standard multiplication to compute standard errors for whatever quantity the user requested.

When I do this, results are very different from what you get internally from TMB.

1. Is there another strategy you would recommend?
2. If not, what caveats should I include in warnings or documentation?

library(glmmTMB)
library(marginaleffects)

m <- glmmTMB(Petal.Width ~ Sepal.Length + (1 | Species), iris)

predict(m, se.fit = TRUE) |> data.frame() |> head()
#         fit     se.fit
# 1 0.2613778 0.02712136
# 2 0.2318260 0.02715171
# 3 0.2022742 0.02852497
# 4 0.1874983 0.02966407
# 5 0.2466019 0.02696368
# 6 0.3057055 0.02955279

predictions(m, vcov = vcov(m)) |> head()
# 
#  Estimate Std. Error     z Pr(>|z|)     S 2.5 % 97.5 %
#     0.261    0.00324  80.7   <0.001   Inf 0.255  0.268
#     0.232    0.00348  66.5   <0.001   Inf 0.225  0.239
#     0.202    0.00941  21.5   <0.001 337.9 0.184  0.221
#     0.187    0.01244  15.1   <0.001 168.0 0.163  0.212
#     0.247    0.00140 176.2   <0.001   Inf 0.244  0.249
#     0.306    0.01218  25.1   <0.001 459.6 0.282  0.330
# 
# Columns: rowid, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, Petal.Width, Sepal.Length, Species 
# Type:  response

[bbolker]

Internally, glmmTMB/TMB also use a delta-method approach (although in this case the predictions are a linear combination of fixed & random effects, so the 'delta method' is exact/just linear variance propagation).

Here's an equivalent form that you could use externally: the issue is whether we're including just the variance due to the fixed effects, or whether we're also including the variances of the conditional modes/BLUPs ("random effects") and their covariances with the fixed effects. You can get the full joint covariance matrix, as shown below, but it will be expensive for models with high dimensions (= large numbers of fixed and random effects/levels of the grouping variables), especially the solve() step. It may be possible to do a little better on the linear algebra, and we can be more efficient in the last step (see lme4:::quad.tdiag, adapted from the emulator package).

pp <- predict(m, se.fit = TRUE) |> data.frame()

s <- TMB::sdreport(m$obj, getJointPrecision = TRUE)
Sigma <- solve(s$jointPrecision)
keepvars <- rownames(Sigma) %in% c("beta", "b")
Sigma <- Sigma[keepvars, keepvars]
XZ <- as.matrix(cbind(getME(m, "X"), getME(m, "Z")))
se <- sqrt(diag(XZ %*% (Sigma %*% t(XZ))))

all.equal(unname(se), pp$se.fit)  ## TRUE

Without digging into your code, I'm not sure what it's doing with vcov(). When I pass Sigma (i.e. the joint covariance matrix of beta and b) to predictions instead of vcov(m), it doesn't complain, but it also doesn't change the results. (The mean prediction obviously is using the random effects ...)

You can warn users that "standard errors take only uncertainty in the fixed effects into account", but I don't know how much that will mean to them; I didn't appreciate the potential magnitude of the issue until the latest round of digging into the guts of the lme4 machinery ... (worth noting, though, that the example you've constructed here is extreme - good for pointing out differences but the differences here may be atypically large)

[vincentarelbundock]

I have another question, @bbolker, but would of course understand if you don’t have time.

Reproducing the predict.glmmTMB result manually

# fit
library(glmmTMB)
dat <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/causaldata/thornton_hiv.csv") |> na.omit()
m <- glmmTMB(got ~ age * hiv2004, data = dat)

# extract
J <- getME(m, "X")
V <- vcov(m)$cond

# standard errors
manu <- sqrt(diag(J %*% (V %*% t(J))))
auto <- predict(m, fast = FALSE, se.fit = TRUE)$se.fit
all.equal(unname(manu), unname(auto))
# [1] TRUE

The same J matrix can be obtained by numeric differentiation

In principle, we should be able to construct an equivalent J by numerically differentiating the predictions with respect to each of the coefficients. This is obviously much slower than getting the model matrix, but it is a more general approach for all the quantities in marginaleffects. In this example, I modify the second parameter by a tiny amount inside the object to make different predictions, and then do forward finite difference:

m_alt <- m
m_alt$fit$parfull[2] <- m_alt$fit$parfull[2] + 1e-6
m_alt$fit$par[2] <- m_alt$fit$par[2] + 1e-6

# forward finite difference of predictions w.r.t. beta[2]
J2 <- (predict(m_alt, fast = FALSE) - predict(m, fast = FALSE)) / 1e-6

all.equal(unname(J[, 2]), unname(J2))
# [1] TRUE

The same strategy fails when there are random effects

Unfortunately, the same strategy does not produce a vector equivalent to the one produced by the code you gave me.

Here is the J matrix you use:

m <- glmmTMB(got ~ age * hiv2004 + (1 | villnum), data = dat)

J <- as.matrix(cbind(getME(m, "X"), getME(m, "Z")))

Now try to get the same with the strategy from the previous section:

m_alt <- m
m_alt$fit$parfull[2] <- m_alt$fit$parfull[2] + 1e-6
m_alt$fit$par[2] <- m_alt$fit$par[2] + 1e-6

J2 <- (predict(m_alt, fast = FALSE) - predict(m, fast = FALSE)) / 1e-6

all.equal(unname(J[, 2]), unname(J2))
# [1] "Mean relative difference: 0.7119016"

Since J is different, my standard errors are also different.

Naive question: is there an equivalent way to produce your XZ matrix by modifying the glmmTMB fit object, and calling predict()?

[bbolker]

I'll see if I can take a look.

[bbolker]

I think this works: it's a modified version of the predict(., fast=TRUE) code, with the important addition that it sets the map element to fix the b values (and prevent them from being modified in the inner loop of the Laplace approximation). I suppose it would be pretty easy to add a fix_b argument? (In fact I think this could/should be done in conjunction with the newparams argument ...)

This code is fairly general, but doesn't necessarily handle all of the possible use cases:

• zero-inflation stuff
• predictions on scales other than link/lin pred scale? Standard errors?
• does it properly restore all of the stored versions of parameter vectors?

m <- glmmTMB(got ~ age * hiv2004 + (1 | villnum), data = dat)
J <- as.matrix(cbind(getME(m, "X"), getME(m, "Z")))
delta <- 1e-6
J_ind <- 2
## predict machinery by hand ('fast' version)
do_pred_val <- 0

ee <- m$obj$env
dd <- ee$data         ## data object

## store copies of anything we
## are modifying so we can restore the environment at the end ...

orig_vals <- dd[c("whichPredict", "doPredict", "ziPredictCode")]
orig_map <- ee$map
orig_lp <- ee$last.par.best
## used in $report() call below
lp <- orig_lp
lp[J_ind] <- lp[J_ind] + delta

dd$ziPredictCode <- 0
dd$whichPredict <- as.numeric(seq(nobs(m)))  ## replace 'whichPredict'
dd$doPredict <- 0  ## no SEs
assign("data",dd, ee) ## stick this in the appropriate environment
b_map <- list(b = factor(rep(NA_integer_, length(ee$parList()[["b"]]))))
map <- orig_map
if (is.null(map)) {
    map <- b_map
} else {
    if (is.null(map[["b"]])) {
        map <- c(map, b_map)
    } else {
        map[["b"]] <- b_map[[1]]
    }
}

pred_dev <- m$obj$report(lp)$mu_predict
## restore original environment
for (i in names(orig_vals)) {
    dd[[i]] <- orig_vals[[i]]
    assign("data",dd, ee)
}
ee$map <- orig_map
pred0 <- predict(m)

J2 <- (pred_dev - pred0) / delta
all.equal(unname(J[, J_ind]), unname(J2))

[vincentarelbundock]

Ah, that's really great! Thanks so much.

Just to be sure: Would this also work with fast=FALSE? I need to use the newdata argument...

Is there anything I can do to help with this?

[bbolker]

Some version of it should work. The key is setting the map element, the rest is just plumbing. To be clear, you'd like to be able to call predict with newparams (that's effectively what you're doing when you perturb the parameters), and newdata, and have b held fixed ... ?

[vincentarelbundock]

Yes, that's exactly right. It would be a dream if all developers implemented a newparams argument. Alas...

[vincentarelbundock]

I am just pasting personal notes to close the duplicate thread. This is not new info.

library(marginaleffects)
library(glmmTMB)
dat <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/causaldata/thornton_hiv.csv")
m <- glmmTMB(got ~ age * hiv2004 + (1 | villnum), data = dat)

get_se_manual <- function(m, random = FALSE) {
    s <- TMB::sdreport(m$obj, getJointPrecision = TRUE)
    Sigma <- solve(s$jointPrecision)
    if (isTRUE(random)) {
        XZ <- as.matrix(cbind(getME(m, "X"), getME(m, "Z")))
        keepvars <- c("beta", "b")
    } else {
        XZ <- as.matrix(getME(m, "X"))
        keepvars <- "beta"
    }
    keepvars <- rownames(Sigma) %in% keepvars
    Sigma <- Sigma[keepvars, keepvars]
    se <- sqrt(diag(XZ %*% (Sigma %*% t(XZ))))
    return(se)
}

se1 <- get_se_manual(m, random = TRUE)
se2 <- predict(m, se.fit = TRUE)$se.fit

V <- vcov(m)$cond
J <- model.matrix(m)
se3 <- sqrt(diag(M %*% (V %*% t(M))))
se4 <- get_se_manual(m, random = FALSE)

all.equal(unname(se1), unname(se2))
all.equal(unname(se3), unname(se4))

# mismatch
se5 <- predictions(m, vcov=vcov(m))$std.error

[DrJerryTAO]

@vincentarelbundock, though consistent in the latest personal notes and previous scripts, for se3 <- sqrt(diag(M %*% (V %*% t(M)))), the parenthesis around V %*% t(M) makes no difference. I reviewed the associativity of matrix multiplication to confirm it, because I was initially worried about my own codes that did not apply the parentheses. I used Jacobian from marginaleffects and vcov() of model coefficients to recover the variance-covariance matrix of predictions and comparisons. I suggest adding the variance-covariance matrix of predictions and comparisons as output of predictions() and comparisons(), as further derivation (such as group comparisons) among these output require beyond solely standard errors but also correlation between groups.

The same strategy fails when there are random effects

The reason is that that the random effects should not be taken differential, as it is just a group indicator. It won't make sense to find response changes if the dummy variable, equal to 1 for observations in the first group, changes by a tiny fraction.

library(glmmTMB)
summary(Model_Temp[[1]] <- glmmTMB(
  got ~ age * hiv2004  + (1 | villnum), data = read.csv(paste0(
  "https://vincentarelbundock.github.io/Rdatasets/csv/", 
  "causaldata/thornton_hiv.csv")) |> 
    na.omit()))
"     AIC      BIC   logLik deviance df.resid 
  3578.8   3614.4  -1783.4   3566.8     2819 
Random effects:
 Groups   Name        Variance Std.Dev.
 villnum  (Intercept) 0.01648  0.1284  
 Residual             0.19878  0.4459  
Number of obs: 2825, groups:  villnum, 119
Dispersion estimate for gaussian family (sigma^2): 0.199 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.6317232  0.0262568  24.059  < 2e-16 ***
age          0.0022487  0.0006390   3.519 0.000433 ***
hiv2004     -0.0600519  0.1236294  -0.486 0.627151    
age:hiv2004  0.0005985  0.0032872   0.182 0.855530    "
getME(Model_Temp[[1]], "X") # matrix
getME(Model_Temp[[1]], "Z") # 'dgTMatrix'
cbind(getME(Model_Temp[[1]], "X"), as.matrix(getME(Model_Temp[[1]], "Z"))) %>%
  as.data.frame()
"   (Intercept) age hiv2004 age:hiv2004 1 2 3 5 6 7 8 9 10 11 12 13 14 15 16 17
1            1  22       0           0 1 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0
3            1  19       0           0 1 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0
5            1  53       0           0 1 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0
6            1  50       0           0 1 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0
11           1  21       0           0 1 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0
12           1  62       0           0 1 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0
15           1  47       1          47 1 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0
17           1  15       0           0 1 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0..."

The logarithms of random effect standard deviation and residual variance are represented by theta and betad in glmmTMB()$fit$par objects. However, only one of them needs estimation, as the other will be determined given one of them. Therefore, insight::get_parameters() reports one of them (4 fixed effect + 1 log residual var, 5 parameters) whereas insight::get_varcov() gives the covariance matrix of all (6 by 6 matrix). Furthermore, I realized that estimates of random intercepts, as numeric values named b in glmmTMB()$fit$parfull of length as the count of (1 | villnum) groups cannot recover the random-effect standard deviation as reported in the model summary. Can someone tell why it is so?

insight::get_parameters(Model_Temp[[1]], component = "all")
"    Parameter      Estimate   Component
1 (Intercept)  0.6317231594 conditional
2         age  0.0022487446 conditional
3     hiv2004 -0.0600518660 conditional
4 age:hiv2004  0.0005984992 conditional
5 (Intercept) -1.6155417453  dispersion
Here exp(-1.6155417453) = 0.198783 is the residual variance
Random-effect parameters are not included here"
Model_Temp[[1]]$fit$par["betad"] # -1.6155417453
exp(Model_Temp[[1]]$fit$par["betad"]) # 0.198783 = residual variance
Model_Temp[[1]]$fit$par["theta"] # -2.052738
exp(Model_Temp[[1]]$fit$par["theta"]) # 0.1283829 = random effect SD
Model_Temp[[1]]$fit$parfull %>%
  `[`(names(.) == "b") %>% sd() # 0.09753274, not 0.1284 as model print
vcov(Model_Temp[[1]]) # only gives fixed-effect part
"              (Intercept)           age       hiv2004   age:hiv2004
(Intercept)  6.894172e-04 -1.342988e-05 -4.046168e-04  1.036128e-05
age         -1.342988e-05  4.082644e-07  1.064802e-05 -3.217706e-07
hiv2004     -4.046168e-04  1.064802e-05  1.528423e-02 -3.908703e-04
age:hiv2004  1.036128e-05 -3.217706e-07 -3.908703e-04  1.080591e-05"
insight::get_varcov(Model_Temp[[1]], component = "all")
"                    (Intercept)           age       hiv2004   age:hiv2004
(Intercept)        6.894172e-04 -1.342988e-05 -4.046168e-04  1.036128e-05
age               -1.342988e-05  4.082644e-07  1.064802e-05 -3.217706e-07
hiv2004           -4.046168e-04  1.064802e-05  1.528423e-02 -3.908703e-04
age:hiv2004        1.036128e-05 -3.217706e-07 -3.908703e-04  1.080591e-05
d~(Intercept)     -6.302971e-06 -6.070752e-08  4.966683e-06 -1.728403e-07
theta_1|villnum.1  1.275741e-04  1.228732e-06 -1.005286e-04  3.498384e-06
                  d~(Intercept) theta_1|villnum.1
(Intercept)       -6.302971e-06      1.275741e-04
age               -6.070752e-08      1.228732e-06
hiv2004            4.966683e-06     -1.005286e-04
age:hiv2004       -1.728403e-07      3.498384e-06
d~(Intercept)      7.398080e-04     -2.905456e-04
theta_1|villnum.1 -2.905456e-04      1.322257e-02
Here sqrt(7.398080e-04) = 0.02719941 about d~(Intercept) should be the  
standard error of betad = log(residual variance), while 
sqrt(1.322257e-02) = 0.1149894 about theta_1|villnum.1 should be the 
standard error of theta = log(random effect SD), 
although neither is confirmed"

You can warn users that "standard errors take only uncertainty in the fixed effects into account"

I agree with @bbolker on this. This will be consistent with how marginaleffects treat models with random effects such as lmer() models. Unlike ggeffects that allows confidence intervals to include variance component of the random effects, marginaleffects focuses on effects on the population mean unless you intend to expand the functionality. To incorporate both parameter uncertainty and random effects, users can simply add variances from both sources and take the square root to find the standard error, because models with random effects have to assume that the random effects are uncorrelated with any fixed-effect predictors. If the individual-specific effects are correlated with other fixed-effect predictors, fixed-effect estimators need to be used by simply including dummy variables of group membership.

In the hiv example, we can find that the squared standard errors from predict.glmmTMB() is larger than those from predictions() by a positive constant. However, I could not find the relationship between this constant and any model estimates of glmmTMB(). Could someone interpret what this difference in squared standard errors is?

data.frame(predict(Model_Temp[[1]], se.fit = TRUE)) %>% head()
"          fit     se.fit
1   0.6246165 0.07009961
2   0.6178703 0.07029221
3   0.6943276 0.07116783
4   0.6875814 0.07082346
5   0.6223678 0.07015805
6   0.7145663 0.07249584"
predictions(Model_Temp[[1]], vcov = vcov(Model_Temp[[2]])) %>% head()
" Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
    0.625    0.00853 73.2   <0.001   Inf 0.608  0.641
    0.618    0.00999 61.9   <0.001   Inf 0.598  0.637
    0.694    0.01495 46.4   <0.001   Inf 0.665  0.724
    0.688    0.01322 52.0   <0.001   Inf 0.662  0.713
    0.622    0.00900 69.2   <0.001   Inf 0.605  0.640
    0.715    0.02036 35.1   <0.001 894.2 0.675  0.754"
(predict(Model_Temp[[1]], se.fit = TRUE)$se.fit^2 - 
   predictions(Model_Temp[[1]], vcov = vcov(Model_Temp[[2]]))$std.error^2) %>%
  head()
"0.004841222 0.004841222 0.004841222 0.004841222 0.004841222 0.004841222
Its square root is 0.06957889"

[vincentarelbundock]

Notes to self:

https://github.com/vincentarelbundock/marginaleffects/issues/970#issuecomment-1840857440

https://github.com/vincentarelbundock/marginaleffects/pull/1023

[DrJerryTAO]

This is an update of current function behavior and fix suggestions regarding glmmTMB::predict() and marginaleffects::predictions().

• predict(Model, re.form = NA, se.fit = T) matches predictions(Model, re.form = NA, vcov = vcov(Model)) exactly. This is to omit any random effects (individual-level predictions) and seek only population-level predictions.
• predict(Model, re.form = NULL, se.fit = T) matches predictions(Model, re.form = NULL, vcov = vcov(Model)) in only point estimates, where re.form = NULL is the default setting to request individual-level predictions with random effects, but the latter reports a smaller standard error due to not incorporating uncertainty in BLUP estimates of random effects of individual groups.
• The current suggestion to "explicitly provide a variance-covariance-matrix for the vcov argument to calculate standard errors" is incomplete. It must also include re.form = NA or ~ 0 to remove random effects and states that it includes uncertainty of fixed effects only. Supplying a larger vcov such as from insight::get_varcov(Model, component = "all") that does include covariance with random effects will not work because only its upper-left square of fixed effects is used.
• To match marginaleffects results with predict() when re.form = NULL, we need to (1) redefine get_coef() and set_coef() to extract not Model$fit$par but Model$fit$parfull that includes BLUP random coefficients of all groups, (2) replace get_vcov() just like get_se_manual() defined in @vincentarelbundock 's last notes but with re.form = NULL as default argument and returns potentially huge-dimension matrix Sigma; (3) redefine Jacobian, currently it only calculates for model parameters including those for betad = log(residual variance) and theta = log(random effect SD) but not including those for b = BLUP random coefficients. We need to remove those for betad, theta and add those for b. I tried with updated get_coef(), set_coef(), get_vcov() but still got previous results in predictions. I think it is because Jacobian and vcov is still on model parameters, neither of which includes those for b = BLUP.
• Once the above fixes are done, there might be still challenges in making adjusted predictions and comparisons. Among the expanded model matrix including random effects by group, only the columns getME(m, "X") corresponding to fixed effects should be adjusted (e.g. from 100 to 101, or from level B to level A). Columns like getME(m, "Z") are simply group indicators and predictors of random effects, which I don't think should be adjusted for new data. Further, estimates of random effects are conditional on fixed effects. When fixed effects are all together adjusted for new data (e.g. every person's age increases by one year), I am not sure if the previous random effects are still trustworthy. Nevertheless, lme4::glFormula() seems the most promising method to construct a model matrix with random effects and new data. The method glmmTMB:::model.matrix.glmmTMB() does not allow random term predictors to be generated.

Test scripts I used

Data <- read.csv(paste0(
  "https://vincentarelbundock.github.io/Rdatasets/csv/", 
  "causaldata/thornton_hiv.csv")) |> na.omit()
library(glmmTMB)
library(marginaleffects)

## RMLE = FALSE
summary(Model <- glmmTMB(
  got ~ age * hiv2004  + (1 | villnum), data = Data))
"     AIC      BIC   logLik deviance df.resid 
  3578.8   3614.4  -1783.4   3566.8     2819 
 Groups   Name        Variance Std.Dev.
 villnum  (Intercept) 0.01648  0.1284  
 Residual             0.19878  0.4459  
Number of obs: 2825, groups:  villnum, 119
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.6317232  0.0262568  24.059  < 2e-16 ***
age          0.0022487  0.0006390   3.519 0.000433 ***
hiv2004     -0.0600519  0.1236294  -0.486 0.627151    
age:hiv2004  0.0005985  0.0032872   0.182 0.855530"
predict(Model, fast = F, se.fit = T) |> data.frame() |> head(10)
predict(Model, fast = T, se.fit = T) |> data.frame() |> head(10)
predict(Model, fast = T, se.fit = T, newdata = Data[1:10, ])
"fast=TRUE is not compatible with newdata/newparams/population-level prediction"
predict(Model, fast = F, se.fit = T, newdata = Data[1:10, ]) |> data.frame()
"         fit     se.fit
1  0.6246165 0.07009961
2  0.6178703 0.07029221
3  0.6943276 0.07116783
4  0.6875814 0.07082346
5  0.6223678 0.07015805
6  0.7145663 0.07249584
7  0.6489127 0.08420340
8  0.6088753 0.07062916
9  0.6223678 0.07015805
10 0.7033226 0.07170418 
identical results as fast = F/T"
predict(Model, re.form = NA, se.fit = T) |> data.frame() |> head(10)
predict(Model, re.form = ~ 0, se.fit = T) |> data.frame() |> head(10)
predict(Model, re.form = NA, se.fit = T, newdata = Data[1:10, ]) |> data.frame()
"         fit     se.fit
1  0.6811955 0.01720763
2  0.6744493 0.01806835
3  0.7509066 0.02031415
4  0.7441604 0.01915960
5  0.6789468 0.01747589
6  0.7711453 0.02436145
7  0.7054917 0.04985585
8  0.6654543 0.01945200
9  0.6789468 0.01747589
10 0.7599016 0.02201958"

predictions(Model) |> head(10)
"By default, standard errors for models of class `glmmTMB` are not
  calculated. For further details, see discussion at
  {https://github.com/glmmTMB/glmmTMB/issues/915}.
  Set `vcov = FALSE` or explicitly provide a variance-covariance-matrix for
  the `vcov` argument to calculate standard errors."
str(vcov(Model)) # a list of 1  $ cond for Conditional model
insight::get_varcov(Model, component = "all")
"                   (Intercept)           age       hiv2004   age:hiv2004
(Intercept)        6.894172e-04 -1.342988e-05 -4.046168e-04  1.036128e-05
age               -1.342988e-05  4.082644e-07  1.064802e-05 -3.217706e-07
hiv2004           -4.046168e-04  1.064802e-05  1.528423e-02 -3.908703e-04
age:hiv2004        1.036128e-05 -3.217706e-07 -3.908703e-04  1.080591e-05
d~(Intercept)     -6.302971e-06 -6.070752e-08  4.966683e-06 -1.728403e-07
theta_1|villnum.1  1.275741e-04  1.228732e-06 -1.005286e-04  3.498384e-06
                  d~(Intercept) theta_1|villnum.1
(Intercept)       -6.302971e-06      1.275741e-04
age               -6.070752e-08      1.228732e-06
hiv2004            4.966683e-06     -1.005286e-04
age:hiv2004       -1.728403e-07      3.498384e-06
d~(Intercept)      7.398080e-04     -2.905456e-04
theta_1|villnum.1 -2.905456e-04      1.322257e-02"
predictions(Model) |> head(10)
predictions(Model, vcov = vcov(Model)) |> head(10)
predictions(Model, vcov = vcov(Model)$cond) |> head(10)
predictions(Model, vcov = insight::get_varcov(Model, component = "all")) |> 
  head(10)
predictions(Model, newdata = Data[1:10, ], vcov = vcov(Model)) |> head(10)
" Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
    0.625    0.00853 73.2   <0.001   Inf 0.608  0.641
    0.618    0.00999 61.9   <0.001   Inf 0.598  0.637
    0.694    0.01495 46.4   <0.001   Inf 0.665  0.724
    0.688    0.01322 52.0   <0.001   Inf 0.662  0.713
    0.622    0.00900 69.2   <0.001   Inf 0.605  0.640
    0.715    0.02036 35.1   <0.001 894.2 0.675  0.754
    0.649    0.04742 13.7   <0.001 139.2 0.556  0.742
    0.609    0.01213 50.2   <0.001   Inf 0.585  0.633
    0.622    0.00900 69.2   <0.001   Inf 0.605  0.640
    0.703    0.01733 40.6   <0.001   Inf 0.669  0.737
EST match predict(Model, re.form = NULL, se.fit = T)
SE are much smaller, because uncertainty from sd(RE) not included.
Does not match any of the above predict() results
Need grouping variables in design matrix AND vcov with expanded dimensions
Only the top-left square for fixed effects of get_varcov() is used"
range(
  predict(Model, se.fit = T)$se.fit^2 - 
    predictions(Model, vcov = vcov(Model))$std.error^2)
"0.001429471 0.014137696"
range(
  predict(Model, se.fit = T, re.form = NA)$se.fit^2 - 
    predictions(Model, vcov = vcov(Model))$std.error^2)
"-0.0008082855  0.0005793517"
predictions(Model, re.form = NA, vcov = vcov(Model)) |> head(10)
" Estimate Std. Error    z Pr(>|z|)      S 2.5 % 97.5 %
    0.681     0.0172 39.6   <0.001    Inf 0.647  0.715
    0.674     0.0181 37.3   <0.001 1010.6 0.639  0.710
    0.751     0.0203 37.0   <0.001  991.2 0.711  0.791
    0.744     0.0192 38.8   <0.001    Inf 0.707  0.782
    0.679     0.0175 38.9   <0.001    Inf 0.645  0.713
    0.771     0.0244 31.7   <0.001  728.1 0.723  0.819
    0.705     0.0499 14.2   <0.001  148.6 0.608  0.803
    0.665     0.0195 34.2   <0.001  849.6 0.627  0.704
    0.679     0.0175 38.9   <0.001    Inf 0.645  0.713
    0.760     0.0220 34.5   <0.001  864.5 0.717  0.803
Match predict(Model, re.form = NA, se.fit = T) exactly
Note that EST changes when switching re.form = NULL to re.form = NA"
str(getME(Model, name = "X"))
" num [1:2825, 1:4] 1 1 1 1 1 1 1 1 1 1 ..."
str(getME(Model, name = "X")[ , 2]) # Named num [1:2825] 22 19 53 50
str(Model$frame$age) # int [1:2825] 22 19 53 50 ...
as.matrix(cbind(getME(Model, "X"), getME(Model, "Z")))[1:5, ]
"A n*(p + groups) design matrix, not Jacobian"
Model$fit$par
"         beta          beta          beta          beta         betad 
 0.6317231594  0.0022487446 -0.0600518660  0.0005984992 -1.6155417453 
        theta 
-2.0527384323
exp(betad) = exp(-1.6155417453) = 0.198783 = residual variance
exp(theta) = exp(-2.0527384323) = 0.1283829 = random effect SD"
Model$fit$parfull # adds BLUP RE estimates as b

# REML = TRUE
summary(Model <- glmmTMB(
  got ~ age * hiv2004  + (1 | villnum), data = Data, REML = TRUE))
"     AIC      BIC   logLik deviance df.resid 
  3612.7   3648.4  -1800.3   3600.7     2823 
 Groups   Name        Variance Std.Dev.
 villnum  (Intercept) 0.01676  0.1295  
 Residual             0.19899  0.4461  
Number of obs: 2825, groups:  villnum, 119
Dispersion estimate for gaussian family (sigma^2): 0.199 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.6317980  0.0263178  24.007  < 2e-16 ***
age          0.0022495  0.0006394   3.518 0.000434 ***
hiv2004     -0.0601110  0.1237032  -0.486 0.627017    
age:hiv2004  0.0006006  0.0032892   0.183 0.855124"
predict(Model, fast = F, se.fit = T) |> data.frame() |> head(10)
predict(Model, fast = T, se.fit = T) |> data.frame() |> head(10)
predict(Model, fast = T, se.fit = T, newdata = Data[1:10, ])
"fast=TRUE is not compatible with newdata/newparams/population-level prediction"
predict(Model, fast = F, se.fit = T, newdata = Data[1:10, ]) |> data.frame()
"         fit     se.fit
1  0.6243777 0.07029107
2  0.6176293 0.07048326
3  0.6941111 0.07135910
4  0.6873627 0.07101509
5  0.6221282 0.07034938
6  0.7143564 0.07268566
7  0.6487295 0.08437872
8  0.6086314 0.07081955
9  0.6221282 0.07034938
10 0.7031090 0.07189487"
predict(Model, re.form = NA, se.fit = T) |> data.frame() |> head(10)
predict(Model, re.form = NA, se.fit = T, newdata = Data[1:10, ]) |> data.frame()
"         fit     se.fit
1  0.6680648 0.01109493
2  0.6612708 0.01235834
3  0.7382701 0.01501680
4  0.7314760 0.01352887
5  0.6658001 0.01149779
6  0.7586522 0.01986940
7  0.6906714 0.04734265
8  0.6522120 0.01425101
9  0.6658001 0.01149779
10 0.7473288 0.01711696"
predictions(Model, vcov = vcov(Model)) |> head(10)
predictions(Model, vcov = vcov(Model)$cond) |> head(10)
predictions(Model, vcov = insight::get_varcov(Model, component = "all")) |> 
  head(10)
predictions(Model, newdata = Data[1:10, ], vcov = vcov(Model)) |> head(10)
"Uncertainty estimates cannot be computed for `glmmTMB` models with the
  `REML=TRUE` option. Either set `REML=FALSE` when fitting the model, or
  set `vcov=FALSE` when calling a `slopes` function to avoid this error"

get_se_manual <- function(m, random = FALSE) {
  s <- TMB::sdreport(m$obj, getJointPrecision = TRUE)
  Sigma <- solve(s$jointPrecision)
  if (isTRUE(random)) {
    XZ <- as.matrix(cbind(getME(m, "X"), getME(m, "Z")))
    keepvars <- c("beta", "b")
  } else {
    XZ <- as.matrix(getME(m, "X"))
    keepvars <- "beta"
  }
  keepvars <- rownames(Sigma) %in% keepvars
  Sigma <- Sigma[keepvars, keepvars]
  se <- sqrt(diag(XZ %*% (Sigma %*% t(XZ))))
  return(se)
}
get_Sigma <- function(m, random = FALSE) {
  s <- TMB::sdreport(m$obj, getJointPrecision = TRUE)
  Sigma <- solve(s$jointPrecision)
  if (isTRUE(random)) {
    # XZ <- as.matrix(cbind(getME(m, "X"), getME(m, "Z")))
    keepvars <- c("beta", "b")
  } else {
    # XZ <- as.matrix(getME(m, "X"))
    keepvars <- "beta"
  }
  keepvars <- rownames(Sigma) %in% keepvars
  Sigma <- Sigma[keepvars, keepvars]
  # se <- sqrt(diag(XZ %*% (Sigma %*% t(XZ))))
  return(Sigma)
}
get_XZ <- function(m, random = FALSE) {
  if (isTRUE(random)) {
    XZ <- as.matrix(cbind(getME(m, "X"), getME(m, "Z")))
  } else {
    XZ <- as.matrix(getME(m, "X"))
  }
  return(XZ)
}

summary(Model <- glmmTMB(
  got ~ age * hiv2004  + (1 | villnum), data = Data))

se1 <- get_se_manual(Model, random = TRUE)
se2 <- predict(Model, se.fit = TRUE)$se.fit
all.equal(unname(se1), unname(se2)) # TRUE

V <- vcov(Model)$cond # fixed effect vcov only
J <- model.matrix(Model) # fixed effect predictors only
se3 <- sqrt(diag(J %*% (V %*% t(J))))
se4 <- get_se_manual(Model, random = FALSE)
se5 <- predict(Model, se.fit = TRUE, re.form = NA)$se.fit
all.equal(unname(se3), unname(se4)) # TRUE
all.equal(unname(se4), unname(se5)) # TRUE

get_Sigma(Model)
get_Sigma(Model, random = TRUE) # 123 x 123 sparse Matrix of class "dgCMatrix"

get_coef.glmmTMB <- function(model, re.form = NULL, ...) {
  if (is.null(re.form)) {
    b <- model$fit$parfull
  } else {
    b <- model$fit$par
  }
  b <- setNames(as.vector(b), row.names(b))
  return(b)
}

set_coef.glmmTMB <- function(model, coefs, ...) {
  out <- model
  out$fit$parfull <- coefs
  return(out)
}

get_vcov.glmmTMB <- function(model, re.form = NULL, ...) {
  if (is.null(re.form)) {
    sdr <- TMB::sdreport(model$obj, getJointPrecision = TRUE)
    vcov <- solve(sdr$jointPrecision)
    keepvars <- c("beta", "b")
  } else {
    vcov <- Model$sdr$cov.fixed
    keepvars <- "beta"
  }
  keep <- rownames(vcov) %in% keepvars
  vcov <- vcov[keep, keep]
  return(vcov)
}

get_predict.glmmTMB <- function(
    model, newdata = insight::get_data(model), type = "response", ...) {
  # glmmTMB:::predict.glmmTMB()
  out <- stats::predict(model, newdata = newdata, type = type, ...)
  if (is.list(out)) out <- out$fit
  out <- data.frame(rowid = seq_len(nrow(newdata)), estimate = out)
  return(out)
}

get_vcov(Model)
dim(get_vcov(Model)) # 123 123
get_coef(Model)
set_coef(Model, coefs = get_coef(Model) + 1e-6)
get_predict(Model) # worked
get_predict(Model, se.fit = T)
"Error in View : Unable to extract predictions of type response from a model 
of class glmmTMB. Please report this problem, along with reproducible code and 
data on Github: https://github.com/vincentarelbundock/marginaleffects/issues"
stats::predict(Model, se.fit = T)

predictions(Model, vcov = as.matrix(get_vcov(Model))) |> head(10)
" Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
    0.625    0.00853 73.2   <0.001   Inf 0.608  0.641
    0.618    0.00999 61.9   <0.001   Inf 0.598  0.637
    0.694    0.01495 46.4   <0.001   Inf 0.665  0.724
    0.688    0.01322 52.0   <0.001   Inf 0.662  0.713
    0.622    0.00900 69.2   <0.001   Inf 0.605  0.640
    0.715    0.02036 35.1   <0.001 894.2 0.675  0.754
    0.649    0.04742 13.7   <0.001 139.2 0.556  0.742
    0.609    0.01213 50.2   <0.001   Inf 0.585  0.633
    0.622    0.00900 69.2   <0.001   Inf 0.605  0.640
    0.703    0.01733 40.6   <0.001   Inf 0.669  0.737
same as before, wrong SE"
View(attributes(predictions(Model, vcov = as.matrix(get_vcov(Model)))))
"Jacobian has 6 columns, not 4 + 119 = 123
vcov is 6 x 6, not 123 by 123 as in get_vcov(Model)"

marginaleffects:::get_predict.glmmTMB()
marginaleffects:::get_predict.glmmPQL()
insight::get_data(Model)
stats::predict(Model, se.fit = T)
lme4::findbars(Model$modelInfo$allForm$combForm)
model.matrix(
  lme4::formula(form), model.frame(object, ...), 
  contrasts.arg = object$modelInfo$contrasts)
Model$call$formula
model.matrix(Model$call$formula, data = Data)
model.matrix(Model)
model.frame(Model)
glmmTMB:::model.matrix.glmmTMB()
str(model.matrix(
  Model, component = c("cond", "zi", "disp"), part = c("fixed", "random")))
model.matrix(
  lme4:::formula.merMod(Model$modelInfo$allForm$formula), model.frame(Model), 
  contrasts = Model$modelInfo$contrasts)
"error trying to get slot frame from an object formula that is not an S4 object"
lme4::glFormula(
  Model$modelInfo$allForm$formula, data = model.frame(Model), 
  contrasts = Model$modelInfo$contrasts)
"list element $X has fixed predictors
$reTrms$Zt is sparse matrix, when transposed into matrix has random effects"
lme4::glFormula(
  got ~ age * hiv2004 + (1 + age | villnum), data = model.frame(Model), 
  contrasts = Model$modelInfo$contrasts)$reTrms$Zt |> t() |> as.matrix()
"First two columns for person 1. Row and column names are group index"

[vincentarelbundock]

@DrJerryTAO Thanks. This is known. I tried similar things here without success: https://github.com/vincentarelbundock/marginaleffects/pull/1023

[vincentarelbundock]

Quick question for @bbolker: Does the newparams argument in predict.glmmTMB() accept the full vector of parameters model$fit$parfull, or just model$fit$par?

My understanding is that I would need to make predictions using the full set of parameters if I want the equivalent of re.form=NULL when computing delta method standard errors for functions of predictions in glmmTMB models.

Background:

marginaleffects computes standard errors for predictions and functions of predictions via the delta method. First, we build a J matrix with the same number of rows as newdata and a number of columns equal to the number of coefficients in the model. Entries in the first column of J are (numerical) derivatives of the model's predictions with respect to the first coefficient, and so on. Then, we take the diagonal of JVJ' as the estimated variance of predictions for each observation.

In another thread, you recommended that I compute the numerical derivatives by feeding different newparams to the predict.glmmTMB() method instead of modifying the values hosted inside the model object.

I tried to do this in PR #1023 but got this warning:

In par[-random] <- par.fixed :
  number of items to replace is not a multiple of replacement length

[DrJerryTAO]

I just tested: predict.glmmTMB() accept only model$fit$par that does not include b for random effects.

predict(Model, newparams = get_coef(Model, re.form = NULL) + 1e-6)
"warning 12: In par[-random] <- par.fixed :
  number of items to replace is not a multiple of replacement length"
predict(Model, newparams = get_coef(Model, re.form = NA) + 1e-6) # Worked

[vincentarelbundock]

Closing in favor of https://github.com/vincentarelbundock/marginaleffects/issues/1064

Please note that there are critical issues in the computation of standard errors for glmmTMB models. After several attempts, I was not able to fix them due to lack of expertise in that model type and code base. Support for that package may be removed in a future release of marginaleffects.

[bwiernik]

I'll follow up in the other issue--happy to spend some time helping you diagnose issues if possible

[saudiwin]

Well this is discouraging. glmmTMB is the only MLE implementation of ordbetareg at the moment.