Support `ordinal::clmm` and `ordinal::clmm2`

[vincentarelbundock]

I may be wrong but ordinal::clmm2 and ordinal::clmm seem to use predict method (maybe it's been updated). Would you mind checking if you are able to include these models? Thank you.

Originally posted by @hongyunsoo in #49 (comment)

Thanks @hongyunsoo I will take a look as soon as I find some time. But "real work" is crazy these days, so I can't promise a quick turnaround on this.

If you want give this a shot yourself, you can read about extensions here: https://marginaleffects.com/vignettes/extensions.html

[vincentarelbundock]

clmm does not have a predict() method.

clmm2 does have a predict method. Adding it to the sanity_models list seems to produce (untested) estimates, but no standard error

## Binomial example with the cbpp data from the lme4-package:
if(require(lme4)) {
    cbpp2 <- rbind(cbpp[,-(2:3)], cbpp[,-(2:3)])
    cbpp2 <- within(cbpp2, {
        incidence <- as.factor(rep(0:1, each=nrow(cbpp)))
        freq <- with(cbpp, c(incidence, size - incidence))
    })

    ## Fit with Laplace approximation:
    fm1 <- clmm2(incidence ~ period, random = herd, weights = freq,
                 data = cbpp2, Hess = 1)
    summary(fm1)

    ## Fit with the adaptive Gauss-Hermite quadrature approximation:
    fm2 <- clmm2(incidence ~ period, random = herd, weights = freq,
                 data = cbpp2, Hess = 1, nAGQ = 7)
    summary(fm2)
}

[vincentarelbundock]

I am interested in this and might eventually work on it, but closing to keep the repo clean. This issue is linked in the main Support New Models thread for future reference.

[hongyunsoo]

Thank you for looking into it! I really appreciate it. Best, Yun Soo From: Vincent Arel-Bundock ***@***.***> Sent: Tuesday, January 23, 2024 9:29 PM To: vincentarelbundock/marginaleffects ***@***.***> Cc: Yun Soo Hong ***@***.***>; Mention ***@***.***> Subject: Re: [vincentarelbundock/marginaleffects] Support `ordinal::clmm` and `ordinal::clmm2` (Issue #1003) External Email - Use Caution clmm does not have a predict() method. clmm2 does have a predict method. Adding it to the sanity_models list seems to produce (untested) estimates, but no standard error ## Binomial example with the cbpp data from the lme4-package: if(require(lme4)) { cbpp2 <- rbind(cbpp[,-(2:3)], cbpp[,-(2:3)]) cbpp2 <- within(cbpp2, { incidence <- as.factor(rep(0:1, each=nrow(cbpp))) freq <- with(cbpp, c(incidence, size - incidence)) }) ## Fit with Laplace approximation: fm1 <- clmm2(incidence ~ period, random = herd, weights = freq, data = cbpp2, Hess = 1) summary(fm1) ## Fit with the adaptive Gauss-Hermite quadrature approximation: fm2 <- clmm2(incidence ~ period, random = herd, weights = freq, data = cbpp2, Hess = 1, nAGQ = 7) summary(fm2) } - Reply to this email directly, view it on GitHub<#1003 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AQZGAZBZTDM235ZHR62E2Q3YQBWXFAVCNFSM6AAAAABCEC4UR2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMBXGI2DKNZVHE>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[DrJerryTAO]

Is no standard error reported in predict a problem? According to the the extension instruction, it does not seem to require standard error methods. I found that clmm2 objects actually use predict.clm2. I think the challenge is a proper manipulation of predict.clm2 for new data.

According to ?predict.clm2, newdata must have a response column otherwise errors, unlike many other prediction methods. If the original response column is supplied, it gives a vector of predicted probability of only the observed category. To produce probabilities of other categories, the help file gives an example of generating all possible response categories for each observation through expand.grid(). It still gives a vector, length equal to levels times observations.

[vincentarelbundock]

It's not a problem as we can compute standard errors in house.

I am still interested in this if someone prepared a pull request, but I don't have time to work on it myself.

[carly-gray-19]

I was able to get this type of result for my clmm2 model using the clarify package, which calculates standard errors using simulations rather than the delta method or bootstrapping. A lot of the clarify package appears to draw from marignaleffects under the hood. However, to work with a clmm2 model, I had to create a custom function that I could then implement via sim_apply() to get my quantities of interest, but this wasn't too tough to do.

I hope that's okay to point to this other package. Maybe there's a way to integrate the functionality of clarify with marginaleffects eventually?

[vincentarelbundock]

@carly-gray-19, thanks for the nudge. I tried to get this to work but it seems tricky.

FWIW, I think you should be very careful in using a framework like clarify or marginaleffects (maybe @ngreifer will tell you otherwise). To see why, install this branch, and then restart your R session completely:

remotes::install_github("vincentarelbundock/marginaleffects@clmm2")

Download data and fit a model:

library(marginaleffects)
library(ordinal)
dat_soup <- get_dataset("soup", "ordinal")
dat_soup <- subset(dat_soup, as.numeric(dat_soup$RESP) <= 24)
dat_soup$SURENESS <- factor(dat_soup$SURENESS, ordered = TRUE)
dat_soup$RESP <- factor(dat_soup$RESP)
dat_soup$RESP <- dat_soup$RESP[drop = TRUE]  # Drop unused levels
dat_soup$PROD <- factor(dat_soup$PROD)  # Ensure PROD is a factor
mod_soup <- clmm2(SURENESS ~ PROD, random = RESP, data = dat_soup,
                  link = "probit", Hess = TRUE, method = "ucminf",
                  threshold = "symmetric")

At first glance, everything seems to work fine:

avg_predictions(mod_soup)

 Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
    0.335     0.0216 15.5   <0.001 177.6 0.293  0.378

avg_comparisons(mod_soup)

 Estimate Std. Error    z Pr(>|z|)    S  2.5 % 97.5 %
    0.142     0.0221 6.42   <0.001 32.7 0.0987  0.185

However, when we look at unit-level predictions, something is wrong:

predictions(mod_soup)

 Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
    0.219         NA    NA       NA   NA    NA     NA
    0.178     0.0169 10.53   <0.001 83.8 0.145  0.211
    0.184         NA    NA       NA   NA    NA     NA
    0.623     0.0637  9.78   <0.001 72.6 0.498  0.748
    0.184         NA    NA       NA   NA    NA     NA
--- 190 rows omitted. See ?print.marginaleffects ---
    0.623     0.0637  9.78   <0.001 72.6 0.498  0.748
    0.219         NA    NA       NA   NA    NA     NA
    0.623     0.0637  9.78   <0.001 72.6 0.498  0.748
    0.178     0.0169 10.53   <0.001 83.8 0.145  0.211
    0.238         NA    NA       NA   NA    NA     NA

Now, note that marginaleffects does support the same kind of simulation-based inference as clarify, via the inferences() function. But here again, we get similar problems, with many confidence intervals of width 0:

predictions(mod_soup) |> inferences(method = "simulation")

 Estimate 2.5 % 97.5 %
    0.219 0.219  0.219
    0.178 0.140  0.202
    0.184 0.184  0.184
    0.623 0.502  0.739
    0.184 0.184  0.184
--- 190 rows omitted. See ?print.marginaleffects ---
    0.623 0.502  0.739
    0.219 0.219  0.219
    0.623 0.502  0.739
    0.178 0.140  0.202
    0.238 0.238  0.238
Type: response

IIUC, clarify uses the set_coef() methods in marginaleffects as one of their key components. Unless they have modified that code, it is quite likely that they are hitting the same problems. I strongly recommend that you conduct due diligence before moving ahead with interpretation.

I would be very happy to get any help on solving this issue. It probably involves fixing the set_coef() method in this branch: #1610

[vincentarelbundock]

FWIW, bootstrapping appears to produce more reasonable output:

predictions(mod_soup) |> inferences(method = "rsample")

But I don't know this package or model well enough to be confident in moving ahead with this implementation. Would love it if an expert could give me pointers.

[ngreifer]

@carly-gray-19 Could you provide some more detail on the function you wrote for sim_apply()? That might help us integrate it better into marginaleffects. I'm not familiar with this type of model, either.

[carly-gray-19]

Hi @vincentarelbundock and @ngreifer , thanks for your responses and diligence in trying to figure this out. I am relatively new to this type of model and to these packages, so I definitely welcome your caution.

Below is the function I wrote for sim_apply(). I should mention that my particular use case is a model for which I created an interaction term using the interact() function, due to a specific level of the interaction being missing by design. You'll see below that I specify a reference level for the interaction in the function arguments and pass a reference level that uses the interact() syntax when I run the function.

set.seed(1234)
sim_agreement_03 <- clarify::sim(agreement_out_03, n=10000) ## simulate

sim_interaction_pred_prob_CI <- function(coefs, sim_model, reference_cat) {  
   ## coefs will automatically get extracted by the simulated model estimated by clarify::sim() when that model is passed in the first position of sim_apply()
   ## sim_model is that same simulated model and needs to be specified again
   ## reference_at is the reference level of the interaction term so that we can get estimated probabilities for that level

  n_theta <- length(sim_model$fit$Theta)
  
  thetas <- c(-1e5,unname(coefs[1:n_theta]), 1e5)
  theta_names <- sim_model$fit$lev

  betas <- c(0,sim_model$fit$beta)
  beta_names <- names(betas)
  beta_names[1] <- reference_cat

  F <- plogis
  p_est <- matrix(NA_real_,
                  nrow = length(betas),
                  ncol = length(theta_names),
                  dimnames = list(beta_names, as.character(theta_names)))

  # Fill matrix
  for (i in seq_along(betas)) {
    for (j in seq_along(theta_names)) {
      p_est[i, j] <- F(thetas[j + 1] - betas[i]) - F(thetas[j] - betas[i])
    }
  }
  out <- as.vector(p_est)
  names(out) <- paste(rep(dimnames(p_est)[[2]], each=nrow(p_est)),
                      rep(dimnames(p_est)[[1]], times=ncol(p_est)),
                      sep = " - ")
  out
}

# simulate average predictions
est <- sim_apply(sim_agreement_03, FUN=sim_interaction_pred_prob_CI, sim_model = sim_agreement_03, reference_cat = "resp_policy_intparent:Full school day")

FWIW, I ran into CIs of zero width a few times when trying to troubleshoot this function. If I recall, the problem arose when I tried to only pass the simulated model (sim_agreement_03) without specifying it twice in the function arguments, presumably because of how the $fit and/or $coefs elements get extracted (see sim_apply() documentation). Once I specified it twice in the function, I was able to get actual confidence intervals. This might be because of the set_coef() method @vincentarelbundock mentioned.

[ngreifer]

Speaking only from the perspective of the clarify author, this is what I see:

Your estimation function uses coefs but not the model fit. set_coef() is never called here. You essentially wrote your own prediction function that is a function only of the coefficients. This is ideal for clarify because you don't need to rely on a pre-written method for extracting predictions from the model object. You mention needing to specify the simulated model "twice"; the model is only automatically supplied to the estimation function supplied to sim_apply() when the argument is called fit, but you called it sim_model; this is explained in the sim_apply() documentation. It doesn't matter here, though; you should only supply the model if you need to use a function on the model itself; here, your function only depends on the coefficients, which can be passed directly as you did. If you wanted to use, e.g., predict() on the model inside sim_interaction_pred_prob_CI(), then you would need to ensure the argument is named fit. This makes sim_apply() use set_coef() to inject the simulated coefficients into the model object.

For marginaleffects, you are relying on the internal functionality to extract predictions from the model and set coefficients into the model object. This is where things can be error-prone and why you might be getting weird issues with marginaleffects.

[carly-gray-19]

Thank you for taking a look at this.

You jogged my memory on why I passed the model to the function twice. The first pass of the model to the function extracts coefficients automatically because I set that argument to coefs, as you identified and is described in the documentation. The second pass of the model, which is to the sim_model argument, gets used to:

1. Extract Theta estimates without also extracting beta estimates to easily identify the number of cutpoints, i.e., n_theta <- length(sim_model$fit$Theta)
2. Extract the names of the Theta cutpoints without having to manipulate the concatenated strings (e.g., "Disagree|Neutral") that get stored incoefs, i.e., theta_names <- sim_model$fit$lev
3. Extract coefficients. This line could have just used coefs -- this was an arbitrary choice. I.e., betas <- c(0,sim_model$fit$beta)

Hope that clarifies why I wrote the function the way I did.

It sounds like you're expecting this function would work as I'm intending it to, meaning that it returns predicted probabilities for each level of the ordinal outcome variable and confidence intervals around these probabilities estimated via simulation? My output is generally what I'd expect based on the raw percentages of each outcome value by the levels of my interaction term, as well as the differences in sample sizes among those groups comprising my interaction term. However, I'd greatly appreciate you letting me know if I'm overlooking anything about the function I wrote or something more fundamental.

[ngreifer]

My only issue is with step 3. It's not arbitrary whether coefs or sim_model$fit$beta is used; coefs contains a draw of simulated coefficients from their distribution, and sim_model$fit$beta contains the coefficients in the original model fit.

If the set_coef() method is working correctly, you should be able to use the following modification of your function (I made some modifications for clarity and performance as well):

sim_interaction_pred_prob_CI <- function(fit, reference_cat) {  
   ## reference_at is the reference level of the interaction term so that we can get estimated probabilities for that level

  thetas <- c(-1e5, unname(fit$Theta), 1e5)
  theta_names <- fit$lev

  betas <- c(0, fit$beta)
  beta_names <- c(reference_cat, names(fit$beta))

  .F <- plogis
  p_est <- matrix(NA_real_,
                  nrow = length(betas),
                  ncol = length(theta_names),
                  dimnames = list(beta_names, as.character(theta_names)))

  # Fill matrix
  for (j in seq_along(theta_names)) {
    p_est[, j] <- .F(thetas[j + 1] - betas) - .F(thetas[j] - betas)
  }
  
  p_est |>
    as.vector() |>
    setNames(paste(rep(colnames(p_est), each = nrow(p_est)),
                   rep(rownames(p_est), times = ncol(p_est)),
                   sep = " - "))
}

This assumes set_coef() correctly injects the simulated coefficients into the model object. I don't know if that's true because it looks like set_coef() ignores the Theta component (but I have no idea what this model is).