# How to simulate interaction between categorical and continuous variables over a binary outcome?

For a power analysis, I need to simulate data to run a glmer model with an interaction between categorical (three levels) and continuous (which is going to be the moderator; high-score participants vs. low-score participants) over a binary outcome (0,1). I don't know how to define different effect size or mean differences for each condition. I expect that high-score participants have higher outcome proportion in level a, moderate proportion in level b, and low proportion in level c, whereas low-score participants have lower proportion in level a and b, but moderate in level c. Or such detail are not relevant for a power analysis?

• You won't be able to simulate a power analysis without knowing the proposed effect size. You need to start with the effect you want to be able to capture, then you work out how to simulate that. Commented Jul 22 at 11:21
• While it won't quite answer your question, it might help you to read this thread: Simulation of logistic regression power analysis - designed experiments. Commented Jul 22 at 11:23
• @gung-ReinstateMonica: your comments look like an answer to me. Want to post as such? Commented Jul 22 at 11:29
• Thanks, @StephanKolassa. I don't have time at the moment, but I can circle back later today. Commented Jul 22 at 11:32
• @Reinstate Monica Thank you. I know the proposed effect sizes, but I don't know how to get around R to make this. Commented Jul 22 at 12:07

Considering this has the tag, the easiest way is by using the same function R uses internally when you enter a model formula: model.matrix.

This function provides us with a design matrix $$\mathbf{X}$$, that we can plug into the formula for logistic regression:

\begin{aligned} \mathbb{E}[y] &\sim \mathsf{binomial}(n, p) \\\\ \log\left( \frac{p}{1-p} \right) &= \mathbf{\eta} \\\\ \mathbf{\eta} &= \mathbf{X\beta} \end{aligned}

I will start with a non-sensically low sample size, so it's easy to see what each step does.

## Generating the design matrix $$\mathbf{X}$$

You did not mention what values the score can take on. For this example I will use uniformly random scores between $$0$$ and $$1$$.

set.seed(2024)
n <- 3 # Sample size per group
k <- 3 # Number of groups

group <- factor(rep(LETTERS[1:k], each = n))
score <- runif(n * k, 0, 1) # uniformly random scores between 0 and 1

X <- model.matrix(~ group * score) # The design matrix X


You can see what's what by looking at X:

#   (Intercept) groupB groupC     score groupB:score groupC:score
# 1           1      0      0 0.8369425    0.0000000    0.0000000
# 2           1      0      0 0.3208675    0.0000000    0.0000000
# 3           1      0      0 0.6803633    0.0000000    0.0000000
# 4           1      1      0 0.6981731    0.6981731    0.0000000
# 5           1      1      0 0.4570092    0.4570092    0.0000000
# 6           1      1      0 0.7014203    0.7014203    0.0000000
# 7           1      0      1 0.4157110    0.0000000    0.4157110
# 8           1      0      1 0.3032091    0.0000000    0.3032091
# 9           1      0      1 0.8765881    0.0000000    0.8765881



The columns of this matrix are the parameters of our model, the same way they would be printed in a standard regression table. The interaction term makes up the last two columns.
(A difference in the slope of score if you belong to group B or C.)

## Generating the effect sizes $$\mathbf{\beta}$$

Now we add coefficients for each of the parameters in the model. Note that these are on the scale of the linear predictor $$\mathbf{\eta}$$.

beta <- numeric(ncol(X))
beta[1] <- -1 # (Intercept):  Group A at a score of 0
beta[2] <-  0 # groupB:       Difference between B and A at a score of 0
beta[3] <-  2 # groupC:       Difference between C and A at a score of 0
beta[4] <-  2 # score:        The effect of score for group A
beta[5] <- -1 # groupB:score: The difference in effect of score between B and A
beta[6] <- -3 # groupC:score: The difference in effect of score between C and A


## Simulating the data

Now we have everything we need to generate the linear predictor $$\mathbf{\eta}$$ and plug it into a binomial random number generator. To go from $$\mathbf{\eta}$$ to $$p$$, you need to use the inverse of the link function:

$$p = \frac{1}{1 + e^{-\mathbf{\eta}}}$$

You could write your own function, but this is also readily available by using binomial()$linkinv: eta <- X %*% beta # The linear predictor y <- rbinom(nrow(X), size = 1, prob = binomial()$linkinv(eta)) # The outcome

# Combine all variables in a data frame and clean up
DF <- data.frame(group, score, y)
rm(n, k, X, beta, eta, y)


## The result

Now I will repeat all steps with a more sensible sample size (e.g., $$n = 500$$ per group) and confirm that it does what we're hoping for:

GLM <- glm(y ~ group * score, family = "binomial", data = DF)
summary(GLM)$coefficients # Estimate Std. Error z value Pr(>|z|) # (Intercept) -1.1966702 0.1962809 -6.096724 1.082644e-09 # groupB 0.1637966 0.2820374 0.580762 5.614009e-01 # groupC 2.3367624 0.2772178 8.429337 3.476284e-17 # score 2.3120185 0.3331853 6.939137 3.945025e-12 # groupB:score -1.4708693 0.4784921 -3.073968 2.112323e-03 # groupC:score -3.6971650 0.4734487 -7.809009 5.763923e-15  Fairly close to what we simulated. Let's see what that looks like in terms of probability of success: require("sjPlot") plot_model(GLM, type = "pred", terms = c("score", "group"))  This assumes a linear effect of time on the scale of the linear predictor, resulting in the characteristic S-shape of logistic regression. For more complex relationships, you could play around with splines, like in the simulation used here. ## The whole script set.seed(2024) n <- 100 k <- 3 group <- factor(rep(LETTERS[1:k], each = n)) score <- runif(n * k, 0, 1) X <- model.matrix(~ group * score) beta <- numeric(ncol(X)) beta[1] <- -1 # (Intercept): Group A at a score of 0 beta[2] <- 0 # groupB: Difference between B and A at a score of 0 beta[3] <- 2 # groupC: Difference between C and A at a score of 0 beta[4] <- 2 # score: The effect of score for group A beta[5] <- -1 # groupB:score: The difference in effect of score between B and A beta[6] <- -3 # groupC:score: The difference in effect of score between C and A eta <- X %*% beta # The linear predictor y <- rbinom(nrow(X), size = 1, prob = binomial()$linkinv(eta))
DF <- data.frame(group, score, y)
rm(n, k, X, beta, eta, y)


The answer by @Frand Rodenburg is completely correct. There is a more flexible way to do this though, using the simDAG R package (full disclosure: I am the developer of that package). First, install the developmental version from github:

devtools::install_github("RobinDenz1/simDAG")


(It is also available on CRAN, but the formula interface is not yet).

You could then use the following code to get the same kind of data that was shown in the previous answer:

library(simDAG)

dag <- empty_dag() +
node("group", type="rcategorical", probs=c(0.3333, 0.3333, 0.3333),
coerce2factor=TRUE, labels=c("A", "B", "C")) +
node("score", type="runif", min=0, max=1) +
node("outcome", type="binomial", formula= ~ -1 + groupB*0 + groupB*2 +
score*2 + groupB:score*-1 + groupC:score*-3,
coerce2numeric=TRUE)

data <- sim_from_dag(dag, n_sim=1000)


You can check that this works by setting n_sim to a very high value to generate a huge dataset and next fitting the corresponding logistic regression model:

mod <- glm(outcome ~ group + score + group*score, data=data, family="binomial")
summary(mod)


By adding some more node() calls to the dag object, it is possible to generate arbitrarily complex data. This is explained in more detail in the documentation of the package and in the associated vignettes.