运行BYM2模型遇数据初始化错误，寻求技术协助

问题背景 #

运行BYM2空间模型时出现变量识别错误，伴随Stan文件格式警告，无法完成采样。

运行代码 #

library("rstan")
library("rstudioapi")
library("parallel")
library("brms")

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores()) # 原代码detectgores()为拼写错误，需改为detectCores()

library(pkgbuild) # load package
find_rtools() # 若安装Rtools 3.5应返回TRUE

# 拟合icar.stan模型到NYC普查区邻域地图
install.packages('tidyverse', dependencies = TRUE)
install.packages('rstanarm', dependencies = TRUE)
library(rstan);
library(tidyverse)
library(rstanarm)

# 注：此处定义的Stan代码字符串未被使用，实际调用外部文件bym2_predictor_plus_offset.stan
"data {
int<lower=0> N;
int<lower=0> N_edges;
int<lower=1, upper=N> node1[N_edges]; // node1[i]与node2[i]相邻
int<lower=1, upper=N> node2[N_edges]; // 且node1[i] < node2[i]

int<lower=0> y[N]; // 计数结果
vector<lower=0>[N] E; // 暴露量
int<lower=1> K; // 协变量数量
matrix[N, K] x; // 设计矩阵

real<lower=0> scaling_factor; // 空间效应方差的缩放因子
}
transformed data {
vector[N] log_E = log(E);
}
parameters {
real beta0; // 截距
vector[K] betas; // 协变量系数

real<lower=0> sigma; // 总标准差
real<lower=0, upper=1> rho; // 非结构化与结构化空间方差的比例

vector[N] theta; // 异质性效应
vector[N] phi; // 空间效应
}
transformed parameters {
vector[N] convolved_re;
// 确保每个分量的方差近似为1
convolved_re = sqrt(1 - rho) * theta + sqrt(rho / scaling_factor) * phi;
}
model {
y ~ poisson_log(log_E + beta0 + x * betas + convolved_re * sigma); // 加入协变量

// phi的先验（比例形式）
target += -0.5 * dot_self(phi[node1] - phi[node2]);

beta0 ~ normal(0.0, 1.0);
betas ~ normal(0.0, 1.0);
theta ~ normal(0.0, 1.0);
sigma ~ normal(0, 1.0);
rho ~ beta(0.5, 0.5);
// phi的软和为零约束
sum(phi) ~ normal(0, 0.001 * N); // 等价于mean(phi) ~ normal(0,0.001)
}
generated quantities {
real logit_rho = log(rho / (1.0 - rho));
vector[N] eta = log_E + beta0 + x * betas + convolved_re * sigma; // 线性预测器
vector[N] mu = exp(eta); // 均值
}"
options(mc.cores = parallel::detectCores())

library(INLA)

source("mungecardata4stan.R")
source("iran_data.R")
y = data$y;
E = data$E;
K = 1;
x = 0.1 * data$x;

nbs = mungeCARdata4stan(data$adj, data$num);
N = nbs$N;
node1 = nbs$node1;
node2 = nbs$node2;
N_edges = nbs$N_edges;
adj.matrix = sparseMatrix(i=nbs$node1,j=nbs$node2,x=1,symmetric=TRUE)
Q= Diagonal(nbs$N, rowSums(adj.matrix)) - adj.matrix
Q_pert = Q + Diagonal(nbs$N) * max(diag(Q)) * sqrt(.Machine$double.eps)
Q_inv = inla.qinv(Q_pert, constr=list(A = matrix(1,1,nbs$N),e=0))
scaling_factor = exp(mean(log(diag(Q_inv))))
scot_stanfit = stan("bym2_predictor_plus_offset.stan", data=list(N,N_edges,node1,node2,y,x,E,scaling_factor), warmup=5000, iter=6000); 

错误与警告信息 #

错误信息 #

Error in new_CppObject_xp(fields$.module, fields$.pointer, …) : 
  Exception: 变量不存在；处理阶段=data initialization；变量名=N；基础类型=int（位于‘string’第3行第2至17列）
未能创建采样器；未完成采样

警告信息 #

In readLines(file, warn = TRUE) :
  ‘C:\Users\Uaer\Downloads\bym2_predictor_plus_offset.stan’文件存在不完整的最后一行

解决方案 #

1. 修复数据列表命名问题
   Stan要求传入的data必须是命名列表，原代码中list(N,N_edges,node1,node2,y,x,E,scaling_factor)未指定变量名，导致Stan无法匹配模型中的变量定义。同时模型需要的K变量也未传入，修改后代码如下：

scot_stanfit = stan("bym2_predictor_plus_offset.stan", 
                    data=list(N=N, 
                              N_edges=N_edges, 
                              node1=node1, 
                              node2=node2, 
                              y=y, 
                              x=x, 
                              E=E, 
                              scaling_factor=scaling_factor,
                              K=K),
                    warmup=5000, iter=6000); 

2. 修复Stan文件格式问题
   打开bym2_predictor_plus_offset.stan文件：

• 检查最后一行是否为完整语句（如闭合的}）
• 在文件末尾添加一个空行，消除“不完整最后一行”的警告
  保存后重新运行代码。

3. 修正拼写错误
   原代码中parallel::detectgores()为拼写错误，需改为parallel::detectCores()，否则无法正确获取可用核心数。

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