如何用R的MARSS包构建多物种多站点联合MARSS模型？

问题背景 #

我正在用R的MARSS包构建模型，目标是同时对多物种多站点的长期监测数据建模。研究设计为：30年间对9个站点的9个物种进行年度采样。目前已实现单站点多物种、多站点单物种的模型运行（代码见下文），但不清楚如何将二者整合，也不确定这种整合是否可行；同时也在考虑是否应该为每个站点单独运行模型。最终希望获取每个站点独立的B矩阵，用于比较不同站点的群落稳定性指标。

现有实现代码 #

library(MARSS)

dat <- data.frame(site = rep(1:2,each = 20),          # 两个独立站点
                  year = rep(c(1:20),2),              # 采样年份（1-20）
                  spp1 = sample(50,40, replace = T),  # 物种1的模拟计数数据
                  spp2 = sample(20,40, replace = T),  # 物种2的模拟计数数据
                  spp3 = sample(30,40, replace = T))  # 物种3的模拟计数数据

# 将计数数据转换为宽格式（行=物种，列=年份）
site1 <- t(dat[1:20,3:5])
site2 <- t(dat[21:40,3:5])

rownames(site1) <- paste0("site1_", rownames(site1))
rownames(site2) <- paste0("site2_", rownames(site2))

y <- rbind(site1,site2)

# 运行单物种双站点模型
y1 <- rbind(y[1,], y[4,])

# Z矩阵格式化
z.model1 <- matrix(c(1,0,0,1),ncol = 2)

# A矩阵格式化
a.model1 <- matrix(list(0),2,1)
a.model1[2,1] <- c("b")

# U矩阵格式化
u.model1 <- matrix(list(0),2,1)
u.model1[1:2,1] <- c("a","b")

# 整合模型各部分
model.1 <- list(Z = z.model1,                
                A = a.model1,                
                Q = "diagonal and unequal", 
                R = "diagonal and equal",                     
                U = u.model1,                
                tinitx = 1) 

# 运行模型
m1.out <- MARSS(y1, model = model.1, method = "BFGS")


# 运行单站点多物种模型
y2 <- y[1:3,]

# Z矩阵格式化
z.model2 <- matrix(0,3,3)
diag(z.model2) <- 1

# U矩阵格式化
u.model2 <- matrix(list(0),3,1)
u.model2[1:3,1] <- c("a","b","c")

# B矩阵格式化
b.model2 <- matrix(c("b11","b12","b13",
                     "b21","b22","b23",
                     "b31","b32","b33"),
                   ncol = 3)

# 整合模型各部分
model.2 <- list(Z = z.model2,                
                B = b.model2,
                Q = "diagonal and unequal",  
                R = "diagonal and equal",                     
                U = u.model2,                
                tinitx = 1) 

# 运行模型
m2.out <- MARSS(y2, model = model.2, method = "BFGS")

解决方案思路 #

1. 分站点独立建模（推荐优先尝试） #

这种方案完全合理，且优势明显：

• 实现简单：直接复用现有单站点多物种模型代码，循环遍历9个站点即可
• 参数易解释：每个站点得到独立的B矩阵、U截距、Q和R参数，完全契合“比较站点间群落稳定性”的目标
• 计算效率高：参数数量远少于整合模型，收敛速度更快，不易出现拟合失败
• 灵活性强：可根据单个站点的数据特征调整模型结构（比如某站点物种数据缺失时，单独调整Z矩阵）

循环处理多站点的代码示例：

# 整理数据为按站点分组的列表
site_data_list <- list(
  site1 = t(dat[dat$site == 1, 3:5]),
  site2 = t(dat[dat$site == 2, 3:5])
)

