Solved – How to simulate new observations to a fitted model object while specifying the covariates in R

rsimulation

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?

Best Answer

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().

Related Solutions

Solved – simulate GLM with square root link in R

I figured out some answers to my questions.

Regarding the mathematical expressions of predict and simulate, these can be obtained by the following code (thanks to a tip from W. van der Elst):

getS3method(c("predict"), class = "glm")
getS3method(c("simulate"), class = "lm")

Regarding the predict function, I incorrectly used the option type=”response”. This does not include the stochastic uncertainty as was my aim. It is the prediction on the backtransformed linear scale. This can be tested by the following graph (run after the initial code stated in the question):

plot(predict(m, type="response"), p_lin^2)

Regarding the simulate function, it seemed to be correct if I look at the plot from the initial question:

plot(p, d[,1], xlab="predicted values", ylab="original data", xlim=xylim, ylim=xylim, col=rgb(0,0,0,alpha=0.1))
points(simulate(m)$sim_1, d[,1], col=rgb(0,1,0,alpha=0.1))

However, if I upscale the dependent variable:

d[,1] <- d[,1]*10000

And recalculate the predictions:

m <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4], family=gaussian(link="sqrt"))
p_lin <- m$coef[1] + m$coef[2]*d[,2] + m$coef[3]*d[,3] + m$coef[4]*d[,4]
p <- rnorm(n=n, mean=p_lin^2, sd=sd(p_lin^2 - d[,1]))

And plot the results:

par(mfrow=c(1,1), mar=c(4,2,2,1), pch=16, cex=0.8, pty="s")
xylim <- c(min(c(d[,1], p)), max(c(d[,1], p)))
plot(p, d[,1], xlab="predicted values", ylab="original data", xlim=xylim, ylim=xylim, col=rgb(0,0,0,alpha=0.1))
points(simulate(m)$sim_1, d[,1], col=rgb(0,1,0,alpha=0.1))

I see different predictions.

This seems to be attributable to a difference between the value of the sd in my general formula approach and the simulate() function approach. Recall the calculation of the sd in the general formulas approach:

sd=sd(p_lin^2 - d[,1])

in the simulate function R calculated the sd differently from the general formulas approach (for a Gaussian family):

vars <- deviance(m)/df.residual(m)
if (!is.null(m$weights)) vars <- vars/m$weights # the m$weights seems similar to the m$fitted.values multiplied by about 4
fitted(m) + rnorm(n, sd = sqrt(vars))

I don’t understand why the sd is calculated and then divided by the m$weights. Why is sd a vector and not a single value? The simulate() help text states: “The methods for linear models fitted by lm or glm(family = "gaussian") assume that any weights which have been supplied are inversely proportional to the error variance.” I can’t seem to grasp the meaning of this. If I run multiple simulations using the simulate() function, the simulations look very much like each other:

plot(simulate(m)$sim_1,simulate(m,2)$sim_2)

If I do not use a sqrt-link function I get simulations that seem to better reflect stochastic uncertainty as they are less identical when I run them twice:

library(MASS)
rm(list=ls())
set.seed(2)
n <- 1500
d <- mvrnorm(n=n, mu=c(0,0,0,0), Sigma=matrix(.7, nrow=4, ncol=4) + diag(4)*.3)
d[,1] <- qgamma(p=pnorm(q=d[,1]), shape=2, rate=2)
d[,1] <- d[,1]*10000
m <- glm(formula=d[,1] ~ d[,2] + d[,3] + d[,4])
plot(simulate(m)$sim_1,simulate(m,2)$sim_2)

What method is correct? What is the difference? (should I post this as a separate question?)

Solved – How to predict a distribution (from a set of predictors) that I can simulate from

You can do this with linear regression! Whether you use the MLE or minimize the mean squared error, you are actually modeling the mean response, conditional on the predictors, and you are implicitly assuming that the response itself is normally distributed around that mean. For other distributions, you would use a generalized linear model.

What you're looking for also arises naturally as the "posterior predictive distribution" in a Bayesian model. It might be helpful to think of maximum-likelihood-based inference as Bayesian inference with flat priors.

Fully nonparametric density estimation is very difficult in general except in low-dimensional cases. I asked an involved question about it about a year ago, but I can't find it now.

Best Answer

Related Solutions

Solved – simulate GLM with square root link in R

Solved – How to predict a distribution (from a set of predictors) that I can simulate from

Related Question