vincentarelbundock/marginaleffects

Ordbeta error in V0.22.03

Closed this issue · 12 comments

I installed the latest Github version of marginaleffects. When doing so a new error appeared that I have never seen before when fitting an ordbeta model with glmmTMB.

Error: The jacobian does not match the variance-covariance matrix.

This error does not appear in the CRAN version but does in the Github version.

Any ideas what is going on?

library(glmmTMB) 
library(marginaleffects) 
library(ordbetareg)

 model_data <- dplyr::select(pew,therm,
                            education="F_EDUCCAT2_FINAL",
                             region="F_CREGION_FINAL",
                             income="F_INCOME_FINAL")
ord_fit <- glmmTMB(formula=therm/100 ~ education + income + (1|region),
                                data=model_data,
                               family = ordbeta))

avg_predictions(ord_fit, variables = "education") 
Error: The jacobian does not match the variance-covariance matrix.

This error doesn't appear when using the CRAN version of marginaleffects.

Sorry. Updated example.

Thanks.

The problem is that get_coef() returns coefficient names with a cond~ prefix, but the names in get_vcov() do not include the same prefix, so there's a misalignment.

The fix should be pretty easy, and involve modifying this file:

https://github.com/vincentarelbundock/marginaleffects/blob/main/R/methods_glmmTMB.R

Unfortunately, I don't have much time right now, so not sure when I'll have time to resolve this issue.

I'll try to give it a shot! Should I just do a pull request?

That would be great!

This error is related to this issue to: #1189.

I attempted to match the columns but still no luck with this error (see below). I think it might have to do with ord beta family including upper and lower cutoff points.

#' @include get_vcov.R
#' @rdname get_vcov
#' @export
get_vcov.glmmTMB <- function(model, ...) {
    # Extract the full covariance matrix
    out <- stats::vcov(model, full = TRUE)

    # Extract the fixed-effect coefficient names from get_coef
    coef_names <- names(get_coef.glmmTMB(model))

    # Handle dispersion and conditional terms
    cleaned_coef_names <- gsub("^cond~", "", coef_names)  # Remove cond~ for conditional terms
    cleaned_coef_names <- gsub("^disp~", "d~", cleaned_coef_names)  # Map disp~ to d~ for dispersion terms

    # The 'upper cutoff' and 'lower cutoff' will remain in both, so no removal

    # Get the current row and column names from the covariance matrix
    current_names <- rownames(out)

    # Match cleaned coef_names with current names in the covariance matrix
    matched_indices <- match(current_names, cleaned_coef_names)

    # Replace row/column names only where there is a valid match
    valid_indices <- which(!is.na(matched_indices))
    
    if (length(valid_indices) > 0) {
        # Apply the correct names from coef_names to matched rows/columns in the covariance matrix
        rownames(out)[valid_indices] <- coef_names[matched_indices[valid_indices]]
        colnames(out)[valid_indices] <- coef_names[matched_indices[valid_indices]]
    } else {
        warning("No matching terms found between the covariance matrix and fixed-effect coefficients.")
    }

    return(out)
}

#' @include get_coef.R
#' @rdname get_coef
#' @export
get_coef.glmmTMB <- function(model, ...) {
    # Extract the fixed-effect coefficients
    out <- unlist(glmmTMB::fixef(model))

    # Apply the gsub logic to rename terms (cond~, disp~, etc.)
    names(out) <- gsub("^(cond|zi|disp)\\.", "\\1~", names(out))

    # No removal of "lower cutoff" and "upper cutoff" - they remain in place
    return(out)
}

@saudiwin maybe you have some insight?

I just tried replacing the methods in the R/methods_glmmTMB.R by the ones you posted above, and then reinstalled the package. It seems to work. What's the problem on your end?

I am still getting the Jacobian matrix error when running avg_predictions/avg_comparisons

What's your sessionInfo(). This is weird. I'm using the latest versions of everything and your code, and get no error.

Okay. My bad. I was working with global objects. I was able to load the package again and no errors! Woot!

Fantastic. Thanks so much for the code contribution!

Merged in main on Github.