使用lme4::nlmer拟合非线性混合模型遇正定矩阵错误求助

问题背景 #

模拟数据后使用lme4::nlmer拟合非线性混合模型时持续报错，nlme::nlme可正常拟合，需解决nlmer的拟合问题以对比两者结果。

报错信息 #

第一种代码执行时报错：

Error in devfun(rho$pp$theta) : Downdated VtV is not positive definite

第二种代码执行时报错：

Error: step factor reduced below 0.001 without reducing pwrss

已尝试的方法 #

• 简化模型，移除部分随机效应
• 更换固定效应初始值
• 调整nlmerControl中的优化器类型与迭代次数

尝试过的代码 #

代码1 #

library(nlme)
library(ggplot2)
library(lme4)
library(gridExtra)
library(MASS)

nparams <- c("delta1", "delta2", "delta3", "Time")
fpl1 <- deriv(as.formula(paste0("~((delta1) * exp(Time - (delta2) )) / (1 + exp((Time - (delta2)) / (delta3)))")), c("delta1", "delta2", "delta3"), function.arg = nparams, envir=environment())

createDataset <- function(nGroups, totalDaysMeasured, d1, d2, d3) {
  residualSD <- 0.939921
  nIndividuals <- 600
  D <- rbind(c(51.88633, -1.0498620, -0.05460614), c(-1.0498620, 15.8465000, -0.04587260), c(-0.05460614, -0.04587260, 0.01362791))
  random_effects <- mvrnorm(n = nIndividuals, mu = rep(0, 3), Sigma = D)
  nDays <- totalDaysMeasured
  groups <- rep(paste0("g",1:nGroups),each=3)
  
  days = seq(0, 25, length.out = totalDaysMeasured)
  data <- data.frame(Day = rep(days,each=nIndividuals),
                     Group = rep(groups,nDays),
                     Individual = rep(1:nIndividuals,nDays))
  
  asymp.RNoiseByIndividual1 <- random_effects[, 1]
  asymp.RNoiseByIndividual2 <- random_effects[, 2] 
  asymp.RNoiseByIndividual3 <- random_effects[, 3] 
  data$X <- with(data, fpl1(d1+asymp.RNoiseByIndividual1[Individual], d2+asymp.RNoiseByIndividual2[Individual], d3+asymp.RNoiseByIndividual3[Individual], Day)+rnorm(nrow(data),mean=0,sd=residualSD))
  data
}


totalDaysMeasured <- 25
delta1 <- 5.051236
delta2 <- 13.86703
delta3 <- 0.8486555

data2Groups <- createDataset(1,totalDaysMeasured, delta1, delta2, delta3)
ggplot(data2Groups,aes(x=Day,y=X,colour=Group))+geom_point()


initialValues2Groups <-  c( delta1=delta1,
                            delta2=delta2,
                            delta3=delta3)

nparams <- c("A","C","D")
fpl <- deriv(as.formula(paste0("~(A*exp(x-C))/(1+exp((x-C)/D))")),
             nparams,
             function.arg=c("x",nparams),
             envir=environment())
attr(fpl,"pnames") <- nparams

tmpstr <- deparse(fpl)
L1 <- grep("^ +\\.value +<-",tmpstr)
L2 <- grep("^ +attr\\(.value",tmpstr)
hej <- deparse(nparams)[1]
tmpstr2 <- c(tmpstr[1:L1],
             paste0(".actualArgs <- as.list(match.call()[",
                    hej,"])"),
             tmpstr[(L1+1):(L2-1)],
             "dimnames(.grad) <- list(NULL, .actualArgs)",
             tmpstr[L2:length(tmpstr)])
fpl <- eval(parse(text=tmpstr2),envir = environment())
nlmerString <- paste0("X ~ fpl(delta1, delta2, delta3, Day) ~ (delta1|Individual)+(delta3|Individual)+(delta3|Individual)")
startTimeNLMER <- Sys.time()
nlmerFit <- do.call("nlmer",list(
  as.formula(nlmerString),
  start=initialValues2Groups,
  data=data2Groups, 
  control = nlmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 100000)
  )))

endTimeNLMER <- Sys.time()
#print(paste0("Time required for the ",deparse(substitute(data)),": ",endTimeNLMER-startTimeNLMER))
#nlmerFit

代码2 #

library("R.utils")

# Packages
library(parallel)
library(nlme)
library(ggplot2)
library(Matrix)
library(MASS)
library(nlraa)
library(lme4)


# Initial settings
num_subjects <- 600
num_time_points <- 25
n_sim <- 1000
N_T = num_subjects * num_time_points

# Define fixed effects
delta1 <- 5.051236
delta2 <- 13.86703
delta3 <- 0.8486555
initial_params <- c(delta1 = 5.051236, delta2 = 13.86703, delta3 = 0.8486555)

