在R中用lme4包的lmer拟合含预定义方差成分的线性混合模型

嗨，我来帮你搞定这个问题！你想用自己预先估算好的方差成分（固定效应方差0.0135、随机效应方差0.0045）来拟合模型，而不是让lme4的lmer()自动估计这些参数对吧？

首先得明确：lme4的lmer()函数本身并没有直接支持固定方差成分的参数——它的设计目标就是通过最大似然/限制最大似然来估计这些值。不过别担心，我们有两种可行的解决方案，一种简单直接，另一种是针对lme4的变通方法：

方法1：用nlme包的lme()函数（首推！） #

nlme包的lme()函数专门支持固定方差参数的操作，语法也更直观。先给你针对你的示例数据写好代码：

首先修正你示例里的小问题（你模型里写的data = d，但数据框实际叫invented），然后加载包准备数据：

library(nlme)
set.seed(123) # 加个随机种子，结果能重复

# 生成你的示例数据
f <- rnorm(10, mean=1.4025, sd=0.101)
m <- rnorm(10, mean=1.9025, sd=0.15)
gain <- c(f,m,f,m)
s <- rep(c("f","m"), each=10)
sex <- rep(s, 4)
anim <- as.factor(rep(1:10, 4))
invented <- data.frame(gain,sex,anim)

接下来就是拟合固定方差的模型：

# 固定随机效应方差为0.0045，残差方差（对应你说的固定效应方差）为0.0135
fixed_model <- lme(
  gain ~ sex,
  random = list(anim = pdDiag(~1)), # 指定单个随机截距的结构
  data = invented,
  # 设置初始参数：直接把我们要固定的方差值放进去
  start = list(
    random = list(anim = 0.0045), # 随机效应方差
    sigma = sqrt(0.0135) # 残差标准差（方差开平方）
  ),
  # 禁止模型迭代优化，直接用我们指定的参数
  control = lmeControl(niterEM = 0, msMaxIter = 0)
)

# 查看结果
summary(fixed_model)

这里的关键点：

• pdDiag(~1)告诉模型我们的随机效应是独立的截距，对应单个方差参数；
• start参数直接指定我们要固定的方差值；
• lmeControl(niterEM = 0, msMaxIter = 0)让模型跳过优化步骤，完全使用我们给定的方差参数，只估计固定效应系数。

方法2：在lme4里用自定义优化逻辑固定参数 #

如果你一定要用lme4的lmer()，那得绕个弯——因为它没有直接固定参数的选项，我们需要手动控制优化过程，把方差参数锁死，只优化固定效应。

思路是先拟合一个普通模型，搞清楚参数结构，然后自定义目标函数，让优化器只调整固定效应，方差参数保持我们指定的值：

library(lme4)

# 先拟合一个普通模型，获取参数初始值
init_lmer <- lmer(gain ~ sex + (1|anim), data = invented)
fixef_init <- getME(init_lmer, "beta") # 固定效应初始值
# 把随机效应方差0.0045转换成lme4用的theta参数（Cholesky变换后的对数形式）
theta_fixed <- log(sqrt(0.0045))
sigma_fixed <- sqrt(0.0135) # 残差标准差

# 定义自定义目标函数：只优化固定效应，方差参数固定
fixed_obj <- function(par) {
  # par是待优化的固定效应参数
  temp_mod <- update(init_lmer, 
                     start = list(beta = par, theta = theta_fixed, sigma = sigma_fixed))
  -logLik(temp_mod) # 返回负对数似然，供优化器最小化
}

# 优化固定效应
opt_result <- optim(fixef_init, fixed_obj, method = "BFGS")
final_fixef <- opt_result$par

# 构建最终模型：用固定的方差参数和优化后的固定效应
final_lmer <- lmer(gain ~ sex + (1|anim), data = invented,
                   start = list(beta = final_fixef, theta = theta_fixed, sigma = sigma_fixed),
                   control = lmerControl(optCtrl = list(maxfun = 0))) # 不迭代，直接用指定参数

summary(final_lmer)

这个方法需要你对lme4的参数逻辑有点了解，比如theta是随机效应方差的Cholesky变换后的对数参数，所以要做相应转换。

一些重要提醒 #

• 如果你只是想固定方差成分，强烈建议用nlme的lme()，操作简单，不容易出错；
• 固定方差参数属于非常规操作，一定要确保你有充分的理由这么做（比如方差值来自已发表的研究、理论推导，或者你需要做敏感性分析），不要随便固定参数哦！

内容的提问来源于stack exchange，提问作者T. Grüter