Calculating Non-Coverage Probability .

#1
Hi , I have a two-level model with one explanatory variable (X) on the respondent level, and one group-level explanatory variable (Z) . I am following this paper http://joophox.net/publist/methodology05.pdf

I am simulating data from the model . Three conditions are varied in the simulation:
(1) number of groups (NG: three conditions, NG = 30, 50, 100),
(2) group size (GS: three conditions, GS = 5, 30, 50),
(3) intraclass correlation (ICC: three conditions, ICC = 0.1, 0.2, 0.3).

So I wrote a function to simulate data from the model:

Code:
# J               = number of groups
# n_j             = group size

simu <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
     N <- sum(rep(n_j,J))  #sample size
     x <- rnorm(N)         #level 1 (individual) level variable
     z <- rnorm(J)         #level 2 (group) level variable
     
     #generate (u_0j,u_1j) from bivariate normal
     mu <- c(0,0)
     sig <- matrix(c(sig2_0,sig01,sig01,sig2_1),ncol=2)
     u   <- rmvnorm(J,mean=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,0.5)
      sim_data <- data.frame(Y=y,X=x,Z=rep(z,each=n_j),group=rep(1:J,each=n_j))
 }  # end of function simu
Next , I wrote a function to estimate the parameters and To assess the accuracy of the standard errors, for each parameter in each simulated data set , the 95% confidence interval was established using the asymptotic standard normal distribution . For each parameter a noncoverage indicator variable was set up that is equal to zero if the true value is in the confidence interval, and equal to one if the true value is outside the confidence interval.

Code:
CI  <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){
         dat      <- simu(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1)
         lmer_obj <- lmer(Y~X+Z+X:Z+(X|group),data=dat)
         res      <- summary(lmer_obj)
         estim    <-  coefficients(res)
         se       <- coefficients(res)
         ci       <-  data.frame(intercept=estim[1] + qnorm(c(.025,.975))*se[1],
                         X=estim[2] + qnorm(c(.025,.975))*se[2],Z=estim[3] + qnorm(c(.025,.975))*se[3],
                         XZ=estim[4] + qnorm(c(.025,.975))*se[4])
         indicator <- ifelse((ci[1,]<c(g00,g10,g01,g11) & ci[2,]>c(g00,g10,g01,g11)),0,1)
   }
For your kind consideration , to make it clear for passing the input in the argument of function CI :

(1) sig01 = 0 , since in the paper , it is mentioned that "To simplify the simulation model, the covariance between the two u-terms is assumed equal to zero."

(2) sig2_0 = sig2_1 , in the paper , it is mentioned that "Busing (1993) shows that the effects for the intercept variance \(\sigma_{00}\) and the slope variance \(\sigma_{11}\) are similar; hence, we chose to set the value of \(\sigma_{11}\) equal to \(\sigma_{00}\) .

Now , I checked for a combination :

Code:
cond4 <- replicate(1000,CI(30,30,1,.3,.3,.3,(1/18),0,(1/18)))
nonc4 <- apply(cond4,1:2,mean)
which produces result

Code:
 intercept X Z XZ
1         0 0 0  0
It is HORRIBLY a wrong answer . The combination that I show you lies in row 4 and the result need to match with that row and column 7,8,9,10 of "Table 4" of the paper .

Where am I doing mistake ?

Any HELP is APPRECIATED .

Many thanks! Regards .