带自定义MAFC logit链接函数的二项式GLMM中emmeans::regrid()反变换失效的解决方法咨询

Great question! The issue here is that emmeans doesn't recognize your custom mafc.logit(2) link function out of the box, so it can't automatically apply the inverse transformation to get back to the [0.5, 1] probability scale. Let's walk through two reliable solutions to fix this:

Solution 1: Manual Inverse Transformation (Most Straightforward) #

Since you know the exact form of your inverse link function (P = 0.5 + 0.5 * plogis(eta)), you can compute the response-scale probabilities, standard errors, and confidence intervals manually using the linear-scale results from emmeans. Here's how:

# 1. Get linear-scale (eta) emmeans and contrasts
mod.emm <- emmeans(mod, ~group|stim_strength, at=list(stim_strength=c(.25,.75)))
contrasts_lin <- contrast(mod.emm, 'pairwise', simple = 'group', combine = TRUE, adjust = 'holm')

# 2. Define inverse link function and its derivative (for SE calculation via delta method)
inv_mafc_logit <- function(eta, m=2) {
  1/m + (1 - 1/m) * plogis(eta)
}

mu_eta <- function(eta, m=2) {
  # Derivative of inverse link w.r.t. eta (needed for SE scaling)
  (1 - 1/m) * plogis(eta) * (1 - plogis(eta))
}

# 3. Convert emmeans to probability scale
emm_df <- as.data.frame(mod.emm)
emm_df$prob <- inv_mafc_logit(emm_df$emmean, m=2)
emm_df$prob_se <- emm_df$SE * mu_eta(emm_df$emmean, m=2)
# Calculate confidence intervals using normal approximation
emm_df$prob_lower <- inv_mafc_logit(emm_df$emmean - qnorm(0.975)*emm_df$SE, m=2)
emm_df$prob_upper <- inv_mafc_logit(emm_df$emmean + qnorm(0.975)*emm_df$SE, m=2)

# View cleaned probability-scale results
cat("Probability-scale emmeans with CIs:\n")
print(emm_df[, c("stim_strength", "group", "prob", "prob_se", "prob_lower", "prob_upper")], digits=3)

# 4. Convert contrasts to response-scale differences
contrast_df <- as.data.frame(contrasts_lin)
# Contrast is (group1 eta) - (group2 eta), so response-scale difference is inv_mafc(eta1) - inv_mafc(eta2)
eta1 <- emm_df$emmean[emm_df$group == 1]
eta2 <- emm_df$emmean[emm_df$group == 2]
contrast_df$prob_diff <- inv_mafc_logit(eta1, m=2) - inv_mafc_logit(eta2, m=2)
# Use delta method for SE of the difference: SE = sqrt( (mu_eta(eta1)*SE1)^2 + (mu_eta(eta2)*SE2)^2 )
se1 <- emm_df$SE[emm_df$group == 1]
se2 <- emm_df$SE[emm_df$group == 2]
contrast_df$prob_diff_se <- sqrt( (mu_eta(eta1)*se1)^2 + (mu_eta(eta2)*se2)^2 )
# Calculate CIs for the difference
contrast_df$prob_diff_lower <- contrast_df$prob_diff - qnorm(0.975)*contrast_df$prob_diff_se
contrast_df$prob_diff_upper <- contrast_df$prob_diff + qnorm(0.975)*contrast_df$prob_diff_se

# View cleaned contrast results
cat("\nResponse-scale group differences with CIs:\n")
print(contrast_df[, c("stim_strength", "prob_diff", "prob_diff_se", "prob_diff_lower", "prob_diff_upper", "p.value")], digits=3)

Solution 2: Make emmeans Recognize Your Custom Link #

If you want emmeans to handle the transformation automatically (e.g., using type='response' or regrid()), you need to ensure your mafc.logit link function is properly defined with a linkinv (inverse link) attribute, and register it with emmeans.

First, double-check your mafc.logit function is defined correctly (this is the standard structure for custom GLM links):

mafc.logit <- function(m) {
  # Link function: maps probability to linear predictor eta
  linkfun <- function(mu) {
    if (any(mu <= 1/m | mu >= 1)) stop("mu must be between 1/m and 1")
    log( (mu - 1/m) / (1 - mu) )
  }
  
  # Inverse link function: maps eta back to probability
  linkinv <- function(eta) {
    1/m + (1 - 1/m) * plogis(eta)
  }
  
  # Derivative of inverse link (needed for model fitting and SE calculations)
  mu.eta <- function(eta) {
    (1 - 1/m) * plogis(eta) * (1 - plogis(eta))
  }
  
  # Validate eta values
  valideta <- function(eta) TRUE
  
  # Structure as a link-glm object with a clear name
  structure(
    list(linkfun = linkfun, linkinv = linkinv, mu.eta = mu.eta, valideta = valideta, 
         name = paste0("mafc.logit(", m, ")")),
    class = "link-glm"
  )
}

Once your link is properly defined, you can tell emmeans how to handle it by adding a transformation method. Run this once in your session:

# Register the inverse transformation for emmeans
emm_options(
  tran = list(
    "mafc.logit(2)" = list(
      linkfun = function(mu) log( (mu - 0.5)/(1 - mu) ),
      linkinv = function(eta) 0.5 + 0.5*plogis(eta),
      mu.eta = function(eta) 0.5*plogis(eta)*(1 - plogis(eta))
    )
  )
)

Now emmeans should recognize the link, and type='response' or regrid() should work as expected:

mod.emm <- emmeans(mod, ~group|stim_strength, at=list(stim_strength=c(.25,.75)))
# Get response-scale CIs
confint(mod.emm, type='response')
# Get response-scale contrasts
contrast(regrid(mod.emm), 'pairwise', simple='group', combine=TRUE, adjust='holm')

This should give you probabilities in the [0.5,1] range directly, without manual calculation.

Key Notes: #

• The delta method used in Solution 1 assumes normality of the linear predictor estimates, which is standard for GLMMs with large sample sizes.
• For Solution 2, if you use m values other than 2, you'll need to register a separate transformation for each (e.g., mafc.logit(3)).

内容的提问来源于stack exchange，提问作者noe