การถดถอยที่ไม่ใช่เชิงเส้นพร้อมเอฟเฟกต์สุ่มและ lsoda

Jan 18 2021

ฉันกำลังเผชิญกับปัญหาที่ฉันไม่สามารถแก้ไขได้ ฉันต้องการใช้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และทำให้เกิดข้อผิดพลาดในจากcheckFunclasoda

ฉันลองใช้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และฟังก์ชั่นของฉันใช้งานได้

คำตอบ

3 tpetzoldt Jan 21 2021 at 01:21

ดูเหมือนว่าโมเดล 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 <- 1resp N <-1000ก่อนการโทร

โมเดลวิ่งผ่านทางนี้ก่อนที่จะทำงานในประเด็นตัวเลข แต่นั่นเป็นอีกเรื่องหนึ่ง

จากนั้นคุณอาจลองใช้ตัวแก้บทกวีอื่น (เช่นvode) ตั้งค่าatolและrtolลดค่าปรับแต่งnmleพารามิเตอร์การเพิ่มประสิทธิภาพใช้ข้อ จำกัด ของกล่อง ... และอื่น ๆ ตามปกติในการปรับให้เหมาะสมแบบไม่เชิงเส้น

1 denis Jan 29 2021 at 20:07

ฉันพบวิธีแก้ไข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