यादृच्छिक प्रभाव और lsoda के साथ गैर रेखीय प्रतिगमन
मैं एक समस्या का सामना कर रहा हूं जिसे मैं हल करने का प्रबंधन नहीं करता हूं। मैं निर्धारित गुणांक (एक नम थरथरानवाला) के साथ एक दूसरे क्रम के अंतर समीकरण के समाधान के रूप में उपयोग करके यादृच्छिक प्रभाव के साथ एक गैर रेखीय प्रतिगमन का उपयोग करना 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सरल उदाहरणों पर काम करने वाले समान सिंटैक्स का उपयोग करना (गैर रेखीय कार्यों का उपयोग नहीं कर रहा है), मैंने कोशिश की:
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, और में checkFuncसे एक त्रुटि का कारण बनता है 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 पर Bne Bolker की टिप्पणी का उपयोग करके कोशिश की
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) में त्रुटि: ऑब्जेक्ट 'y' 'पाया गया
यहाँ मुझे समझ नहीं आ रहा है कि त्रुटि कहाँ से आती है, न ही इसे कैसे हल किया जाए।
प्रशन
- क्या आप समस्या को पुन: उत्पन्न कर सकते हैं?
- किसी को भी इस समस्या को हल करने के लिए एक विचार है,
nlmeया तो का उपयोग करnlmODE? - यदि नहीं, तो क्या अन्य पैकेज का उपयोग करके कोई समाधान है? मैंने देखा
nlmixr(https://cran.r-project.org/web/packages/nlmixr/index.html), लेकिन मुझे यह पता नहीं है, वृत्ति जटिल है और इसे हाल ही में सीआरएएन से हटा दिया गया था
संपादित करता
@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
...
तो nlme5 बार (व्यक्तिगत की संख्या) समय वेक्टर को बांधता है और इसे फ़ंक्शन में पास करता है, मापदंडों के साथ समय की समान संख्या दोहराई जाती है। जो निश्चित रूप से lsodaमेरे काम करने के तरीके के अनुकूल नहीं है ।
जवाब
ऐसा लगता है कि ऑड मॉडल को एक गलत तर्क के साथ कहा जाता है, ताकि इसे 2 के बजाय 2000 राज्य चर के साथ एक वेक्टर मिल जाए। समस्या को देखने के लिए निम्नलिखित प्रयास करें:
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कई 100k समीकरणों के साथ सिस्टम को हल करने में तेज है, 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कॉल से पहले।
मॉडल संख्यात्मक मुद्दों में चलने से पहले इस तरह से चलता है, लेकिन यह एक और कहानी है ...
फिर आप एक और ode solver (जैसे vode) का उपयोग करने की कोशिश कर सकते हैं , सेट कर सकते हैं atolऔर rtolनिचले मानों के लिए, nmle's ऑप्टिमाइज़ेशन पैरामीटर, बॉक्स की कमी का उपयोग कर सकते हैं ... और इसी तरह, नॉनलाइनियर ऑप्टिमाइज़ेशन में।
मुझे एक समाधान हैकिंग nlmeव्यवहार मिला : जैसा कि मेरे संपादन में दिखाया गया है, समस्या इस तथ्य से आती है कि nlmeनिंडलियर फ़ंक्शन के लिए निंडलिप्लेक्सनपॉइंट्स के एक वेक्टर को पास करता है, यह मानते हुए कि फ़ंक्शन सहयोगी प्रत्येक समय के लिए एक मान इंगित करता है। लेकिन 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