# Simulate data based on output from a generalized linear model in R

I have a generalized linear model for a data set and I want to use the model to simulate data that matches the data used to generate the model. I am not wanting to resample the original data, but rather use the model output to generate data that matches original data.

For example, suppose I have a binomial data set in which survivorship is a function of size.

size <- c(3, 5, 10, 5, 17, 7, 17, 5, 18, 20, 9, 4, 17, 14, 8, 19, 11,
3, 12, 14, 4, 11, 10, 22, 10, 12, 10, 7, 3, 10, 9, 5, 4, 23,
8, 7, 7, 11, 19, 12, 3, 9, 6, 6, 8, 14, 18, 10, 9, 11)

surv <- c(1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1,
0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0,
1, 0, 1, 1, 0, 1, 0, 1, 0)

dat=data.frame(surv,size)
out.glm = glm(surv ~ size, family=binomial, data = dat)


I can graph the model predictions and standard errors, but how do I use the output of the glm to simulate data that should theoretically match the original dataset? Specifically, should I use the standard error of the mean as the input to rbinom() across the range of probabilities and then resample from the output?

• Check ?simulate -- it's a generic functions from simulating data from given R object, it is implemented for many R methods among other S3 methods available for them (i.e. predict, summary). – Tim Jan 26 '17 at 14:37

A generalized linear model specification consists of a distribution family and a link function (depending on the family, it may also include a free dispersion parameter).

The model's parameter estimates for the mean allow us to estimate the conditional mean, $\hat{\mu}(\mathbf{x})=E[Y|\mathbf{x}]$.

Assuming you want the estimated parameters as the population parameters (to simulate just from the fitted model as if it were the population model), the general steps would be:

• At each predictor combination (at each $\mathbf{x}_i$), compute $\hat{\mu}(\mathbf{x}_i)$. There will typically be a function to do this computation of fitted value automatically.

• obtain an estimate (or identify) the dispersion parameter if needed (e.g. in a gamma GLM you'll require an estimate of dispersion, whether from the GLM output or from other means). You won't need one for an ordinary binomial GLM (i.e. skip this step in that case).

• If necessary reparameterize from those parameter estimates to a convenient parameterization for random generation (supplied generators won't usually be in mean-dispersion form) - these will be different for every observation. (again, you won't need this for the binomial)

• use random number generating routines to simulate from that conditional distribution for each observation (predictor combination $\mathbf{x}_i$).

[In terms of your original example problem the entire set of steps above is doable in a single line of R code - rbinom(length(surv),1,fitted(out.glm,method="response"))]

Such simulations are random and don't match the original data set (the responses will be different in general), but should have similar characteristics

If you want to allow for the fact that the parameter estimates themselves are uncertain it's more complicated, but can at least be done approximately.

• Thanks for the detailed explanation. What I am asking is how to actually simulate the data using R. I think that I need to create a sample set using rbinom(), and then resample, but just not sure how to code it. – user44796 Jan 26 '17 at 14:32
• If you're really just asking for R code your post is most likely off topic here. Nonetheless I've included a suggestion which generates a sample from the fit. If you need more detail than that, you should probably ask it elsewhere. – Glen_b Jan 26 '17 at 14:34
• You hinted at my question, and that is whether it is better capture the uncertainty in the parameter estimates, or simply use the parameter estimates directly. I will test both ways and see (1) if I can recover the parameters and (2) if it makes a difference. – user44796 Jan 26 '17 at 15:14
• Including uncertainty in the mean (which of itself won't be an input to rbinom) will not make the simulated data more like the original data.. It's not clear what you mean by "better" above ... better at what, exactly? – Glen_b Jan 26 '17 at 15:20
• Sorry @Glen_b, I am not being clear. I was thinking about the binomial distribution like a normal distribution. As you clearly stated, since the binomial distribution doesn't have a dispersion parameter, using the parameter estimates alone is fine. If I was simulating a normal distribution, I would need to include the sd from the original data. – user44796 Jan 26 '17 at 15:37