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
Then I tried to find the relative biases according to this table.
So I ran the following code :
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 ?
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 .
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
Code:
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]
So I ran the following code :
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
I ran the commands for several combinations of the parameters but still the results seem to differ far away. Why ?
Code:
#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
Thanks ! Regards .
Last edited: