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

【问题标题】：How to extract result from lme() function from multiple groups and then combine in R?如何从多个组中提取 lme() 函数的结果，然后在 R 中组合？
【发布时间】：2020-01-13 00:07:34
【问题描述】：

首先，将以下数据根据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() 函数的运行模型，然后将 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.

【问题讨论】：

theta[g]=...?对了，sl不是有4组，不是2组吗？
有两组，根据拆分，第一组有三个sl，第二组有一个。所以对于第一组，我将有三个thetas 用于 group-1 和一个`theta` 用于 group-2
请参阅idx。它取值 1 （用于 group-1）或 2（用于 group-2）。创建rg 是为了识别sl 中的哪一个应该用于group-g(1 or 2)
知道了。很微妙。通常对于组，它预计会反映在数据中。那么我之前提出的用g 索引theta 的问题不能解决您的问题吗？
不，theta[g] 不工作 :(

标签： r extract resultset

【解决方案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

}

【讨论】：