3

I have a dataset that looks something like this

d<–structure(list(groupid = c(2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 3L, 3L, 3L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 
3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 
1L, 2L, 2L, 2L, 1L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 3L, 3L, 
3L, 3L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L), participant_id = c(1L, 
1L, 1L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 
7L, 7L, 7L, 8L, 8L, 9L, 9L, 9L, 10L, 10L, 10L, 11L, 11L, 11L, 
12L, 12L, 13L, 13L, 13L, 14L, 14L, 14L, 15L, 15L, 15L, 16L, 16L, 
17L, 17L, 17L, 18L, 18L, 19L, 19L, 19L, 20L, 20L, 20L, 21L, 21L, 
21L, 22L, 22L, 22L, 23L, 23L, 24L, 24L, 24L, 25L, 25L, 26L, 26L, 
26L, 27L, 27L, 28L, 28L, 28L, 29L, 29L, 29L, 30L, 30L, 31L, 31L, 
31L, 32L, 32L, 32L, 33L, 33L, 34L, 34L, 34L, 35L, 35L, 35L, 36L, 
36L, 36L, 37L, 37L, 37L, 38L, 38L, 38L, 39L, 39L, 39L, 40L, 40L, 
40L, 41L, 41L, 41L, 42L, 42L, 42L, 43L, 43L, 43L, 44L, 44L, 45L, 
45L, 46L, 46L, 47L, 47L, 47L, 48L, 48L, 49L, 49L, 50L, 50L), 
    attrib1_A = c(0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 
    0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 
    0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 
    0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 
    0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 
    1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 
    0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 
    0, 0, 1, 1, 0), attrib1_B = c(1, 0, 0, 0, 0, 0, 0, 1, 1, 
    0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 
    0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 
    0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 
    0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 
    0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 
    0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 
    0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1)), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -134L), groups = structure(list(
    participant_id = 1:50, .rows = structure(list(1:3, 4:5, 6:8, 
        9:11, 12:14, 15:17, 18:20, 21:22, 23:25, 26:28, 29:31, 
        32:33, 34:36, 37:39, 40:42, 43:44, 45:47, 48:49, 50:52, 
        53:55, 56:58, 59:61, 62:63, 64:66, 67:68, 69:71, 72:73, 
        74:76, 77:79, 80:81, 82:84, 85:87, 88:89, 90:92, 93:95, 
        96:98, 99:101, 102:104, 105:107, 108:110, 111:113, 114:116, 
        117:119, 120:121, 122:123, 124:125, 126:128, 129:130, 
        131:132, 133:134), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), row.names = c(NA, -50L), class = c("tbl_df", 
"tbl", "data.frame"), .drop = TRUE))

# Groups:   participant_id [4]
   groupid participant_id attrib1_A attrib1_B
     <int>          <int>     <dbl>     <dbl>
 1       2              1         0         1
 2       2              1         1         0
 3       2              1         0         0
 4       1              2         0         0
 5       1              2         1         0
 6       2              3         1         0
 7       2              3         0         0
 8       2              3         0         1
 9       2              4         0         1
10       2              4         1         0

Where groupid indicates the cluster of the participant_id. attrib1 and attrib2 are regressors I want to use for the DGP of the variable I would like to create.

I would like to generate my binary outcome variable y that follows the Data Generating Process (DGP) specified below.

$y=a+factor(attrib1)+ factor(attrib2)$

Where a is the constant: the probability of y=1 when attrib1 and attrib2 correspond to reference categories. The betas are attrib1_A= 0.3 attrib1_B= -0.5

Finally, I would like that the variable y would be created flowing a specified intracluster correlation (e.g. 0.05).

Intracluster correlation is the between-cluster variability divided by the sum of the within-cluster and between-cluster variabilities. In our case clusters are groupid.

Does anyone know how I can generate a variable with a specified DGP and a specified intracluster correlation?

thanks a lot in advance for your help

Alex
  • 1,207
  • 9
  • 25

1 Answers1

1

Ok, so this answer is based on a simple formula of the ICC for binomial glmm, where the variance is fixed at pi^2 / 3

The first block of code in your question didn't work for me, and convergence is problematic for an ICC of .05, as you can imagine. But here is a setup for the problem. The answer uses the broom.mixed, arm, performance, and lme4 packages

library(broom.mixed)

b1 <- .3 # first coef
b2 <- -.5 # second coef
# select second stage variance and check ICC
second_stage_var <- 10 # this is the sd
icc <- second_stage_var^2 /(second_stage_var^2 + ( (pi ^ 2) / 3))
icc 

### generate data
id <- rep(letters, 50) #ids 
id_effect <- rnorm(length(letters), 0, second_stage_var)  #id effects
x1 <- rnorm(length(id), 10, 10) # first_covariate <
x2 <- rnorm(length(id), 30, 10) # second covariate <

df <- data.frame(id, id_effect, x1, x2) 

probs <- arm::invlogit(1 + b1*x1 + b2*x2 + id_effect)
df$y <- rbinom(length(id), 1, prob = probs)

mod1 <- lme4::glmer(data = df, formula = y ~ x1 + x2 + (1|id), family =binomial())
summary(mod1)

## check icc from the performance package
performance::icc(mod1)

## check icc by "hand"
est <- tidy(mod1)

total <- est$estimate[[4]]^2 +  (pi^2/3)
est$estimate[[4]]^2 / total
SushiChef
  • 588
  • 3
  • 12