# 循环运行模型
model_results <- lapply(site_data_list, function(site_y) {
  n_spp <- nrow(site_y)
  # 构建观测矩阵Z（对角矩阵，每个状态对应一个物种观测）
  z_mat <- diag(n_spp)
  # 构建截距矩阵U（每个物种对应独立截距）
  u_mat <- matrix(paste0("u_", 1:n_spp), nrow = n_spp, ncol = 1)
  # 构建物种相互作用矩阵B（全参数化）
  b_mat <- matrix(paste0("b_", rep(1:n_spp, each = n_spp), 1:n_spp), 
                  nrow = n_spp, ncol = n_spp)
  # 定义模型
  site_model <- list(
    Z = z_mat,
    B = b_mat,
    U = u_mat,
    Q = "diagonal and unequal",
    R = "diagonal and equal",
    tinitx = 1
  )
  # 运行模型
  MARSS(site_y, model = site_model, method = "BFGS")
})

# 提取每个站点的B矩阵
site_B_matrices <- lapply(model_results, function(res) {
  coef(res, type = "matrix")$B
})

2. 整合式多站点多物种模型（可选） #

若希望同时估计所有站点的参数（比如检验站点间B矩阵的差异是否显著），可以构建整合模型。核心思路是将每个站点的物种状态作为独立的状态组，状态矩阵B为分块对角矩阵（每个块对应一个站点的物种相互作用矩阵）。

关键实现步骤：

• 数据整理：将所有站点的物种观测按行堆叠（如现有代码中的y矩阵）
• Z矩阵：构建分块对角矩阵，每个块对应一个站点的观测-状态映射（即每个站点的Z是对角矩阵）
• B矩阵：构建分块对角矩阵，每个块是对应站点的物种相互作用矩阵，块间元素为0（不同站点的物种状态无相互作用）
• U矩阵：每个站点的每个物种对应独立的截距项

2站点3物种的整合模型代码示例：

n_sites <- 2
n_spp_per_site <- 3
total_spp <- n_sites * n_spp_per_site

# 构建分块对角Z矩阵
z_model <- matrix(0, nrow = total_spp, ncol = total_spp)
for (s in 1:n_sites) {
  start_idx <- (s-1)*n_spp_per_site + 1
  end_idx <- s*n_spp_per_site
  z_model[start_idx:end_idx, start_idx:end_idx] <- diag(n_spp_per_site)
}

# 构建分块对角B矩阵（每个块对应站点的物种相互作用）
b_model <- matrix(0, nrow = total_spp, ncol = total_spp)
for (s in 1:n_sites) {
  start_idx <- (s-1)*n_spp_per_site + 1
  end_idx <- s*n_spp_per_site
  block_labels <- paste0("site", s, "_b_", rep(1:n_spp_per_site, each = n_spp_per_site), 1:n_spp_per_site)
  b_model[start_idx:end_idx, start_idx:end_idx] <- matrix(block_labels, nrow = n_spp_per_site)
}

# 构建U矩阵（每个站点的每个物种对应独立截距）
u_model <- matrix(paste0("site", rep(1:n_sites, each = n_spp_per_site), "_u_", rep(1:n_spp_per_site, n_sites)),
                  nrow = total_spp, ncol = 1)

# 定义整合模型
integrated_model <- list(
  Z = z_model,
  B = b_model,
  U = u_model,
  Q = "diagonal and unequal",  # 每个状态的过程误差独立
  R = "diagonal and equal",    # 所有观测的测量误差相同（可按需调整）
  tinitx = 1
)

# 运行模型
integrated_out <- MARSS(y, model = integrated_model, method = "BFGS")

# 提取每个站点的B矩阵
site_B_list <- list()
for (s in 1:n_sites) {
  start_idx <- (s-1)*n_spp_per_site + 1
  end_idx <- s*n_spp_per_site
  site_B_list[[paste0("site", s)]] <- coef(integrated_out, type = "matrix")$B[start_idx:end_idx, start_idx:end_idx]
}

现有代码的小修正 #

• 单物种双站点模型中，A和U矩阵使用了重复的参数标签"b"，会导致参数混淆，建议给不同站点的参数设置唯一标签（比如u.model1[1:2,1] <- c("site1_u", "site2_u")）
• 单站点多物种模型中，全参数化的B矩阵适合数据量充足的情况，若数据量较小，可约束B矩阵（比如只估计对角项或关键相互作用）避免过拟合

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