# Define the upper triangular matrix for the covariance matrix D
D <- rbind(c(51.88633, -1.0498620, -0.05460614), c(-1.0498620, 15.8465000, -0.04587260), c(-0.05460614, -0.04587260, 0.01362791))


# Parallel ---------------------------------------------------------------------
# Model function
nparams <- c("delta1", "delta2", "delta3", "eta1i", "eta2i","eta3i", "Time")
fpl <- deriv(as.formula(paste0("~((delta1+eta1i) * exp(Time - (delta2+eta2i) )) / (1 + exp((Time - (delta2+eta2i)) / (delta3+eta3i)))")), c("eta1i", "eta2i","eta3i"), function.arg = nparams, envir=environment())


nlmerString <- paste0("value ~ fpl(delta1, delta2, delta3, eta1i, eta2i,eta3i, Time) ~ (eta1i|Subject)+(eta2i|Subject)+(eta3i|Subject)")

random_effects <- mvrnorm(n = num_subjects, mu = rep(0, 3), Sigma = D)
eps <- c()

# Create an empty list to store the data
data_list <- list()

# Generate data for each subject
for (subject_id in 1:num_subjects) {
  # Get the random effects for the subject
  eta1i <- random_effects[subject_id, 1]
  eta2i <- random_effects[subject_id, 2]
  eta3i <- random_effects[subject_id, 3]
  eps_i <- rnorm(num_time_points, mean = 0, sd = 0.939921)
  eps <- append(eps, eps_i)
  
  # Create a data frame for the subject
  subject_data <- data.frame(
    Subject = factor(subject_id),
    Time = seq(0, 25, length.out = num_time_points),
    eps = eps_i,
    eta1i = eta1i,
    eta2i = eta2i,
    eta3i = eta3i
  )
  
  # Calculate delta values for the subject
  delta1i <- delta1 + eta1i
  delta2i <- delta2 + eta2i
  delta3i <- delta3 + eta3i
  
  # Generate data based on the specified model
  subject_data$value <- with(subject_data, (delta1i * exp(Time - delta2i)) / (1 + exp((Time - delta2i) / delta3i)) + eps_i)
  
  # Add the subject's data to the list
  data_list[[subject_id]] <- subject_data
}

simulated_data_normal <- do.call("rbind", data_list)

nlmerFit <- do.call("nlmer",list(
  as.formula(nlmerString),
  start = initial_params,
  data=simulated_data_normal)
)

解决方案 #

1. 修正随机效应结构定义 #

第一种代码中随机效应部分重复定义(delta3|Individual)，导致冗余参数与数值不稳定。应改为与数据生成逻辑匹配的完整随机效应结构：

nlmerString <- paste0("X ~ fpl(delta1, delta2, delta3, Day) ~ (delta1 + delta2 + delta3 | Individual)")

2. 简化模型公式，避免自定义导数函数问题 #

自定义deriv函数易出现参数映射错误，可直接在nlmer公式中使用非线性表达式，让lme4自动处理导数计算：

nlmer_form <- X ~ (delta1 * exp(Day - delta2)) / (1 + exp((Day - delta2)/delta3)) ~ (delta1 + delta2 + delta3 | Individual)

3. 利用nlme结果提供精准初始值 #

仅提供固定效应初始值不足以支撑nlmer的数值优化，可从nlme拟合结果中提取随机效应协方差的初始参数：

# 先拟合nlme模型获取初始信息
nlme_fit <- nlme(X ~ (delta1 * exp(Day - delta2)) / (1 + exp((Day - delta2)/delta3)),
                 data = data2Groups,
                 fixed = delta1 + delta2 + delta3 ~ 1,
                 random = delta1 + delta2 + delta3 ~ 1 | Individual,
                 start = initialValues2Groups)
# 提取随机效应协方差的初始theta值
theta_start <- getME(nlme_fit, "theta")

# 用完整初始值拟合nlmer
nlmerFit <- nlmer(nlmer_form,
                  start = c(initialValues2Groups, theta = theta_start),
                  data = data2Groups,
                  control = nlmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 200000)))

4. 修正第二种代码的模型逻辑 #

第二种代码错误地将随机效应作为数据列输入，nlmer会将其视为固定效应。正确做法是移除数据中的eta1i、eta2i、eta3i列，让nlmer自动估计随机效应：

# 生成数据时移除eta列
subject_data <- data.frame(
  Subject = factor(subject_id),
  Time = seq(0, 25, length.out = num_time_points),
  eps = eps_i
)

# 修正模型公式
nlmerString <- paste0("value ~ (delta1 * exp(Time - delta2)) / (1 + exp((Time - delta2)/delta3)) ~ (delta1 + delta2 + delta3 | Subject)")

总结 #

核心问题在于随机效应结构定义错误、自定义导数函数的参数映射偏差，以及初始值信息不足。通过修正随机效应结构、简化模型公式、利用nlme结果提供精准初始值，即可解决nlmer的拟合报错问题。

内容的提问来源于stack exchange，提问作者Sandra Sörensen