нелинейная регрессия со случайным эффектом и lsoda

Question

Я столкнулся с проблемой, которую не могу решить. Я хотел бы использовать nlmeили nlmODEвыполнить нелинейную регрессию со случайным эффектом, используя в качестве модели решение дифференциального уравнения второго порядка с фиксированными коэффициентами (затухающий осциллятор).

Мне удается использовать nlmeс простыми моделями, но кажется, что использование deSolveдля генерации решения дифференциального уравнения вызывает проблемы. Ниже пример и проблемы, с которыми я столкнулся.

Данные и функции

Вот функция для генерации решения дифференциального уравнения с использованием deSolve:

library(deSolve)
ODE2_nls <- function(t, y, parms) {
  S1 <- y[1]
  dS1 <- y[2]
  dS2 <- dS1
  dS1 <- - parms["esp2omega"]*dS1  - parms["omega2"]*S1 + parms["omega2"]*parms["yeq"]
  res <- c(dS2,dS1)
  list(res)}

solution_analy_ODE2 = function(omega2,esp2omega,time,y0,v0,yeq){
  parms  <- c(esp2omega = esp2omega,
              omega2 = omega2,
              yeq = yeq)
  xstart = c(S1 =  y0, dS1 = v0)
  out <-  lsoda(xstart, time, ODE2_nls, parms)
  return(out[,2])
}

Я могу сгенерировать решение для заданного периода и коэффициента затухания, как, например, здесь период 20 и небольшое затухание 0,2:

# small example:
time <- 1:100
period <- 20 # period of oscillation
amort_factor <- 0.2
omega <- 2*pi/period # agular frequency
oscil <- solution_analy_ODE2(omega^2,amort_factor*2*omega,time,1,0,0)
plot(time,oscil)

Теперь я создаю панель из 10 человек со случайной начальной фазой (то есть с другой начальной позицией и скоростью). Цель состоит в том, чтобы выполнить нелинейную регрессию со случайным влиянием на начальные значения.

library(data.table)
# generate panel
Npoint <- 100 # number of time poitns
Nindiv <- 10 # number of individuals
period <- 20 # period of oscillation
amort_factor <- 0.2
omega <- 2*pi/period # agular frequency
# random phase
phase <- sample(seq(0,2*pi,0.01),Nindiv)
# simu data:
data_simu <- data.table(time = rep(1:Npoint,Nindiv), ID = rep(1:Nindiv,each = Npoint))

# signal generation
data_simu[,signal := solution_analy_ODE2(omega2 = omega^2,
                                         esp2omega = 2*0.2*omega,
                                         time = time,
                                         y0 = sin(phase[.GRP]),
                                         v0 = omega*cos(phase[.GRP]),
                                         yeq = 0)+ 
            rnorm(.N,0,0.02),by = ID]

Если мы посмотрим, у нас есть правильный набор данных:

library(ggplot2)
ggplot(data_simu,aes(time,signal,color = ID))+
  geom_line()+
  facet_wrap(~ID)

Проблемы

Использование nlme

Используя nlmeаналогичный синтаксис, работая над более простыми примерами (нелинейные функции, не использующие deSolve), я попробовал:

fit <- nlme(model = signal ~ solution_analy_ODE2(esp2omega,omega2,time,y0,v0,yeq), 
     data = data_simu,
     fixed = esp2omega + omega2 + y0 + v0 + yeq ~ 1,
     random = y0 ~ 1 ,
     groups = ~ ID, 
     start = c(esp2omega = 0.08, 
               omega2 = 0.04,
               yeq = 0,
               y0 = 1,
               v0 = 0))

Я получаю:

Ошибка в checkFunc (Func2, times, y, rho): количество производных, возвращаемых функцией func () (2), должно равняться длине вектора начальных условий (2000)

Отслеживание:

12. stop(paste("The number of derivatives returned by func() (", length(tmp[[1]]), ") must equal the length of the initial conditions vector (", length(y), ")", sep = ""))
11. checkFunc(Func2, times, y, rho)
10. lsoda(xstart, time, ODE2_nls, parms)
9. solution_analy_ODE2(omega2, esp2omega, time, y0, v0, yeq)
.
.

Похоже, я nlmeпытаюсь передать вектор начального условия solution_analy_ODE2и вызывает ошибку в checkFuncfrom lasoda.

Я пробовал использовать nlsList:

