rastgele etki ve lsoda ile doğrusal olmayan regresyon

Question

Çözemediğim bir sorunla karşı karşıyayım. Sabit katsayılı (sönümlü bir osilatör) ikinci dereceden diferansiyel denklemin çözümünü model olarak kullanarak rastgele efektli doğrusal olmayan bir regresyon kullanmak nlmeveya nlmODEgerçekleştirmek istiyorum .

nlmeBasit modellerle kullanmayı başarıyorum , ancak görünen o ki deSolvediferansiyel denklemin çözümünü üretmek için kullanımı bir soruna neden oluyor. Aşağıda bir örnek ve karşılaştığım sorunlar.

Veriler ve işlevler

Diferansiyel denklemin çözümünü oluşturmak için aşağıdaki fonksiyon şudur 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])
}

Belirli bir süre ve sönümleme faktörü için bir çözüm üretebilirim, örneğin burada 20'lik bir süre ve 0,2'lik hafif bir sönümleme:

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

Şimdi rastgele bir başlangıç ​​fazına (yani farklı başlangıç ​​pozisyonu ve hız) sahip 10 kişilik bir panel oluşturuyorum. Amaç, başlangıç ​​değerleri üzerinde rastgele etki ile doğrusal olmayan bir regresyon gerçekleştirmektir.

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]

Bir bakarsak, uygun bir veri kümesine sahibiz:

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

Problemler

Nlme kullanma

Daha nlmebasit örnekler üzerinde çalışan benzer sözdizimi kullanarak (deSolve kullanmayan doğrusal olmayan fonksiyonlar), şunu denedim:

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

Ben şunu elde ederim:

CheckFunc'da hata (Func2, times, y, rho): func () (2) tarafından döndürülen türevlerin sayısı, başlangıç ​​koşulları vektörünün (2000) uzunluğuna eşit olmalıdır

İzleme:

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

Ben benziyor nlmebaşlangıç koşulu bir vektör geçmesi için çalışıyor solution_analy_ODE2ve hataya neden oluyor checkFuncdan lasoda.

Kullanmayı denedim 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

Doğrusal olmayan uyumun bireysel sinyallerde iyi çalıştığını görebiliriz. Şimdi, rasgele etkilerle veri kümesinin bir regresyonunu gerçekleştirmek istersem, sözdizimi şöyle olmalıdır:

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

Ama aynı hata mesajını alıyorum.

Daha sonra nlmODEBne Bolker'ın birkaç yıl önce sorduğum benzer bir soru hakkındaki yorumundan sonra kullanmayı denedim.

nlmODE kullanarak

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

Şu hatayı alıyorum:

Resnlmeode hatası (esp2omega, omega2, yeq, time, ID): 'yhat' nesnesi bulunamadı

Burada hatanın nereden geldiğini ve nasıl çözüleceğini anlamıyorum.

Sorular

• Sorunu yeniden oluşturabilir misin?
• nlmeVeya kullanarak bu sorunu çözmek için bir fikri olan var mı nlmODE?
• Değilse, başka bir paket kullanarak bir çözüm var mı? Ben testere nlmixr(https://cran.r-project.org/web/packages/nlmixr/index.html), ancak bilmiyorum, kurulum karmaşık ve kısa süre önce CRAN'dan kaldırıldı

---

Düzenlemeler

@tpetzoldt nlme, davranışlarda hata ayıklamanın güzel bir yolunu önerdi ve bu beni çok şaşırttı. İşte doğrusal olmayan bir işleve sahip çalışan bir örnek, burada bireyler arasında değişen rastgele bir parametre ile 5 kişilik bir grup oluşturuyorum:

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]

Fonksiyondaki kediler burada şunları verir:

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

Şimdi şununla yapıyorum 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))

Alırım:

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

Böylece nlme, zaman vektörünü 5 kez (bireysel sayısı) bağlar ve aynı sayıda tekrar eden parametrelerle bunu işleve aktarır. Elbette ki bu yöntemle uyumlu değil lsodave işlevim çalışıyor.

Answer 1

Öde modelinin yanlış bir bağımsız değişkenle çağrıldığı anlaşılıyor, bu nedenle 2 yerine 2000 durum değişkenli bir vektör elde ediyor. Sorunu görmek için aşağıdakileri deneyin:

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

Düzenleme : Analitik fonksiyonun işe yaradığını düşünüyorum, çünkü vektörize edildiğinden, ode modelini yineleyerek veya (daha iyisi) dahili olarak durum değişkenleri olarak vektörleri kullanarak ode fonksiyonunu vektörleştirmeyi deneyebilirsiniz. Gibi odebirçok 100k denklemlerle sistemlerini çözümünde hızlı 2000 uygulanabilir olmalıdır.

Sanırım hem durumların hem de parametrelerin nlmevektörler olarak aktarılıyor. Ode modelinin durum değişkeni daha sonra bir "uzun" vektördür, parametreler bir liste olarak uygulanabilir.

İşte bir örnek (şimdi parametrelerle birlikte liste halinde düzenlenmiştir):

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

Daha sonra denklem sayısını ayarlayın (veya hesaplayın), örneğin N <- 1resp. N <-1000aramalardan önce.

Model, sayısal konulara girmeden önce bu şekilde çalışır, ancak bu başka bir hikaye ...

Daha sonra , doğrusal olmayan optimizasyonda her zamanki gibi başka bir ode çözücü kullanmayı (örn. vode) Deneyebilir , değerleri ayarlayabilir atolve rtolazaltabilir, nmleoptimizasyon parametrelerinde ince ayar yapabilir , kutu kısıtlamalarını kullanabilirsiniz ... vb.

Answer 2

Bir çözüm hackleme nlmedavranışı buldum : Düzenlememde gösterildiği gibi, sorun nlme, fonksiyonun her zaman için bir değeri ilişkilendirdiğini varsayarak, Nind IndividualxNpoints vektörünün doğrusal olmayan fonksiyona aktarılmasından kaynaklanıyor. Ancak lsodabunu yapmayın, çünkü zaman boyunca bir denklemi bütünleştirir (yani, bir değer üretmek için belirli bir zaman aralığına kadar tüm zamana ihtiyacı vardır).

Çözümüm, nlmefonksiyonuma geçen parametreleri ayrıştırmak, hesaplamayı yapmak ve bir vektörü yeniden oluşturmaktan ibarettir :

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 zaman vektörünü tek zaman vektörleri tanımlayıcısına ayırın:

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

Ve sonra, her bir birey üzerinde sayısal entegrasyon döngüsünü gerçekleştiren işlev ve ortaya çıkan vektörleri birbirine bağlar:

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

İşleve nlmegeçenlere benzer bir vektörü geçtiğim bir test yapıyorum :

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

İşe yarıyor. @Tpetzoldt yöntemi ile çalışmaz, çünkü zaman vektörü 10'dan 0'a geçer, bu da entegrasyon sorunlarına neden olur. Burada gerçekten nlnmeişlerin yolunu kesmem gerekiyor . Şimdi:

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

tıkır tıkır çalışıyor

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