如何从多个组的lme()函数提取结果，然后在R中合并? [英] How to extract result from lme() function from multiple groups and then combine in R?

问题描述

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

First, the following data are split randomly into two groups according to the sl variable and then run the model for both groups using the for loop shown below the data set

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 [随机拦截项]存储为向量

Runnign model using lme() function for both group and then store the beta coefficients as matrix and theta[ random intercept term ] as vector

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的大小为1. 任何帮助表示赞赏.

I am expecting a theta vector of length m. Unfortunately, theta comes as the size of one. Any help is appreciated.

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值的方式有关

It's only about how you store theta values

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

}

这篇关于如何从多个组的lme()函数提取结果，然后在R中合并?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！