4
$\begingroup$

This pdf illustrates nicely how is to test the random effect of multilevel model . But I am simulating data from a two-level model and estimating the parameters of the model for various combination of the parameters. For each condition , I generated 1000 simulated data sets. I have used R for both simulation and estimation. The codes are following :

simfun <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
     N <- sum(rep(n_j,J))  

     x <- rnorm(N)         
     z <- rnorm(J)         

     mu <- c(0,0)
     sig <- matrix(c(sig2_0,sig01,sig01,sig2_1),ncol=2)
     u   <- MASS::mvrnorm(J,mu=mu,Sigma=sig)


     b_0j <- g00 + g01*z + u[,1]
     b_1j <- g10 + g11*z + u[,2]

      y <- rep(b_0j,each=n_j)+rep(b_1j,each=n_j)*x + rnorm(N,0,sqrt(0.5))
     sim_data <- data.frame(Y=y,X=x,Z=rep(z,each=n_j),group=rep(1:J,each=n_j))

  } 

fit <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
    dat <- simfun(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1)
    model <- lmer(Y~X+Z+X:Z+(X||group),data=dat)
 }

c4 <- replicate(1000,fit(30,30,1,.3,.3,.3,(1/18),0,(1/18)))

Now I want to test the both fixed effect and random effects . But not understanding how is to test all of them simultaneously .

To test both fixed effect and random effects , I modified the last function fit as :

fit <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
    dat <- simfun(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1)
    full <- lmer(Y~X+Z+X:Z+(X||group),data=dat,control=lmerControl(optCtrl=list(maxfun=20000)))

  #Testing significance of random intercept
  null.U0 <- update(full, .~.-(1 | group))
  dev.U0 <- as.numeric(2*(logLik(full)-logLik(null.U0)))
  p.U0 <- 0.5*(1-pchisq(dev.U0,1))

 #Testing significance of random slope
  null.U1 <- update(full, .~.-(0 + X | group))
  dev.U1 <- as.numeric(2*(logLik(full)-logLik(null.U1)))
  p.U1 <- 0.5*(1-pchisq(dev.U1,1))

 #Testing significance of intercept of fixed part
 null.int <- update(full, .~.-1)
 dev.int <- as.numeric(2*(logLik(full)-logLik(null.int)))
  p.int <- (1-pchisq(dev.U1,1))

 #Testing significance of X of fixed part
 null.x <- update(full, .~.-X)
 dev.x <- as.numeric(2*(logLik(full)-logLik(null.x)))
  p.x <- (1-pchisq(dev.x,1))

#Testing significance of Z of fixed part
 null.z <- update(full, .~.-X)
 dev.z <- as.numeric(2*(logLik(full)-logLik(null.z)))
  p.z <- (1-pchisq(dev.z,1))

#Testing significance of interaction part
 null.xz <- update(full, .~.-X:Z)
 dev.xz <- as.numeric(2*(logLik(full)-logLik(null.xz)))
  p.xz <- (1-pchisq(dev.xz,1))

 pvals <- data.frame(p.U0=p.U0,p.U1=p.U1,p.int=p.int,p.x=p.x,p.z=p.z,p.xz=p.xz)
}

c1 <- replicate(7,fit(30,5,1,.3,.3,.3,(1/18),0,(1/18)))
colMeans(apply(c1,2,unlist))

But I think I am introducing the multiplicity error. How can I test the random and fixed effect part by Bonferroni Correction ?

Improve this question
$\endgroup$
2
  • $\begingroup$ I cannot re-run your code."Error in model.frame.default(data = dat, drop.unused.levels = TRUE, formula = Y ~ : variable lengths differ (found for 'X || group')" $\endgroup$
    – Deep North
    Commented Aug 9, 2015 at 10:16
  • $\begingroup$ I think your in code "model <- lmer(Y~X+Z+X:Z+(X||group),data=dat)" the variables need to change to lower case. I am wondering how your program still can run. Before you run your codes, i think you may add rm(list=ls()) $\endgroup$
    – Deep North
    Commented Aug 11, 2015 at 12:02

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.