0

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.
Uddin
  • 754
  • 5
  • 18
  • `theta[g]=...`? By the way, doesn't *sl* have 4 groups, not 2? – Parfait Jan 13 '20 at 00:13
  • there are two groups, according to the split, group 1 has three `sl`, and group 2 has one. So for group one, I will have three `theta`s for group-1and one` theta` for the group-2 – Uddin Jan 13 '20 at 00:19
  • please see `idx`. it takes value either `1 ` (for group-1) or `2` (for group-2). The `rg` is created to identify which of the `sl` should go for group-g`(1 or 2)` – Uddin Jan 13 '20 at 01:09
  • Got it. Very subtle. Usually for groups, it is expected to be reflected in data. So does my earlier question to index `theta` with `g` not resolve your issue? – Parfait Jan 13 '20 at 01:27
  • No, `theta[g]` is not working :( – Uddin Jan 13 '20 at 01:28
  • Here is the warning message `In theta[g] <- c(unname(lme_mod$coefficients$random$sl)) : number of items to replace is not a multiple of replacement length` – Uddin Jan 13 '20 at 01:40
  • Issue may be right-hand side as it may return `NULL` or more than one item. Wrap with `tryCatch` and check NAs in result *theta* vector: `theta[g] <- tryCatch(c(unname(lme_mod$coefficients$random$sl)), warning = function(w) NA, error = function(e) NA)` – Parfait Jan 13 '20 at 01:50

1 Answers1

1

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

}
Ronak Shah
  • 377,200
  • 20
  • 156
  • 213