0
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • $\begingroup$ @Glen_b no, need this for work. $\endgroup$ – spore234 Oct 12 '15 at 7:41
5
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ as far as I understand this you describe how I can generate a gamma outcome given some predictors X. But I want the predictors X given some gamma outcome. $\endgroup$ – spore234 Oct 12 '15 at 8:02
  • 1
    $\begingroup$ I'm sorry, but that's not how GLMs work (or indeed, regression models). GLMs can't have constant mean response with significant predictors. The mean of the gamma distribution is supposed to vary with the predictors -- that's really the point of a GLM (to estimate how the mean varies as a function of predictors). Unless you set up your predictors so they're perfectly collinear (in such a way that $X\beta$ is constant, and in which case you can't estimate the model you used to generate the data), then it's literally impossible to make the mean constant. You won't be able to have significance. $\endgroup$ – Glen_b Oct 12 '15 at 8:07
  • 1
    $\begingroup$ More to that, why will you want to generate predictors given a response?. $\endgroup$ – Chamberlain Foncha Oct 12 '15 at 8:13
  • $\begingroup$ @ChamberlainFoncha I have some real data where the outcome is gamma and I encountered some problems. I want to see what changes when the covariables are more predictive. $\endgroup$ – spore234 Oct 12 '15 at 9:45
  • $\begingroup$ @Glen_b is this still the case when I exchange GLM with Random Forest? $\endgroup$ – spore234 Oct 12 '15 at 9:46
2
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ $y$ generated in that fashion will be a shifted gamma, which is not what you model with a GLM. $\endgroup$ – Glen_b Oct 12 '15 at 8:02
  • $\begingroup$ The shift is absorbed in the intercept. $\endgroup$ – Chamberlain Foncha Oct 12 '15 at 8:14
  • $\begingroup$ Not in a gamma glm it isn't, and you're not shifting them by 2, but by 2+3*x1+4*x2+5*x3. Note, for example, (i) gamma variates go from 0 to infinity; yours doesn't-- it goes from $X\beta$ to infinity; (ii) with a gamma glm the variance is proportional to the square of the mean, but if you shift your gammas, you move the mean but the variance doesn't change. You're generating data that doesn't follow the GLM model. If you were generating an ordinary regression model that approach works, but in GLMs only the gaussian case is a location family (i.e. amenable to just being shifted) $\endgroup$ – Glen_b Oct 12 '15 at 8:16
  • $\begingroup$ If it's still not clear, consider what happens when you use this strategy with a logistic regression model. Many of the underlying issues there apply to essentially every other GLM but the Gaussian case. However, if you use a log-link with the gamma, you can adapt your approach - add the log of the pre-existing gamma to $X\beta$ and then exponentiate; in that case it all works, with only the intercept being shifted from the one used with $X\beta$. $\endgroup$ – Glen_b Oct 12 '15 at 8:42
  • $\begingroup$ @Glen_b as I want to use a log link this solution sounds very plausible and easy to adapt. $\endgroup$ – spore234 Oct 12 '15 at 9:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.