# How can I simulate new observations to a fitted model object while specifying the covariates in R?

I have a model that is fitted with glm and I want to be able to calculate the response of this model when I specify the independent variables, while simulating the error terms. For example, if my formula is y~x1*x2, then this corresponds to:

$y_{i} = \mu + \beta_{1} x_{1i} + \beta_{2} x_{2i} + \beta_{3} x_{1i} x_{2i} + \epsilon_{i}$

where $\epsilon_{i}$ is an error term.

I want to be able to simulate

$y^{\star}_{i} = \hat{\mu} + \hat{\beta}_{1} x^{\star}_{1i} + \hat{\beta}_{2} x^{\star}_{2i}+ \hat{\beta}_{3} x^{\star}_{1i}x^{\star}_{2i} + \epsilon_{i}$,

with $\hat{\mu}$, $\hat{\beta}_{1}$, $\hat{\beta}_{2}$, and $\hat{\beta}_{3}$ being the values estimated by glm and $x^{\star}_{1 i}$ and $x^{\star}_{2 i}$ being the specified covariates / independent variables.

This is similar to what predict() does, however, it sets $\epsilon_{i}=0$.

It looks like what I want to do is similar to the functionality provided by simulate() and predict(). Although it doesn't look like simulate() allows you to specify the covariates.

What is the best way to do this?

-
It's not clear to me what you mean by "specify the parameters"; aren't they part of the fitted model already? Can you be more clear about why simulate isn't what you need? (PS. You might also look at sim in the arm package which simulates the uncertainty in the fitted model as well.) –  Aaron Feb 17 '12 at 1:52
I think Jonathan wants to pre-specify the values of the covariates for simulation. If you set x to c(1,2,3) and mod<-lm(c(1,2,3)~x) then simulate(mod) always gives 1,2,3 and it would be nice if it would be possible to pre-specify x as c(7,9,23) and get 7, 9, 23 as result of the simulation but leaving the mod-object as it is. Jonathan, is this your question? –  psj Feb 17 '12 at 15:55
@psj Yes, I want to be able to specify the covariates / independent variables, but I still want the error terms to be randomly drawn. I guess "simulate" is the wrong term? I'll work on clarifying my question. –  Jonathan Feb 17 '12 at 16:08
@Aaron, @psj I edited the question to make it clearer. I looked at the arm package and didn't see a way to specify the covariates (however, I didn't make it clear that that was what I was trying to do when you commented). –  Jonathan Feb 17 '12 at 16:33

Use the predict() command's newdata argument to obtain mean responses at given covariate value, then generate data using the rnorm(), rbinom(), rpois() function appropriate for your GLM. For Gaussian data you'll also need to specify a dispersion parameter, but this is part of the glm() output, or its summary().

-
How would I get the standard deviation for these functions? How would I use the dispersion? –  Jonathan Feb 17 '12 at 16:50
I guess sd(model$residuals) should be used? – Jonathan Feb 17 '12 at 17:05 That's one way, which I agree is less parametric. Being more parametric you could you use e.g. lm1 <- lm(y~x) then get an appropriate estimate of sd from summary(lm1)$sigma –  guest Feb 18 '12 at 5:06

I ended up doing the following:

model <- glm(formula,data)
pred <- predict(model, newdata)
y <- pred + sample( model\$residuals, N=dim(newdata)[1] )


The benefit of this over rnorm(), etc. is that it is less parametric.

-