整合含周期性平滑与缺失值估计的JAGS模型技术求助

问题背景 #

需要构建一个JAGS模型，同时实现两个核心功能：

• 对week变量加入周期性平滑项（使用循环立方样条bs="cc"）
• 在模型内部估计含缺失值的协变量（age、雌雄计数），并纳入性别比（sr）的不确定性

目前能分别实现单功能模型，但无法整合：

• 用mgcv::jagam能生成带平滑项的JAGS代码，但会自动删除含缺失值的观测，且未考虑性别比的抽样不确定性
• 自定义JAGS代码能处理缺失值和性别比，但没有周期性平滑项

示例数据集 #

set.seed(123) # 固定随机种子便于复现
data <- data.frame(
  y = rbinom(100, 1, 0.5),
  week = runif(100, 1, 52),
  age = floor(runif(100, 1, 10)),
  females = floor(runif(100, 60, 100))
)
data$males <- floor(runif(100, 1, data$females))

# 加入缺失值
data$age[c(1,10,60)] <- NA
data$males[c(3,50,72)] <- NA
data$females[c(3,50,72)] <- NA

现有单功能模型 #

1. 带周期性平滑项的模型（jagam生成） #

此模型无法处理缺失值，且直接计算sr=males/females忽略了其不确定性：

library(mgcv)
data$sr <- data$males/data$females # 直接计算，未考虑抽样误差

# 生成带平滑项的JAGS代码
jd <- jagam(
  y ~ s(week, bs = "cc", k = 10, by = age) + s(week, bs = "cc", k = 10, by = sr),
  data = data,
  file = "model_jagam.txt",
  family = binomial,
  sp.prior = "gamma" # 平滑参数用gamma先验
)

2. 处理缺失值与性别比的模型（自定义） #

此模型无平滑项，仅用线性协变量：

model_linear <- function() {
  # 线性回归参数
  beta0 ~ dnorm(0, 0.1)
  beta_age ~ dnorm(0, 0.1)
  beta_sr ~ dnorm(0, 0.1)
  
  # 线性预测与响应模型
  for(obs in 1:n_obs) {
    mu[obs] <- beta0 + beta_age * age[obs] + beta_sr * sr[obs]
    p[obs] <- exp(mu[obs]) / (1 + exp(mu[obs]))
    y[obs] ~ dbin(p[obs], 1)
  }
  
  # 缺失age的先验模型
  for(obs in 1:n_obs){
    age[obs] ~ dnorm(age_mu, age_tau)
  }
  age_mu ~ dnorm(0, 0.0001)
  age_tau <- pow(age_sd, -2)
  age_sd ~ dunif(0, 100)
  
  # 性别比与缺失雌雄计数的模型
  for(obs in 1:n_obs){
    males[obs] ~ dbinom(sr[obs], females[obs])
    females[obs] ~ dpois(females_lambda)
    sr[obs] ~ dbeta(sr_alpha, sr_beta)
  }
  sr_alpha ~ dunif(.0001, 100)
  sr_beta ~ dunif(0.0001, 100)
  females_lambda ~ dunif(0, 100)
}

整合方案与完整代码 #

核心思路：提取jagam生成的平滑项结构（随机效应+罚项），嵌入到自定义的缺失值模型中，替换原线性协变量部分。

步骤说明 #

1. 手动构建平滑项的设计矩阵（避免jagam自动删除缺失值）
2. 在自定义JAGS模型中加入平滑项的参数（b为平滑系数，sp为平滑参数）
3. 将平滑项的预测值加入线性预测器mu[obs]
4. 保留原缺失值、性别比的模型部分

整合后的完整JAGS模型 #

# 手动构建平滑项的设计矩阵与罚项矩阵
library(mgcv)
library(Matrix)

# 构建带by变量的循环样条设计矩阵（保留所有观测）
smooth_age <- smoothCon(s(week, bs="cc", k=10, by=age), data=data, knots=list(week=c(1,52)))[[1]]
smooth_sr <- smoothCon(s(week, bs="cc", k=10, by=sr), data=data, knots=list(week=c(1,52)))[[1]]

X_age <- predict(smooth_age, data=data)
X_sr <- predict(smooth_sr, data=data)
X <- cbind(X_age, X_sr) # 合并设计矩阵

# 构建块对角罚项矩阵
S <- bdiag(smooth_age$S, smooth_sr$S)

# 准备JAGS输入数据
jags_data <- list(
  y = data$y,
  X = X,
  S = as.matrix(S),
  n_obs = nrow(data),
  n_b = ncol(X)
)

# 整合后的JAGS模型代码
model_combined <- function() {
  ### 平滑项参数
  # 平滑系数先验：服从多元正态分布，协方差由平滑参数与罚项矩阵控制
  b ~ dmnorm(zero[], prec_b[])
  prec_b <- sp[1] * S[1:(n_b/2), 1:(n_b/2)] + sp[2] * S[(n_b/2+1):n_b, (n_b/2+1):n_b]
  # 平滑参数的gamma先验
  for (k in 1:2) {
    sp[k] ~ dgamma(1, 0.001)
  }
  
  ### 线性预测与响应模型
  beta0 ~ dnorm(0, 0.1) # 截距项
  for(obs in 1:n_obs) {
    mu[obs] <- beta0 + inprod(X[obs, ], b[]) # 加入平滑项预测值
    p[obs] <- exp(mu[obs]) / (1 + exp(mu[obs]))
    y[obs] ~ dbin(p[obs], 1)
  }
  
  ### 缺失age的模型
  for(obs in 1:n_obs){
    age[obs] ~ dnorm(age_mu, age_tau)
  }
  age_mu ~ dnorm(0, 0.0001)
  age_tau <- pow(age_sd, -2)
  age_sd ~ dunif(0, 100)
  
  ### 性别比与缺失雌雄计数的模型
  for(obs in 1:n_obs){
    males[obs] ~ dbinom(sr[obs], females[obs])
    females[obs] ~ dpois(females_lambda)
    sr[obs] ~ dbeta(sr_alpha, sr_beta)
  }
  sr_alpha ~ dunif(.0001, 100)
  sr_beta ~ dunif(0.0001, 100)
  females_lambda ~ dunif(0, 100)
}

关键注意事项 #

• 平滑项设计矩阵需手动构建，避免jagam自动剔除含缺失值的观测
• 运行模型时，需将age、males、females、sr、b、sp等变量加入待估计参数列表
• 需为缺失变量提供合理初始值（比如用均值填充缺失的age，用观测值均值初始化sr）

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