含随机效应的回归模型方程重构问题求助

我拟合了一个以Age（因子变量）为固定效应、Year（因子变量）为随机效应的混合效应模型，模拟数据与brms拟合代码如下：

set.seed(123)
# Packages
library(brms)
#install.packages("rstan")
library(rstan)

# Data simulation
age <- 10:20
year <- 1990:2000
n1 <-  length(age)
n2 <-  length(year) 
n <- n1*n2
y <- rnorm(n = n,mean = 50,sd = 30 )

beta_0 <- 0.25
beta_age <- 0.25
beta_year <- - 0.85

y <-beta_0+beta_age*age+beta_year*year

# Final Dataframe
d <- data.frame(y=y,
     Age=as.factor(rep(age,n2)),
     Year=as.factor(rep(year,n1))
     )

fit <- brm(y ~ 0+Age+(1|Year), d, iter = 50, warmup = 10, chains = 1)

模型输出的固定效应和随机效应结果如下：

固定效应（fixef(fit)） #

Estimate Est.Error      Q2.5      Q97.5
Age10 -1.5333147 1.3878804 -3.839673  1.0539617
Age11 -0.1285437 0.9321314 -1.221775  1.9525542
Age12 -2.1895797 1.2388774 -4.829345 -1.0598727
Age13  3.7247485 2.4157734  1.408428  8.4137437
Age14 -3.0134557 2.1435395 -8.148383 -1.1134809
Age15 -0.6631218 1.2442192 -2.122722  1.8919087
Age16 -1.2382832 1.6596973 -4.145329  0.3973843
Age17  1.7116087 2.4962776 -1.389026  5.3579521
Age18 -1.1881552 1.3884096 -4.155563  0.2962675
Age19 -0.6924640 1.0659378 -1.558155  1.3804710
Age20 -0.4457888 0.6615983 -1.763389  0.8225215

随机效应（ranef(fit)） #

$Year
, , Intercept

Estimate Est.Error       Q2.5    Q97.5
1990 -1.3702433 1.7680814 -3.2623355 0.885315
1991 -2.2592311 2.3086167 -4.6791703 1.036105
1992  1.7228745 2.5290462 -0.4616315 4.225590
1993 -0.8166271 1.2671842 -2.1690639 1.809690
1994  1.7928488 2.2368656 -0.6837609 4.464356
1995  0.7169620 1.8096562 -0.5032051 3.002480
1996 -0.1112231 0.7797995 -0.7787910 1.614591
1997  2.1406982 2.2829444 -0.2656490 4.370186
1998 -0.3399042 1.0476426 -1.3453246 2.147565
1999  1.7134314 1.9736618 -0.5384686 4.063977
2000 -0.5178907 1.0331019 -1.5312055 1.854528

现在我已得到5年期的Year随机效应预测值Year_F（示例值：c(1.0342, 0.9514, 0.9234, 0.8345, 0.7863)，长度小于原Year的长度），请问如何重构预测方程 y_F ~ 0 + Age + (1 | Year_F)？

要重构含新随机效应Year_F的预测方程，核心是利用已估计的固定效应系数，结合新的随机效应值计算预测值，具体步骤如下：

1. 提取固定效应系数 #

从拟合好的模型中提取Age的固定效应估计值，优先选用Estimate列的数值：

age_coefs <- fixef(fit)[, "Estimate"]

2. 构建预测数据集 #

创建包含所有Age水平和新Year_F的组合数据集，确保每个Age与每个Year_F都能匹配：

# 示例Year_F值
Year_F <- c(1.0342, 0.9514, 0.9234, 0.8345, 0.7863)
# 为Year_F添加命名（可选，方便后续识别）
names(Year_F) <- paste0("Year_F_", 1:length(Year_F))

# 生成全组合的预测数据框
pred_data <- expand.grid(
  Age = levels(d$Age),
  Year_F = names(Year_F)
)
# 匹配对应Year_F的随机效应值
pred_data$Year_F_effect <- Year_F[pred_data$Year_F]

3. 计算预测值y_F #

原模型不含截距项（~0+Age），因此预测值直接等于对应Age的固定效应加上Year_F的随机效应值：

# 匹配Age对应的固定效应系数
pred_data$age_effect <- age_coefs[pred_data$Age]
# 计算最终预测值y_F
pred_data$y_F <- pred_data$age_effect + pred_data$Year_F_effect

4. （可选）封装为可复用函数 #

如果需要重复进行类似预测，可以将上述步骤封装成函数：

predict_yF <- function(fit_model, year_f_effects) {
  age_coefs <- fixef(fit_model)[, "Estimate"]
  # 为新随机效应生成命名
  year_f_names <- paste0("Year_F_", 1:length(year_f_effects))
  names(year_f_effects) <- year_f_names
  
  # 构建预测数据集
  pred_data <- expand.grid(
    Age = levels(fit_model$data$Age),
    Year_F = year_f_names
  )
  pred_data$Year_F_effect <- year_f_effects[pred_data$Year_F]
  pred_data$age_effect <- age_coefs[pred_data$Age]
  # 计算预测值
  pred_data$y_F <- pred_data$age_effect + pred_data$Year_F_effect
  
  return(pred_data)
}

# 函数使用示例
result <- predict_yF(fit, c(1.0342, 0.9514, 0.9234, 0.8345, 0.7863))

关键说明 #

• 原模型无截距项，因此无需额外添加截距值，直接计算固定效应与随机效应的和即可。
• 如果Year_F本身带有年份标签（如c("2021"=1.0342, "2022"=0.9514)），可以直接使用原标签作为Year_F的水平，无需重新命名。
• 若要考虑系数的不确定性，可以结合fixef(fit)的置信区间和Year_F的预测区间，进一步计算y_F的置信区间。

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