1
$\begingroup$

I want to simulate data from multilevel logistic regression .

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}$

I prefer R and started to write codes for the simulation as :

 set.seed(36)

 x <- rnorm(1000)  ### individual level variable , x_ij

 z <- rnorm(1000)  ### group level variable , z_j

If I have initial value for $\gamma_{00}=-1 , \gamma_{01}=0.3,\gamma_{10}=0.3,\gamma_{11}=0.3$ , how can I generate $\pi_{0j},\pi_{1j}$ since there is $u_{0j},u_{1j}$ in equation (2) and (3) ?

$\endgroup$
2
$\begingroup$

In your case, basically $\boldsymbol u = (u_0, u_1)$ is distributed as bi-variate normal. You can use R package mvtnorm to generate $u_{0j}, u_{1j}$. Then you can get $\pi_{0j}, \pi_{1j}$ easily. (You would wish to do all those in terms of random vector manipulation.)

The following code generates $\boldsymbol u = (u_0, u_1)$ for you.

require(mvtnorm)
set.seed(1234)
## need variance values as input
s2_0 <- 2
s2_1 <- 3
s01 <- 0.5

## generate bi-variate normal rv for u0, u1
mean <- c(0, 0)
sigma <- matrix(c(s2_0, s01, s01, s2_1), ncol = 2)

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

Here I assume you know the values of the hyperparameters. And you get something like following.

> dim(u)
[1] 1000    2
> head(u, 10)
             [,1]       [,2]
 [1,]  0.03563353 -2.0906992
 [2,] -2.93442022  1.8783072
 [3,]  0.82448029  0.7432658
 [4,] -0.92269067 -0.9954788
 [5,] -1.39514484 -0.9776595
 [6,] -1.51995711 -0.8265220
 [7,] -0.13484601 -1.3445112
 [8,]  0.12429872  1.6618925
 [9,] -1.40900867 -0.8850944
[10,]  3.10290171 -1.4500239
$\endgroup$
  • $\begingroup$ Could you please tell me what is hyperparameter ? $\endgroup$ – ABC Jun 1 '15 at 4:16
  • $\begingroup$ @ABC, I mean $\sigma_0^2, \sigma_1^2$ and $\sigma_{01}$ in your case. Basically a hypterparameter is the one for your prior distribution, not the one you wish to estimate and make inference about. That is mainly a term in Bayesian statistics. I might misuse it here. Apologize for misleading. $\endgroup$ – Aaron Zeng Jun 1 '15 at 4:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.