3

I have written R codes for simulating data from Multilevel logistic regression model .

I focus on the following multilevel logistic model with one explanatory variable at level 1 (individual level) and one explanatory variable at level 2 (group level) :

$$\text{logit}(p_{ij})=\pi_{0j}+\pi_{1j}x_{ij}\ldots (1)$$ $$\pi_{0j}=\gamma_{00}+\gamma_{01}z_j+u_{0j}\ldots (2)$$ $$\pi_{1j}=\gamma_{10}+\gamma_{11}z_j+u_{1j}\ldots (3)$$

where , $u_{0j}\sim N(0,\sigma_0^2)$ , $u_{1j}\sim N(0,\sigma_1^2)$ , $\text{cov}(u_{0j},u_{1j})=\sigma_{01}$

In this paper in equation (2) , they assumed $\text{cov}(u_{0j},u_{1j})=\sigma_{01}$ , that is not independent . But also they mentioned in the methodology section that :

The group random components $u_{0j}$ and $u_{1j}$ are "independent" normal variables with mean zero and standard deviations $σ_0$ and $σ_1$.

So I assumed $\text{cov}(u_{0j},u_{1j})=0$ .

R code :

## Simulating data from multilevel logistic regression 

set.seed(1234)
x <- rnorm(1000) ### individual level variable
z <- rnorm(1000) ### group level variable

##fixed effect parameter
g_00 <- -1
g_01 <- 0.3
g_10 <- 0.3
g_11 <- 0.3

g <- matrix(c(g_00,g_01,g_10,g_11),ncol=1)

require(mvtnorm)

##need variance values as input 
s2_0 <- 0.36
s2_1 <- 1
s01 <- 0

##generate bi-variate normal rv for u0, u1

avg <- c(0,0) ##mean
sigma <- matrix(c(s2_0,s01,s01,s2_1),ncol=2)

u <- rmvnorm(1000,mean=avg,sigma=sigma,method="chol")

pi_0j <- g_00 +g_01*z + as.vector(u[,1])
pi_1j <- g_10 +g_11*z + as.vector(u[,2])
p <- exp(pi_0j+pi_1j*x)/(1+exp(pi_0j+pi_1j*x))

y <- rbinom(1000,1,p)

But i am not understanding where is to consider the group ? If i select number of groups to be $100$ $(j=1,2,\ldots, 100)$, then will I assign the groups randomly against each $y_{i,j}$ ?

  • Have i correctly simulated data from Multilevel Logistic Distribution ?
ABC
  • 1,367
  • 3
  • 13
  • 31
  • Since you said the paper assume $\text{cov}(u_{0j},u_{1j})=\sigma_{01}$ and you assume $\text{cov}(u_{0j},u_{1j})=0$ how can you produce the results of the papere? I think you will produce different data. – Deep North Jun 01 '15 at 06:52
  • @DeepNorth But I assumed $\text{cov}(u_{0j},u_{1j})=0$ because in the paper , they told that "The group random components u0j and u1j are `independent` normal variables with mean zero and standard deviations σ0 and σ1." And if two variable are independent , their covariance is also 0 . – ABC Jun 01 '15 at 13:35

1 Answers1

5

Have I correctly simulated data from Multilevel Logistic Distribution?

Well, maybe, but I don't think so. Think of it this way: You want to simulate the case where you have a sample of individuals, and each individual belongs to a group. The individuals have covariate values and the groups have covariate values.

The article is a little confusing at first. If you check out Maas and Joop (2004), you will see that Moineddin et al (2007) are following a very similar line of development. That might clarify some of their decisions.

Maas and Joop (2004) examine the situation in the multilevel regression situation. The model is quite similar, with the same exact structure for the $u_{0j}$ and $u_{1j}$. And there on page 131 we find the following statement:

To simplify the simulation model, without loss of generality, the covariance between the two u-terms is assumed equal to zero.

So, Moineddin et al (2007) are just following along. There is no reason you have to assume that $\sigma_{01} = 0$. Of course, to reproduce the results in the paper, you would.

They also are following the same line of development with respect to evaluating the effect of the intraclass correlation coefficient (ICC) on estimation. In the logistic model, the ICC is computed a little differently than in the linear regression situation. But, the bottom line is the same --- specifying an ICC equates to specifying a variance for $\sigma_0^2$.

Here is a revision of your code:

## Simulating data from multilevel logistic regression

require(mvtnorm)

set.seed(1234)

# Set up various parameters.

J   <- 30          ### number of groups
n_j <- rep(5, J)   ### number of individuals in jth group

# What is the sample size in this case?

N <- sum(n_j)

# Set the fixed effect parameters.

g_00 <- -1
g_01 <- 0.3
g_10 <- 0.3
g_11 <- 0.3

# Set the variance values for the group level random effect.

s2_0 <- 0.13    ### variance corresponding to specific ICC
s2_1 <- 1       ### variance standardized to 1
s01  <- 0       ### covariance assumed 0

# Simulate the covariate values for this sample size.

z <- rnorm(J)   ### group-level covariate values
x <- rnorm(N)   ### individual-level covariate values

At this point, you are all set to simulate. Notice that your original variance value was actually the standard deviation ($0.36^2 \approx 0.13$). That variance value was driven by the desire to have the ICC $=0.04$.

# Generate (u_0j, u_1j) from a bivariate normal.

mu    <- c(0, 0)
sigma <- matrix(c(s2_0, s01, s01, s2_1), ncol=2)

u <- rmvnorm(J, mean=mu, sigma=sigma, method="chol")

# Now form the linear predictor.  First, the group-level
# predictor.  There are J of these.

pi_0 <- g_00 + g_01 * z + as.vector(u[, 1])
pi_1 <- g_10 + g_11 * z + as.vector(u[, 2])

# Next, the individual level predictor.  There are
# N of these.

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)

Maas, Cora JM, and Joop J. Hox (2004) Robustness issues in multilevel regression analysis. Statistica Neerlandica 58.2: 127-137.

Moineddin, Rahim, Flora I. Matheson, and Richard H. Glazier (2007) A simulation study of sample size for multilevel logistic regression models. BMC Medical Research Methodology 7.1: 34.

jvbraun
  • 2,239
  • 13
  • 19