test <- nlsList(model = signal ~ solution_analy_ODE2(omega2,esp2omega,time,y0,v0,yeq) | ID, 
        data = data_simu, 
        start = list(esp2omega = 0.08, omega2 = 0.04,yeq = 0,
                     y0 = 1,v0 = 0),
        control = list(maxiter=150, warnOnly=T,minFactor = 1e-10), 
        na.action = na.fail, pool = TRUE)
head(test)

Call:
  Model: signal ~ solution_analy_ODE2(omega2, esp2omega, time, y0, v0, yeq) | ID 
   Data: data_simu 

Coefficients:
   esp2omega     omega2           yeq         y0          v0
1  0.1190764 0.09696076  0.0007577956 -0.1049423  0.30234654
2  0.1238936 0.09827158 -0.0003463023  0.9837386  0.04773775
3  0.1280399 0.09853310 -0.0004908579  0.6051663  0.25216134
4  0.1254053 0.09917855  0.0001922963 -0.5484005 -0.25972829
5  0.1249473 0.09884761  0.0017730823  0.7041049  0.22066652
6  0.1275408 0.09966155 -0.0017522320  0.8349450  0.17596648

Мы видим, что нелинейная аппроксимация хорошо работает для отдельных сигналов. Теперь, если я хочу выполнить регрессию набора данных со случайными эффектами, синтаксис должен быть:

fit <- nlme(test, 
     random = y0 ~ 1 ,
     groups = ~ ID, 
     start = c(esp2omega = 0.08, 
               omega2 = 0.04,
               yeq = 0,
               y0 = 1,
               v0 = 0))

Но я получаю точно такое же сообщение об ошибке.

Затем я попытался использовать nlmODE, следуя комментарию Бне Болкера по аналогичному вопросу, который я задал несколько лет назад

используя nlmODE

library(nlmeODE)
datas_grouped <- groupedData( signal ~ time | ID, data = data_simu, 
                              labels = list (x = "time", y = "signal"), 
                              units = list(x ="arbitrary", y = "arbitrary"))

modelODE <- list( DiffEq = list(dS2dt = ~ S1,
                                dS1dt = ~ -esp2omega*S1  - omega2*S2 + omega2*yeq),
                  ObsEq = list(yc = ~ S2),
                  States = c("S1","S2"),
                  Parms = c("esp2omega","omega2","yeq","ID"), 
                  Init = c(y0 = 0,v0 = 0))

resnlmeode = nlmeODE(modelODE, datas_grouped) 
assign("resnlmeode", resnlmeode, envir = .GlobalEnv)
#Fitting with nlme the resulting function
model <- nlme(signal ~ resnlmeode(esp2omega,omega2,yeq,time,ID), 
              data = datas_grouped, 
              fixed = esp2omega + omega2 + yeq + y0 + v0  ~ 1, 
              random = y0 + v0 ~1,
              start = c(esp2omega = 0.08, 
                        omega2 = 0.04,
                        yeq = 0,
                        y0 = 0,
                        v0 = 0)) # 

Я получаю сообщение об ошибке:

Ошибка в resnlmeode (esp2omega, omega2, yeq, time, ID): объект yhat не найден

Здесь я не понимаю, откуда взялась ошибка, и как ее решить.

Вопросов

• Вы можете воспроизвести проблему?
• Есть ли у кого-нибудь идея решить эту проблему, используя либо nlmeили nlmODE?
• Если нет, есть ли решение с использованием другого пакета? Я видел nlmixr(https://cran.r-project.org/web/packages/nlmixr/index.html), но я этого не знаю, установка сложная и недавно была удалена из CRAN

---

Правки

@tpetzoldt предложил хороший способ отладки nlmeповедения, и это меня очень удивило. Вот рабочий пример с нелинейной функцией, где я генерирую набор из 5 человек со случайным параметром, варьирующимся между людьми:

reg_fun = function(time,b,A,y0){
  cat("time : ",length(time)," b :",length(b)," A : ",length(A)," y0: ",length(y0),"\n")
  out <- A*exp(-b*time)+(y0-1)
  cat("out : ",length(out),"\n")
  tmp <- cbind(b,A,y0,time,out)
  cat(apply(tmp,1,function(x) paste(paste(x,collapse = " "),"\n")),"\n")
  return(out)
}

time <- 0:10*10
ramdom_y0 <- sample(seq(0,1,0.01),10)
Nid <- 5
data_simu <- 
data.table(time = rep(time,Nid),
           ID = rep(LETTERS[1:Nid],each = length(time)) )[,signal := reg_fun(time,0.02,2,ramdom_y0[.GRP]) + rnorm(.N,0,0.1),by = ID]

Кошки в функции дают здесь:

time :  11  b : 1  A :  1  y0:  1 
out :  11 
0.02 2 0.64 0 1.64 
 0.02 2 0.64 10 1.27746150615596 
 0.02 2 0.64 20 0.980640092071279 
 0.02 2 0.64 30 0.737623272188053 
 0.02 2 0.64 40 0.538657928234443 
 0.02 2 0.64 50 0.375758882342885 
 0.02 2 0.64 60 0.242388423824404 
 0.02 2 0.64 70 0.133193927883213 
 0.02 2 0.64 80 0.0437930359893108 
 0.02 2 0.64 90 -0.0294022235568269 
 0.02 2 0.64 100 -0.0893294335267746
.
.
.

Сейчас занимаюсь nlme:

nlme(model = signal ~ reg_fun(time,b,A,y0), 
     data = data_simu,
     fixed = b + A + y0 ~ 1,
     random = y0 ~ 1 ,
     groups = ~ ID, 
     start = c(b = 0.03, A = 1,y0 = 0))

Я получил:

time :  55  b : 55  A :  55  y0:  55 
out :  55 
0.03 1 0 0 0 
 0.03 1 0 10 -0.259181779318282 
 0.03 1 0 20 -0.451188363905974 
 0.03 1 0 30 -0.593430340259401 
 0.03 1 0 40 -0.698805788087798 
 0.03 1 0 50 -0.77686983985157 
 0.03 1 0 60 -0.834701111778413 
 0.03 1 0 70 -0.877543571747018 
 0.03 1 0 80 -0.909282046710588 
 0.03 1 0 90 -0.93279448726025 
 0.03 1 0 100 -0.950212931632136 
 0.03 1 0 0 0 
 0.03 1 0 10 -0.259181779318282 
 0.03 1 0 20 -0.451188363905974 
 0.03 1 0 30 -0.593430340259401 
 0.03 1 0 40 -0.698805788087798 
 0.03 1 0 50 -0.77686983985157 
 0.03 1 0 60 -0.834701111778413 
 0.03 1 0 70 -0.877543571747018 
 0.03 1 0 80 -0.909282046710588 
 0.03 1 0 90 -0.93279448726025 
 0.03 1 0 100 -0.950212931632136 
 0.03 1 0 0 0 
 0.03 1 0 10 -0.259181779318282 
 0.03 1 0 20 -0.451188363905974 
 0.03 1 0 30 -0.593430340259401 
 0.03 1 0 40 -0.698805788087798 
 0.03 1 0 50 -0.77686983985157 
 0.03 1 0 60 -0.834701111778413 
 0.03 1 0 70 -0.877543571747018 
 0.03 1 0 80 -0.909282046710588 
 0.03 1 0 90 -0.93279448726025 
 0.03 1 0 100 -0.950212931632136 
 0.03 1 0 0 0 
 0.03 1 0 10 -0.259181779318282 
 0.03 1 0 20 -0.451188363905974 
 0.03 1 0 30 -0.593430340259401 
 0.03 1 0 40 -0.698805788087798 
 0.03 1 0 50 -0.77686983985157 
 0.03 1 0 60 -0.834701111778413 
 0.03 1 0 70 -0.877543571747018 
 0.03 1 0 80 -0.909282046710588 
 0.03 1 0 90 -0.93279448726025 
 0.03 1 0 100 -0.950212931632136 
 0.03 1 0 0 0 
 0.03 1 0 10 -0.259181779318282 
 0.03 1 0 20 -0.451188363905974 
 0.03 1 0 30 -0.593430340259401 
 0.03 1 0 40 -0.698805788087798 
 0.03 1 0 50 -0.77686983985157 
 0.03 1 0 60 -0.834701111778413 
 0.03 1 0 70 -0.877543571747018 
 0.03 1 0 80 -0.909282046710588 
 0.03 1 0 90 -0.93279448726025 
 0.03 1 0 100 -0.950212931632136 
 
time :  55  b : 55  A :  55  y0:  55 
out :  55 
0.03 1 0 0 0 
 0.03 1 0 10 -0.259181779318282 
 0.03 1 0 20 -0.451188363905974 
 0.03 1 0 30 -0.593430340259401 
 0.03 1 0 40 -0.698805788087798 
 0.03 1 0 50 -0.77686983985157 
 0.03 1 0 60 -0.834701111778413 
 0.03 1 0 70 -0.877543571747018 
 0.03 1 0 80 -0.909282046710588 
 0.03 1 0 90 -0.93279448726025 
 0.03 1 0 100 -0.950212931632136 
 0.03 1 0 0 0 
 0.03 1 0 10 -0.259181779318282 
 0.03 1 0 20 -0.451188363905974 
 0.03 1 0 30 -0.593430340259401 
 0.03 1 0 40 -0.698805788087798 
 0.03 1 0 50 -0.77686983985157 
 0.03 1 0 60 -0.834701111778413 
 0.03 1 0 70 -0.877543571747018 
 0.03 1 0 80 -0.909282046710588 
 0.03 1 0 90 -0.93279448726025 
 0.03 1 0 100 -0.950212931632136 
 0.03 1 0 0 0 
 0.03 1 0 10 -0.259181779318282 
 0.03 1 0 20 -0.451188363905974 
 0.03 1 0 30 -0.593430340259401 
 0.03 1 0 40 -0.698805788087798 
 0.03 1 0 50 -0.77686983985157 
 0.03 1 0 60 -0.834701111778413 
 0.03 1 0 70 -0.877543571747018 
 0.03 1 0 80 -0.909282046710588 
 0.03 1 0 90 -0.93279448726025 
 0.03 1 0 100 -0.950212931632136 
...

Таким образом, nlmeпривязывает 5 раз (число отдельных) к вектору времени и передает его функции, при этом параметры повторяются одинаковое количество раз. Что, конечно, несовместимо с тем, как lsodaработает моя функция.

Answer 1

Кажется, что модель ode вызывается с неправильным аргументом, поэтому она получает вектор с 2000 переменных состояния вместо 2. Попробуйте следующее, чтобы увидеть проблему:

ODE2_nls <- function(t, y, parms) {
  cat(length(y),"\n") # <----
  S1 <- y[1]
  dS1 <- y[2]
  dS2 <- dS1
  dS1 <- - parms["esp2omega"]*dS1  - parms["omega2"]*S1 + parms["omega2"]*parms["yeq"]
  res <- c(dS2,dS1)
  list(res)
}

Изменить : я думаю, что аналитическая функция сработала, потому что она векторизована, поэтому вы можете попытаться векторизовать функцию ode, либо перебирая модель ode, либо (лучше) внутренне используя векторы в качестве переменных состояния. Поскольку odeрешение систем с несколькими сотнями тысяч уравнений является быстрым, 2000 должно быть возможным.

Я предполагаю, что и состояния, и параметры nlmeпередаются как векторы. Переменная состояния модели оды тогда является «длинным» вектором, параметры могут быть реализованы в виде списка.

Вот пример (отредактированный, теперь с параметрами в виде списка):

ODE2_nls <- function(t, y, parms) {
  #cat(length(y),"\n")
  #cat(length(parms$omega2)) ndx <- seq(1, 2*N-1, 2) S1 <- y[ndx] dS1 <- y[ndx + 1] dS2 <- dS1 dS1 <- - parms$esp2omega * dS1  - parms$omega2 * S1 + parms$omega2 * parms$yeq
  res <- c(dS2, dS1)
  list(res)
}

solution_analy_ODE2 = function(omega2, esp2omega, time, y0, v0, yeq){
  parms  <- list(esp2omega = esp2omega, omega2 = omega2, yeq = yeq)
  xstart = c(S1 =  y0, dS1 = v0)
  out <-  ode(xstart, time, ODE2_nls, parms, atol=1e-4, rtol=1e-4, method="ode45")
  return(out[,2])
}

Затем задайте (или вычислите) количество уравнений, например, N <- 1соотв. N <-1000перед звонками.

Модель проходит этот путь, прежде чем запускать числовые проблемы, но это уже другая история ...

Затем вы можете попытаться использовать другой решатель од (например vode), установить atolи rtolна более низкие значения, настроить nmleпараметры оптимизации, использовать ограничения блока ... и так далее, как обычно при нелинейной оптимизации.

Answer 2

Я нашел способ взлома решения nlme: как показано в моем редактировании, проблема возникает из-за того, что nlmeвектор NindividualxNpoints передается нелинейной функции, предполагая, что функция связывает для каждой временной точки значение. Но lsodaне делайте этого, поскольку он интегрирует уравнение по времени (т.е. ему нужно все время до определенного момента времени, чтобы получить значение).

Мое решение состоит в разложении параметров, которые nlmeпередаются моей функции, вычислении и воссоздании вектора:

detect_id <- function(vec){
  tmp <- c(0,diff(vec))
  out <- tmp
  out <- NA
  out[tmp < 0] <- 1:sum(tmp < 0)
  out <- na.locf(out,na.rm = F)
  rleid(out)
}

detect_id разложить временной вектор на единичный идентификатор временных векторов:

detect_id(rep(1:10,3))
[1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3

Затем функция выполняет цикл числового интегрирования по каждому индивидууму и связывает полученные векторы вместе:

solution_analy_ODE2_modif = function(omega2,esp2omega,time,y0,v0,yeq){
  tmp <- detect_id(time)
  
  out <- lapply(unique(tmp),function(i){
    idxs <- which(tmp == i)
    parms  <- c(esp2omega = esp2omega[idxs][1],
                omega2 = omega2[idxs][1],
                yeq = yeq[idxs][1])
    
    xstart = c(S1 =  y0[idxs][1], dS1 = v0[idxs][1])
    out_tmp <-  lsoda(xstart, time[idxs], ODE2_nls, parms)
    out_tmp[,2]
  }) %>% unlist()
  
  return(out)
}

Я провожу тест, в котором передаю вектор, похожий на то, что nlmeпередается функции:

omega2vec <- rep(0.1,30)
eps2omegavec <- rep(0.1,30)
timevec <- rep(1:10,3)
y0vec <- rep(1,30)
v0vec <- rep(0,30)
yeqvec = rep(0,30)
solution_analy_ODE2_modif(omega2 = omega2vec,
                          esp2omega = eps2omegavec,
                          time = timevec,
                          y0 = y0vec,
                          v0 = v0vec,
                          yeq = yeqvec)
 [1]  1.0000000  0.9520263  0.8187691  0.6209244  0.3833110  0.1321355 -0.1076071 -0.3143798
 [9] -0.4718058 -0.5697255  1.0000000  0.9520263  0.8187691  0.6209244  0.3833110  0.1321355
[17] -0.1076071 -0.3143798 -0.4718058 -0.5697255  1.0000000  0.9520263  0.8187691  0.6209244
[25]  0.3833110  0.1321355 -0.1076071 -0.3143798 -0.4718058 -0.5697255

Оно работает. Это не будет работать с методом @tpetzoldt, потому что вектор времени переходит от 10 до 0, что может вызвать проблемы интеграции. Здесь мне действительно нужно взломать способ nlnmeработы. Сейчас :

fit <- nlme(model = signal ~ solution_analy_ODE2_modif (esp2omega,omega2,time,y0,v0,yeq), 
     data = data_simu,
     fixed = esp2omega + omega2 + y0 + v0 + yeq ~ 1,
     random = y0 ~ 1 ,
     groups = ~ ID, 
     start = c(esp2omega = 0.5, 
     omega2 = 0.5,
     yeq = 0,
     y0 = 1,
     v0 = 1))

работает как шарм

summary(fit)


Nonlinear mixed-effects model fit by maximum likelihood
  Model: signal ~ solution_analy_ODE2_modif(omega2, esp2omega, time, y0,      v0, yeq) 
 Data: data_simu 
        AIC       BIC   logLik
  -597.4215 -567.7366 307.7107

Random effects:
 Formula: list(y0 ~ 1, v0 ~ 1)
 Level: ID
 Structure: General positive-definite, Log-Cholesky parametrization
         StdDev     Corr  
y0       0.61713329 y0    
v0       0.67815548 -0.269
Residual 0.03859165       

Fixed effects: esp2omega + omega2 + y0 + v0 + yeq ~ 1 
              Value  Std.Error  DF   t-value p-value
esp2omega 0.4113068 0.00866821 186  47.45002  0.0000
omega2    1.0916444 0.00923958 186 118.14876  0.0000
y0        0.3848382 0.19788896 186   1.94472  0.0533
v0        0.1892775 0.21762610 186   0.86974  0.3856
yeq       0.0000146 0.00283328 186   0.00515  0.9959
 Correlation: 
       esp2mg omega2 y0     v0    
omega2  0.224                     
y0      0.011 -0.008              
v0      0.005  0.030 -0.269       
yeq    -0.091 -0.046  0.009 -0.009

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.2692477 -0.6122453  0.1149902  0.6460419  3.2890201 

Number of Observations: 200
Number of Groups: 10