How can I generate data for a GLM that explains my outcome well?

Question

I want to test several glm methods of an outcome that follows a gamma distribution.

I can generate this outcome like this:

y <- rgamma(100, shape=0.5, rate=1)

I now need three variables x1, x2 and x3 that explain these data in a very good way. They should be equally significant and not one that explains everything.

How can I do this?

I only know R. It should also be possible to do this for other continuous distributions for y.

Answer 1

For a GLM you want the mean of the response, $\mu$ to change as a function of the predictors.

See the overview of GLMs at Wikipedia, which says:

$\operatorname{E}(\mathbf{Y}) = \boldsymbol{\mu} = g^{-1}(\mathbf{X}\boldsymbol{\beta})$

For a gamma distribution, the mean is shape/rate or equivalently shape x scale. The shape parameter is help constant in a gamma GLM. (So in this case, $\mu$ is proportional to scale or inversely proportional to rate.)

That is, you'll have some link function $g^{-1}$, such that $g^{-1}(\mu_i) = \eta_i = x_i\beta$, where $x_i$ is a row-vector from your matrix of predictors.

So construct some betas (including for the intercept), and stick your predictors together into an X-matrix. You'll also need to choose a shape parameter for your gamma (if you don't have any experience of what sort of gamma you might need, something in the vicinity of 2-3 is usually skew enough to be interestingly non-normal without having the mode at zero). I call that shape parameter gshape below.

• So you could calculate, for example, eta = beta0 + X %*% beta,
• then mu=g(eta), $\:$ (e.g. exp for a log-link)
• then gscale = mu/gshape,
• then call rgamma with scale=gscale or rate=1/gscale (they're the same)
  (note that the n in the call must be the same as the length of the scale parameter.)

You may have to try several betas before they have about the properties you want.

The same approach works for other glm models, and in other packages (with appropriate substitutions of code).

The usual way to achieve this would be to make the scale parameter (or the inverse of the rate parameter) depend on those predictor variables.

Answer 2

A very simple model that assumes that $y$ is related $x_1$, $x_2$ and $x_3$ in a linear fashion, with no interactions is this

x1=runif(100) x2=runif(100) x3=runif(100) x4=runif(100) y <- 2+3*x1+4*x2+5*x3 +rgamma(100, shape=0.5, rate=1) The random variable $y$ now has gamma distribution with $E(y)=0.5+2+3x_1+4x_2+5x_3$ and $Var(y)=0.5$. Hence $y$ has shape, $\alpha$, $\alpha=E^2(y)/Var(y)$ and rate, $\beta$, $\beta=E(y)/Var(y)$. This is a very simple linear model and a good place to start. Of course the values multiplying $x$ you can change them as you wish.