使用mgcv拟合带权重的logistic GAM后如何正确绘制QQ图？

在R中使用mgcv包拟合带权重的logistic GAM模型时，调用qq.gam()生成QQ图会触发警告，提示非整数二项式分母导致分位数计算错误。复现代码如下：

library(mgcv)

set.seed(123)
n <- 100
x <- runif(n, 0, 1)
y <- rbinom(n, size = 1, prob = plogis(2 * x - 1))
y <- cbind(y, ifelse(y == 1, 0, 1)) # 转换为两列响应矩阵
weights <- runif(n, 0, 1)

# 拟合带权重的二项式GAM
mod <- gam(y ~ s(x), family = binomial, weights = weights/mean(weights))

# 使用mgcv的qq.gam生成QQ图
qq.gam(mod)

执行后出现的警告信息：

1: In fam$qf(U, object$fitted.values, object$prior.weights, object$sig2) :
non-integer binomial denominator : incorrect quantiles

问题根源在于mgcv内置的fam$qf函数会将模型的prior.weights作为qbinom()的size参数，但此处使用的是经weights/mean(weights)归一化后的非整数权重（Wood等人建议该归一化方式以避免改变对数似然的整体幅度，详见gam帮助文档）。该函数代码如下：

function (p, mu, wt, scale) 
{
    if (all.equal(wt, ceiling(wt)) != TRUE) {
        wt <- ceiling(wt)
        warning("non-integer binomial denominator: quantiles incorrect")
    }
    qbinom(p, wt, mu)/(wt + as.numeric(wt == 0))
}

针对带非整数权重的二项式GAM，可通过以下两种方法生成正确的QQ图：

方法1：手动基于加权残差绘制QQ图 #

核心思路是利用残差的渐近正态性，直接使用皮尔逊残差或偏差残差，结合正态理论分位数构建QQ图，无需依赖mgcv内置的分位数计算逻辑：

library(mgcv)
library(ggplot2)

# 拟合模型（同原代码）
set.seed(123)
n <- 100
x <- runif(n, 0, 1)
y <- rbinom(n, size = 1, prob = plogis(2 * x - 1))
y <- cbind(y, ifelse(y == 1, 0, 1))
weights <- runif(n, 0, 1)
mod <- gam(y ~ s(x), family = binomial, weights = weights/mean(weights))

# 提取皮尔逊残差
pearson_resid <- residuals(mod, type = "pearson")

# 计算理论正态分位数（采用Blom变换避免极端值）
p <- (seq_along(pearson_resid) - 0.375) / (length(pearson_resid) + 0.25)
theoretical_q <- qnorm(p)

# 绘制QQ图
qq_data <- data.frame(Theoretical = theoretical_q, Pearson_Residuals = pearson_resid)
ggplot(qq_data, aes(x = Theoretical, y = Pearson_Residuals)) +
  geom_point(size = 2, alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0, color = "darkred", linetype = "dashed") +
  labs(title = "QQ图：加权二项式GAM的皮尔逊残差",
       x = "理论正态分位数",
       y = "皮尔逊残差") +
  theme_bw()

方法2：修改二项式族的分位数函数适配非整数权重 #

如果希望继续使用qq.gam()的原生接口，可以临时替换二项式族的qf方法，针对非整数权重使用正态近似计算分位数（大样本下该近似合理）：

library(mgcv)

# 定义适配非整数权重的二项式族
weighted_binom <- binomial()
weighted_binom$qf <- function(p, mu, wt, scale) {
  # 基于正态近似计算分位数：二项分布均值=mu*wt，方差=wt*mu*(1-mu)
  mean_val <- mu * wt
  var_val <- wt * mu * (1 - mu)
  # 将分位数归一化到原响应尺度（与原生函数输出一致）
  qnorm(p, mean = mean_val, sd = sqrt(var_val)) / (wt + as.numeric(wt == 0))
}

# 用修改后的族重新拟合模型
set.seed(123)
n <- 100
x <- runif(n, 0, 1)
y <- rbinom(n, size = 1, prob = plogis(2 * x - 1))
y <- cbind(y, ifelse(y == 1, 0, 1))
weights <- runif(n, 0, 1)
mod_weighted <- gam(y ~ s(x), family = weighted_binom, weights = weights/mean(weights))

# 调用qq.gam生成QQ图，无警告
qq.gam(mod_weighted)

方法说明 #

• 方法1无需修改内置函数，直接基于残差的统计性质，稳定性强，适合大多数分析场景。
• 方法2保留了qq.gam()的使用习惯，但需注意：当权重与整数偏差较大时，正态近似的可靠性依赖于样本量足够大的前提。

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