我在我的代码中模拟了一些数据，并希望使用lme4：：nlmer拟合非线性混合模型，但我一直得到以下错误：

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

字符串
我的代码是：

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

型
有人有解决这个问题的办法吗？
我也试过：

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)
)

型
在那里我得到了

Error: step factor reduced below 0.001 without reducing pwrss

型
如果我尝试调试它并“跳过”给出此错误的行，则会得到另一个错误。
为了使模型更简单，我尝试移除两个随机效应，多次更改初始值，并尝试在“nlmerControl”中使用其他优化器。此外，当我尝试使用nlme：：nlme拟合非线性混合模型时，它工作得很好。但是，我想比较两个函数的结果。