使用R中mgcv包拟合GAM时的后验预测检查问题

问题分析 #

你遇到的incorrect number of dimensions错误，是因为performance包的check_model函数目前对mgcv中ocat（有序类别）族的GAM模型支持不完善——该函数调用simulate时，无法正确处理ocat模型返回的矩阵格式模拟结果。

手动实现后验预测检查可视化 #

我们可以直接用mgcv的simulate函数生成模拟数据，结合tidyverse和ggplot2手动实现类似performance包的可视化效果，步骤如下：

1. 加载包并拟合模型

library(mgcv)
library(tidyverse)
library(ggplot2)

# 整理数据
df <- mtcars |>
  tibble::as_tibble() |>
  dplyr::mutate(carb = as.integer(carb))

# 拟合有序分类GAM
model <- mgcv::gam(
  carb ~ s(wt, k=5),
  data = df,
  family = mgcv::ocat(R=8)
)

2. 生成并整理模拟数据

# 生成100次模拟（次数越多，分布越稳定）
nsim <- 100
simulated_data <- mgcv::simulate(model, nsim = nsim) |>
  tibble::as_tibble() |>
  tidyr::pivot_longer(cols = dplyr::everything(), 
                      names_to = "simulation_id", 
                      values_to = "predicted_carb")

# 整理观测数据
observed_data <- df |>
  dplyr::select(carb) |>
  dplyr::mutate(type = "Observed")

# 合并模拟与观测数据
plot_data <- simulated_data |>
  dplyr::mutate(type = "Simulated") |>
  dplyr::bind_rows(observed_data)

3. 绘制观测vs模拟的分布对比图

这个版本贴近performance的示例效果，用抖动点展示观测值，密度曲线展示模拟数据的分布：

ggplot2::ggplot(plot_data, ggplot2::aes(x = carb)) +
  # 模拟数据的密度曲线（半透明填充）
  ggplot2::geom_density(data = dplyr::filter(plot_data, type == "Simulated"),
                        ggplot2::aes(x = predicted_carb),
                        alpha = 0.3, fill = "#619CFF") +
  # 观测数据的抖动点（避免重叠）
  ggplot2::geom_jitter(data = dplyr::filter(plot_data, type == "Observed"),
                       ggplot2::aes(y = 0),
                       width = 0.2, height = 0.01, color = "#F8766D", size = 2) +
  ggplot2::labs(x = "Number of Carburetors (carb)", 
                y = "Density", 
                title = "Posterior Predictive Check") +
  ggplot2::theme_minimal()

4. 替代方案：计数对比图

如果更关注每个类别下的观测与模拟计数差异，可以用以下代码：

# 计算观测值的计数
obs_counts <- df |>
  dplyr::count(carb, name = "observed_count")

# 计算模拟值的平均计数（除以模拟次数）
sim_counts <- simulated_data |>
  dplyr::count(predicted_carb, name = "simulated_total") |>
  dplyr::mutate(simulated_avg = simulated_total / nsim)

ggplot2::ggplot() +
  # 模拟值的平均计数
  ggplot2::geom_point(data = sim_counts,
                      ggplot2::aes(x = predicted_carb, y = simulated_avg),
                      color = "#619CFF", size = 2) +
  # 观测值的计数
  ggplot2::geom_point(data = obs_counts,
                      ggplot2::aes(x = carb, y = observed_count),
                      color = "#F8766D", size = 3, shape = 18) +
  ggplot2::labs(x = "Number of Carburetors (carb)", 
                y = "Count", 
                title = "Observed vs Simulated Counts") +
  ggplot2::theme_minimal()

关于与performance结果不一致的说明 #

performance包对支持的模型，可能会采用参数化预测或不同的参数抽样逻辑（比如基于贝叶斯后验的抽样），而mgcv::simulate默认是基于模型拟合的协方差矩阵抽样参数后生成模拟数据。如果需要对齐逻辑，可以尝试：

• 调整mgcv::simulate的参数（比如设置parametric = TRUE或use.u = FALSE）
• 增加模拟次数（比如nsim = 500）来稳定模拟分布

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