How to draw a random sample from distribution of prediction?

Question

For my microsimulation, I want to use R to predict values and draw a random sample based on this prediction.

To clarify my point: I want to simulate the number of chronic conditions people suffer from ($y_t$) at a certain point in time. I have a few waves of panel data available to estimate a relation between age, sex, number of chronic conditions in the previous observation period (plus some others that I might include in later stages).

Suppose my regression model is $y_t = β_0 + β_1 age + β_2 sex + β_3y_{t-1} +u $.

Since R provides me with coefficients for the betas, it is easy to predict $y_t$ given the independent variables. However, this is not what I want to do. Instead, I want my population of about 1300 individuals to resemble the variance in the possible outcomes of $y_t$ (otherwise after a few steps my simulated population won’t include those unlucky ones with much more chronic conditions than the average).

I believe what I have to do is to draw a random sample from the distribution of the predicted value $y_t$, conditional on the independent variables. I further believe this can be done by drawing random numbers with mean $β$ and variance $var(β)$, multiplied by the actual values of the independent variables.

So my question is: Is this the correct approach? Will this produce reliable values? Or do I need to take possible covariation of the independent variables etc. into consideration?

Edit: Another point came to my mind. Does it make a difference whether $y_t$ or $y_t - y_{t-1}$ is my left hand side variable?

Thank you for your ideas.

Answer 1

Don't really agree with Macro. To me, it seems like what you're asking is how to perform a Bayesian analysis, in which you specify prior distributions over your $\beta_i$ and combine your observed data to obtain a posterior ("predictive") distribution, which you can sample from. This not only has benefits in terms of avoiding overfitting the regression, but also is a more natural way to handle uncertainty within your model (IMHO).

In order to do this, I'd suggest reading up on the Gibbs Sampling and Metropolis-Hastings algorithms. The basic idea is that you formulate conditional distributions over each of your parameters in terms of the other parameters in your model, and take draws from each parameter in turn. You record every $k$th observation in the chain and the samples will be drawn from the posterior distribution (thanks to some beautiful mathematics). You can use this to estimate moments, quantiles, etc.

Answer 2

I think you want to simulate from the predictor distribution values ${\rm age}_{i}$, ${\rm sex}_{i}$ and error terms $u_{i}$ and calculate

$$ (y_{t})_{i} = \hat{\beta}_{0} + \hat{\beta}_{1}{\rm age}_{i} + \hat{\beta}_{2} {\rm sex}_{i} + \hat{\beta}_{3}(y_{t-1})_{i} + u_{i} $$

to generate a sample $y_{1}, ..., y_{t}$. The predictor values can either be resampled from your data set, or generated from something similar to the empirical distribution of your predictors. The errors should be generated from the parametric distribution assumed when you fit the model, with variance estimated by the model. How you choose $(y_{1})_{i}$ is largely arbitrary.

This is essentially the same as the parametric bootstrap, except leaving off the final step where you then re-estimate the model to characterize the sampling distribution of $\hat{\beta}$, which leads me to say - I'm not completely sure why you want to do this process - if it's to see what kind of variation you can expect in the response values, I don't think this is useful for that, since I'm pretty sure the resulting variance will be about the same as the observed variance from your original data set.