能否用spaMM、glmmTMB或R-INLA拟合分站点行列一阶自回归空间模型？

问题背景 #

我有一个包含两个站点（站点1和站点2）的数据集，每个站点的数据按行列布局排列（例如3行×4列）。核心关注因子为trt，row、col和site为区组因子，数据集示例代码如下：

nsite <- 2
nrow <- 3
ncol <- 4
nperSite <- nrow*ncol
site <- factor(rep(1:nsite, each = nperSite))
d1 <- c(1, 2, 3, 4, 2, 3, 4, 1, 3, 4, 1, 2)
d2 <- c(1, 3, 4, 2, 3, 4, 2, 1, 4, 2, 1, 3)
trt  <- factor(c(d1,d2))
row  <- c(rep(1:nrow,each=ncol),rep(1:nrow,each=ncol))
col  <- c(rep(1:ncol,nrow),rep(1:ncol,nrow))
y    <- matrix(0,nrow=nsite*nperSite,1)
data <- data.frame(trt, row, col, site, y)

需求为构建空间模型，满足：

• 每个站点拥有独立的随机误差
• 每个站点内，对行和列的空间相关性采用**一阶自回归（AR(1)）**模型建模（即站点1对应行AR(r1)、列AR(c1)，站点2对应行AR(r2)、列AR(c2)）

各R包的实现方案 #

spaMM #

spaMM支持自定义相关结构，可轻松实现多站点独立的行列AR(1)模型，同时支持站点独立误差：

1. 构建站点与行、列的交互变量，区分不同站点的行/列组：

data$site_row <- interaction(data$site, data$row)
data$site_col <- interaction(data$site, data$col)

2. 使用fitme()拟合模型：

library(spaMM)
model_spaMM <- fitme(y ~ trt + 
                     corAR1(form = ~ 1 | site_row) + 
                     corAR1(form = ~ 1 | site_col) + 
                     (1 | site), 
                   data = data)

• corAR1(form = ~1 | site_row)：为每个站点的每一行拟合独立的AR(1)过程
• corAR1(form = ~1 | site_col)：为每个站点的每一列拟合独立的AR(1)过程
• (1 | site)：实现站点间独立的随机误差方差

glmmTMB #

glmmTMB支持原生的AR(1)相关结构，结合分组变量可实现多站点的独立行列AR(1)建模：

1. 确保数据按站点、行、列排序（保证AR(1)结构的有序性）：

data <- data[order(data$site, data$row, data$col), ]

2. 拟合模型：

library(glmmTMB)
model_glmmTMB <- glmmTMB(y ~ trt + 
                          (1 | site) +
                          ar1(row | site) + 
                          ar1(col | site), 
                        data = data)

• ar1(row | site)：在每个站点内，行维度的AR(1)空间相关
• ar1(col | site)：在每个站点内，列维度的AR(1)空间相关
• (1 | site)：实现站点间独立的随机误差方差

R-INLA #

R-INLA通过结构化随机效应实现贝叶斯框架下的空间AR(1)建模，适合复杂的多站点场景：

1. 为每个站点的行、列创建唯一索引变量：

data$row_idx <- as.integer(interaction(data$site, data$row))
data$col_idx <- as.integer(interaction(data$site, data$col))

2. 拟合模型：

library(INLA)
model_inla <- inla(y ~ trt +
                    f(row_idx, model = "ar1", hyper = list(theta = list(prior = "pc.prec", param = c(1, 0.01)))) +
                    f(col_idx, model = "ar1", hyper = list(theta = list(prior = "pc.prec", param = c(1, 0.01)))) +
                    f(site, model = "iid"),
                  data = data)

• f(row_idx, model="ar1")：为每个站点的行拟合独立的AR(1)过程（row_idx为站点-行的唯一标识）
• f(col_idx, model="ar1")：为每个站点的列拟合独立的AR(1)过程
• f(site, model="iid")：实现站点间独立的随机误差

总结 #

三个包均支持实现需求中的空间模型：

• spaMM和glmmTMB语法更贴近传统混合模型，操作便捷
• R-INLA适合贝叶斯框架下的复杂空间建模，可自定义先验分布

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