lme4中lmer函数的参数指定与分季节k值估计问题

数据集与模型公式 #

我们有一份植物生长纵向数据集，每行代表单株植物在多季节的测量值，变量包括id（植株ID）、season（季节）、td（测量间隔时长）、Lt（当期株高）、Lt_1（上期株高），数据集构建代码如下：

library(lme4)

df <- data.frame(
  id = c("a", "a", "a", "a", "a", "a", "b", "b", "b", "b", "b", "c", "c", "c", "c", "c", "d", "d", "d", "d", "d", "d", "e", "e", "e", "e", "e", "e"),
  season = c("s1", "s1", "s2", "s3", "s3", "s4", "s1", "s1", "s2", "s3", "s3", "s1", "s1", "s2", "s3", "s3", "s1", "s1", "s2", "s3", "s3", "s4", "s1", "s1", "s2", "s3", "s3", "s4"),
  td = c(59, 95, 210, 134, 62, 169, 60, 96, 208, 134, 62, 60, 94, 210, 133, 63, 59, 99, 206, 134, 62, 169, 59, 122, 182, 134, 62, 171),
  Lt = c(40, 44, 72, 84, 87, 101, 40, 44, 66, 75, 77, 43, 47, 67, 79, 80, 38, 42, 64, 75, 76, 80, 44, 47, 69, 79, 80, 99),
  Lt_1 = c(35.5, 40, 44.5, 72, 84.2, 87, 35, 40.1, 43.8, 65.9, 75.4, 40.8, 43, 47, 66.8, 78.7, 34, 37.5, 41.9, 63.6, 75.2, 76.5, 39.2, 43.5, 46.9, 68.7, 79.4, 80.5)
)

拟采用的非线性生长模型公式为：

formula = Lt ~ 150*((1 + ((150/Lt_1)^(1/p)-1)*exp(-k*td/365))^(-p))

问题1：用lme4::nlmer实现带随机效应的非线性混合模型 #

lmer仅支持线性混合模型，lme4包中对应的非线性混合模型工具是nlmer，可以实现和nlme类似的功能。给参数k添加id水平随机效应的代码如下：

library(lme4)
# 封装非线性生长模型为自定义函数
growth_model <- function(Lt_1, td, k, p) {
  150 * ((1 + ((150 / Lt_1)^(1/p) - 1) * exp(-k * td / 365))^(-p))
}

# 拟合非线性混合模型
model_nlmer <- nlmer(
  Lt ~ growth_model(Lt_1, td, k, p) ~ (k | id),
  data = df,
  start = list(fixed = c(k = 0.2, p = 1))
)

关键说明： #

• 先将非线性模型封装为函数，便于nlmer调用
• 公式中~ (k | id)表示参数k在id组间存在截距型随机效应
• start参数需以list(fixed = c(...))格式指定固定效应的初始值

问题2：分季节估计参数k #

要为每个季节单独估计k的固定效应值，需将season作为分类变量引入模型，为不同季节分配独立的k参数，实现代码如下：

library(dplyr)
# 定义支持分季节k的模型函数
growth_model_season <- function(Lt_1, td, k_s1, k_s2, k_s3, k_s4, p) {
  # 根据季节匹配对应的k值
  k <- case_when(
    season == "s1" ~ k_s1,
    season == "s2" ~ k_s2,
    season == "s3" ~ k_s3,
    season == "s4" ~ k_s4
  )
  150 * ((1 + ((150 / Lt_1)^(1/p) - 1) * exp(-k * td / 365))^(-p))
}

# 拟合分季节k的非线性混合模型（保留id水平的随机效应）
model_season <- nlmer(
  Lt ~ growth_model_season(Lt_1, td, k_s1, k_s2, k_s3, k_s4, p) 
  ~ (0 + k_s1 + k_s2 + k_s3 + k_s4 | id),  # 每个季节的k都带随机效应，可按需调整
  data = df,
  start = list(fixed = c(k_s1 = 0.2, k_s2 = 0.2, k_s3 = 0.2, k_s4 = 0.2, p = 1))
)

关键说明： #

• 用case_when根据season取值匹配对应的k参数
• 若不需要每个季节的k都带随机效应，可修改随机效应部分，比如改为~ (1 | id)仅保留整体截距随机效应
• 初始值需为每个季节的k单独设置，确保模型收敛

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