如何从多个组中提取 lme() 函数的结果，然后在 R 中组合?

Question

首先，根据sl变量将以下数据随机分成两组，然后使用数据集下方所示的 for 循环为两组运行模型

mydata
y  x sl
1  5.297967  1  1
2  3.322833  2  1
3  4.969813  3  1
4  4.276666  4  1
5  5.972807  1  2
6  6.619440  2  2
7  8.045588  3  2
8  7.377759  4  2
9  6.907755  5  2
10 8.672486  6  2
11 8.283999  7  2
12 8.455318  8  2
13 7.414573  9  2
14 8.634087 10  2
15 7.356355  1  3
16 6.606247  2  3
17 6.396930  9  3
18 6.579251 10  3
19 5.521110  1  4
20 2.224221  2  4
21 6.742881  3  4
22 6.709304  4  4
23 6.875232  5  4
24 8.476371  6  4
25 7.360104  7  4

对两个组使用 lme(( 函数的 Runnign 模型，然后将beta系数存储为矩阵，theta[ 随机截距项 ] 作为向量

sl.no=unique(mydata$sl)
m=length(unique(mydata$sl))
ngrp=2
set.seed(125)
idx=sample(1:ngrp, size=m, replace = T)
beta=matrix(NA, nrow = ngrp, ncol=3, byrow=T) #null matrix to store coefficients from both groups 
theta=rep(0,m) #null vector to store intercepts from both groups
library(nlme)
for ( g in 1:ngrp){
rg=sl.no[idx==g]
mydata_rG=mydata[mydata$sl %in% rg,] #Data set belongs to group-g

lme_mod=lme(y~x+I(x^2),random = ~ 1|sl,
data = mydata_rG, method = "ML") #mixed effect model for each group

beta[g,]=c(unlist(lme_mod$coefficients[1])[[1]],
unlist(lme_mod$coefficients[1])[[2]],
unlist(lme_mod$coefficients[1])[[3]])
theta=c(unname(lme_mod$coefficients$random$sl))

}

我期待一个长度为 m 的theta向量，不幸的是，theta的大小只有一个。 任何帮助，不胜感激。

beta和theta的结果

beta
[,1]       [,2]        [,3]
[1,] 4.895805  0.7954474 -0.05602771
[2,] 6.423533 -1.7441753  0.32049662
theta
[1] 4.264366e-21 #it should be length of m.

Answer 1

它只是关于你如何存储theta值

sl.no=unique(mydata$sl)
m=length(unique(mydata$sl))
ngrp=2
set.seed(125)
idx=sample(1:ngrp, size=m, replace = T)
beta=matrix(NA, nrow = ngrp, ncol=3, byrow=T) 
theta=numeric() #Change here
library(nlme)
for ( g in 1:ngrp){
rg=sl.no[idx==g]
mydata_rG=mydata[mydata$sl %in% rg,] 
lme_mod=lme(y~x+I(x^2),random = ~ 1|sl,
data = mydata_rG, method = "ML") 

beta[g,]=c(unlist(lme_mod$coefficients[1])[[1]],
unlist(lme_mod$coefficients[1])[[2]],
unlist(lme_mod$coefficients[1])[[3]])
theta=c(theta, unname(lme_mod$coefficients$random$sl)) #Change here
}