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
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))
}
I am expecting a theta vector of length m. Unfortunately, theta comes as the size of one.
Any help is appreciated.
results of beta and 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.
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
}