随机效应模型中optim函数极大似然估计不收敛问题求助

我有一份包含不同季节重复测量的植物生长纵向数据集，需要估计以下生长方程的参数：

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

其中：

• Lt：测量时的株高
• Lt_1：前一次测量的株高
• td：两次测量的时间间隔

需要估计的参数为k和p，其中k随季节变化（每个季节对应一个特定k值），同时模型需加入个体变异作为随机效应I。

我用optim函数编写了参数估计代码，但运行耗时极长且无法收敛。实际数据中个体数量更多（n=100），计算难度会进一步提升，现寻求高效且有效的参数估计方法建议或替代方案。

模拟数据 #

# generate dummy data
n = 10
id = as.factor(1:n)
season = c("s1","s2","s3","s4")
L0 = rnorm(n,40,5)
td = rnorm(n, 180, 10)

growth <- function(td, Lt_1){
  return(function(k,p){
    150*((1 + ((150/Lt_1)^(1/p)-1)*exp(-k*td/365))^(-p))
  })
}
L1 = growth(Lt_1 = L0, td = td)(p = 1.2, k = 0.55)
L2 = growth(Lt_1 = L1, td = td)(p = 1.2, k = 1.09)
L3 = growth(Lt_1 = L2, td = td)(p = 1.2, k = 0.62)
L4 = growth(Lt_1 = L3, td = td)(p = 1.2, k = 0.97)

df <- data.frame(id = id, L0 = L0, L1 = L1, L2 = L2, L3 = L3, L4 = L4) %>% 
  pivot_longer(cols = L0:L4, names_to = "Ltype", values_to = "Lobs") %>% 
  mutate(Lobs = round(Lobs)) %>% 
  group_by(id) %>% 
  mutate(Lt_1 = lag(Lobs)) %>% 
  filter(!is.na(Lt_1)) %>% 
  mutate(td = round(rnorm(4, 180, 10))) %>% 
  mutate(season = season) %>% 
  ungroup() %>% 
  select(id, season, td, Lt_1, Lobs)

head(df)
# A tibble: 6 x 5
  id    season    td  Lt_1  Lobs
  <fct> <chr>  <dbl> <dbl> <dbl>
1 1     s1       183    36    44
2 1     s2       190    44    65
3 1     s3       193    65    77
4 1     s4       182    77    96
5 2     s1       191    43    53
6 2     s2       177    53    75

原负对数似然计算代码 #

# define the function to calculate the sum of negative log likelihood
len_id = length(unique(df$id))

sum_negloglike <- function(par){
  # parameters to estimate
  p = par[1]
  sigma = par[2]
  k <- par[3:6]
  I = par[seq(7, 6+length(unique(df$id)), by = 1)]
  
  # individuals and seasons to loop
  indv <- unique(df$id)
  season <- unique(df$season)
  
  # define the dataframe
  new_df <- data.frame()
  
  # assign k and I
  ## for individual i
  for(i in 1:length(indv)){
    indv_i = indv[i]
    I_i = I[which(indv == indv_i)]
    
    df_i = df %>% 
      filter(id == indv_i) %>% 
      mutate(I = I_i)
    
    # for season j in indvidual i
    for (j in 1:length(season)){
      season_j = season[j]
      k_j = k[which(season == season_j)]
    
      df_i_j = df_i %>% 
        filter(season == season_j) %>% 
        mutate(k = k_j) %>% 
        # estimate the growth
        mutate(Lpred = growth(td, Lt_1)(k+I, p))
    
      # add it to the dataframe
      new_df <- bind_rows(new_df, df_i_j)
    }
  }
  
  # calculate the sum of negative log-likelihood
  Lobs <- new_df$Lobs
  Lpred <- new_df$Lpred
  negloglike_sum = sum_negloglike(Lobs)(Lpred, sigma)
  
  return(negloglike_sum)
}

原极大似然估计代码 #

# set initial parameters
p = 1.2
sigma = 1.88
k <- c(0.55, 1.09, 0.62, 0.97)
I = rnorm(len_id, 0, 0.01)
par = c(p, sigma, k, I)

# calculate the sum of likelihood with initial parameter values
sum_negloglike(par) 

# estimate the parameters
optim(par, sum_negloglike)

优化方案与替代方法 #

1. 重构似然计算代码，消除冗余循环 #

原代码的嵌套循环+反复绑定数据框是核心性能瓶颈，改为向量化计算可大幅提升速度：

sum_negloglike_optimized <- function(par) {
  p <- par[1]
  sigma <- par[2]
  # 给k值命名，匹配季节
  k_vals <- par[3:6]
  names(k_vals) <- unique(df$season)
  # 给I值命名，匹配个体
  I_vals <- par[7:(6 + length(unique(df$id)))]
  names(I_vals) <- unique(df$id)
  
  # 直接匹配每个观测对应的k和I
  df$k <- k_vals[df$season]
  df$I <- I_vals[df$id]
  
  # 向量化计算预测值
  Lt_1 <- df$Lt_1
  td <- df$td
  k_plus_I <- df$k + df$I
  term1 <- (150 / Lt_1)^(1/p) - 1
  term2 <- exp(-k_plus_I * td / 365)
  Lpred <- 150 * ((1 + term1 * term2)^(-p))
  
  # 计算负对数似然（假设误差服从正态分布）
  negloglike <- sum(dnorm(df$Lobs, mean = Lpred, sd = sigma, log = TRUE))
  return(-negloglike)
}

2. 使用专门的混合模型框架 #

针对带随机效应的非线性纵向数据，nlme或lme4包的优化算法比手动optim更稳定高效：

library(nlme)
# 确保数据按个体排序
df <- df %>% arrange(id)

# 定义生长模型
growth_formula <- function(Lt_1, td, p, k_s1, k_s2, k_s3, k_s4, I) {
  # 匹配季节对应的k值
  k <- case_when(
    df$season == "s1" ~ k_s1,
    df$season == "s2" ~ k_s2,
    df$season == "s3" ~ k_s3,
    df$season == "s4" ~ k_s4
  )
  term1 <- (150 / Lt_1)^(1/p) - 1
  term2 <- exp(-(k + I) * td / 365)
  150 * ((1 + term1 * term2)^(-p))
}

# 拟合非线性混合效应模型
fit <- nlme(Lobs ~ growth_formula(Lt_1, td, p, k_s1, k_s2, k_s3, k_s4, I),
            data = df,
            fixed = list(p + k_s1 + k_s2 + k_s3 + k_s4 ~ 1),
            random = pdDiag(I ~ 1), # 个体随机效应
            start = c(p = 1.2, k_s1 = 0.55, k_s2 = 1.09, k_s3 = 0.62, k_s4 = 0.97),
            control = nlmeControl(maxIter = 1000, msMaxIter = 1000))

3. 优化初始值设置 #

• 先固定随机效应，仅优化固定效应（p、各季节k），得到稳定估计后再加入随机效应
• 避免用随机数作为随机效应初始值，可先假设I=0，得到固定效应后再更新

4. 参数约束与模型简化 #

• 对k、p施加非负约束（比如用指数变换k = exp(par)），避免优化过程出现无意义参数值
• 若部分季节k的差异不显著，可合并参数减少优化维度

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