การถดถอยที่ไม่ใช่เชิงเส้นพร้อมเอฟเฟกต์สุ่มและ 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
กับไวยากรณ์ที่คล้ายกันซึ่งทำงานกับตัวอย่างที่ง่ายกว่า (ฟังก์ชันที่ไม่ใช่เชิงเส้นที่ไม่ใช้ 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, ครั้ง, 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
เราจะเห็นว่า te non linear fit ทำงานได้ดีกับแต่ละสัญญาณ ตอนนี้ถ้าฉันต้องการทำการถดถอยของชุดข้อมูลด้วยเอฟเฟกต์แบบสุ่มไวยากรณ์ควรเป็น:
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, เวลา, 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
และฟังก์ชั่นของฉันใช้งานได้
คำตอบ
ดูเหมือนว่าโมเดล 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
รวดเร็วในการแก้ระบบที่มีสมการ 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
resp N <-1000
ก่อนการโทร
โมเดลวิ่งผ่านทางนี้ก่อนที่จะทำงานในประเด็นตัวเลข แต่นั่นเป็นอีกเรื่องหนึ่ง
จากนั้นคุณอาจลองใช้ตัวแก้บทกวีอื่น (เช่นvode
) ตั้งค่าatol
และrtol
ลดค่าปรับแต่งnmle
พารามิเตอร์การเพิ่มประสิทธิภาพใช้ข้อ จำกัด ของกล่อง ... และอื่น ๆ ตามปกติในการปรับให้เหมาะสมแบบไม่เชิงเส้น
ฉันพบวิธีแก้ไขnlme
พฤติกรรมการแฮ็ก: ดังที่แสดงในการแก้ไขของฉันปัญหามาจากการที่nlme
ส่งเวกเตอร์ของ NindividualxNpoints ไปยังฟังก์ชัน nonlinear โดยสมมติว่าฟังก์ชันเชื่อมโยงกันในแต่ละครั้งที่ชี้ค่า แต่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