Can we generate simulated data using result of Logistic regression? 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.
 A: 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")
summary(model)

This gives (output truncated):
Coefficients:
            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

summary(glm(sim~xx_new,family="binomial"))

The result is quite close to the coefficients we found above and used for the simulation:
Coefficients:
            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 ***

