I have an outcome y (binary) and two variables x1 and x2 (continuous and independent). I have fitted logistic regression y~x1 and y~x2 in two different datasets, but I don't have a dataset that records both x1, x2 and y. I wish to generate a simulated dataset that follows the two fitted logistic regression (i.e. when I fit y~x1 and y~x2 separately in the simulated dataset, the coefficients should be close to the coefficients I obtained in the two real datasets). I have no clue how to achieve this using R and I wish to learn from you. Any help would be appreciated.


1 Answer 1


Recall that logistic regression with a single predictor $x$ models the probability for $y$ to be TRUE (or 1) as

$$ P(y=1) = \frac{1}{1+\exp(-\beta_0+\beta_1x)}. $$

We can now approach your question in two ways. Either you have your estimated coefficients $\hat{\beta}_0$ and $\hat{\beta}_1$. Or you have a fitted model object from a previous call to glm(). Let's do both ways, by first doing the first, then fitting a model to this first simulation, finally simulating from this fitted model.

n_sim <- 100           # for the initial dataset
set.seed(2)            # for replicability
xx <- runif(n_sim)     # predictor values
coefficients <- c(1,1) # my assumption
prob <- 1/(1+exp(-(coefficients[1]+coefficients[2]*xx)))

This last line corresponds to the formula for $P(y=1)$ above. We now draw $y$ according to these probabilities:

yy <- runif(n_sim)<prob

This is the core of the simulation.

We can now fit a model to these data:

model <- glm(yy~xx,family="binomial")

This gives (output truncated):

            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   0.9410     0.4706   2.000   0.0455 *
xx            1.2604     0.9246   1.363   0.1728

This is pretty close to the (1,1) coefficients we specified above. If we re-run this multiple times with different seeds, we find that the fitted coefficients vary quite wildly.

Finally, we can simulate according to the model object that contains the fitted model. This is just a question of calling predict.glm() to obtain the new predicted_probabilities and then simulating, as above:

n_sim_new <- 1000
xx_new <- runif(n_sim_new)
predicted_probabilities <- predict.glm(model,newdata=data.frame(xx=xx_new),type="response")
sim <- runif(n_sim_new)<predicted_probabilities


The result is quite close to the coefficients we found above and used for the simulation:

            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.7383     0.1458   5.064 4.11e-07 ***
xx_new        1.5356     0.2850   5.387 7.17e-08 ***

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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