Relative bias of estimates of multilevel logistic regression models

[Cynderella]

The post is little bit lengthy , but it will give you thorough information and you can run the code easily . Trust me , I believe the answer is really very basic and not lengthy at all .

I want to calculate relative bias of estimates of multilevel logistic model with one explanatory variable at level 1 (individual level) and
one explanatory variable at level 2 (group level) to reproduce the result of this article

library(mvtnorm)
library(lme4)

set.seed(1234)

#Function for Simulating data from multilevel logistic distribution 
simu <- function(J,n,g00,g10,g01,g11,s2_0,s2_1,s01){
     n_j <- rep(n,J)     ## number of individuals in jth group
     N <- sum(n_j)       ## sample size

     #Simulate the covariate value for this sample size .
     z <- rnorm(J)
     x <- rnorm(N)

    #Generate (u_0j,u_1j) from a bivariate normal .
     mu <- c(0,0)
     sig <- matrix(c(s2_0,s01,s01,s2_1),ncol=2)
     u <- rmvnorm(J,mean=mu,sigma=sig,method="chol")

    #Now form the linear predictor .
     
    pi_0 <- g00 +g01*z + u[,1]
    pi_1 <- g10 + g11*z + u[,2]
    eta <- rep(pi_0,n_j)+rep(pi_1,n_j)*x

   #Transform back to the probability scale .
    p <- exp(eta)/(1+exp(eta))

    #Simulate a bernoulli from each individual distribution .
    y <- rbinom(N,1,p)

    sim_data_mat <- matrix(c(y,x,rep(z,n_j),rep(1:30,n_j)),ncol=4)
    sim_data <- data.frame(sim_data_mat)
    colnames(sim_data) <- c("Y","X","Z","cluster")
    return(sim_data)  
   }

# Function for Fixed effect estimates 

coeff <- function(J,n,g00,g10,g01,g11,s2_0,s2_1,s01){
                  sim_data <- simu(J,n,g00,g10,g01,g11,s2_0,s2_1,s01)
                  summary(glmer(Y~X*Z+(1|cluster),family=binomial,data=sim_data,nAGQ=10))$coef
             }

#estimated fixed effect parameter :
est_fix_par <- replicate(1000,coeff(30,5,-1,.3,.3,.3,.13,1,0))
avg_est_fix_par       <- apply(est_fix_par ,MARGIN=1:2,FUN=mean)
hat_fixed<- avg_est_fix_par [,1]

Then I tried to find the relative biases according to this table.

So I ran the following code :

#true value of parameter :
g00=-1;g10=.3;g01=.3;g11=.3
fixed <-c(g00,g10,g01,g11)

##Relative bias of Fixed Effect :
rel_bias_fixed <- ((hat_fixed-fixed)/fixed)*100

#(Intercept)           X           Z         X:Z 
#  -7.223689  -13.324452   -6.877180  -16.767887

which doesn't match with the result of the table .

I ran the commands for several combinations of the parameters but still the results seem to differ far away. Why ?

#Combination 2 : (J,n,g00,g10,g01,g11,s2_0,s2_1,s01)=(30,5,-1,.3,.3,.3,.67,1,0)

#estimated fixed effect parameter :
est_fix_par <- replicate(1000,coeff(30,5,-1,.3,.3,.3,.67,1,0))
avg_est_fix_par       <- apply(est_fix_par ,MARGIN=1:2,FUN=mean)
hat_fixed<- avg_est_fix_par [,1]

##Relative bias of Fixed Effect :
rel_bias_fixed <- ((hat_fixed-fixed)/fixed)*100

#(Intercept)           X           Z         X:Z 
# -9.275858  -18.063167   -6.449318  -14.032816

Why doesn't my result match with the table ? Is this for their random number from Normal distribution for individual level predictor (X) and group level predictor (Z) differ from me ? Or , where am I doing mistake ?

Thanks ! Regards .