6
$\begingroup$

I'm trying to take the approach for understanding how certain concepts work, by trying to generate data for them and checking how the output behaves. Currently, I thus realized I don't quite get what's going on with GLM-s.

Here is my little code:

N = 10000
e = rnorm(N,0,1)
x1 = runif(N,10,30)
y = exp(5*x1+ 10 + e)
mod1 = glm (y ~ x1,family=gaussian(link="log"))
mod2 = lm(log(y) ~ x1)

Calling summary of the models quickly reveals, that mod2 is a nice fit, while mod1 is bonkers. I tried brushing up on the topic, and many pages talk of transforming the mean of y, because y is not normally distributed, but I never really understood this, since the assumption of normality is for the residuals, not $y$, which is quite logical otherwise $y = mx +b $ would never work with $x$ sampled uniformly.

So I have two questions:

  • What am I not getting here?
  • How would I generate data that is valid for the above GLM?

EDIT

I reformulated the code, to reflect more closely the mathematical background (based on Glen_b's answer I realized my way of adding the error doesn't work for all cases).

x = seq(from = 1,to = 15,by = 0.1)
N = length(x)
eta = 5*x + 10

# original
set.seed(5671)

y = exp(eta) + rnorm(N,0,1)
mod = glm(y ~ x,gaussian(link = "log"))

# new
set.seed(5671)

inverse_link = function(x){exp(x)}
means = sapply(eta,function(x){inverse_link(x)})
y = sapply(means,function(x){rnorm(1,mean=x,sd=1)})
mod = glm(y ~ x,gaussian(link = "log"))

The results can be compared to be the same in both cases. Based on this my expectation was that the following code can fit my parameters properly:

x = seq(from = 1,to = 15,by = 0.1)
N = length(x)
eta = 5*x + 10

set.seed(5671)

inverse_link = function(x){1/x}
scale = 1
shapes = sapply(eta,function(x){inverse_link(x)/scale})

y = sapply(shapes,function(x){rgamma(1,shape=x,scale=scale)})
mod = glm (y ~ x,family=Gamma(link="inverse"))

My reasoning was that $\mu = k\theta = 1/\eta$ so I need a gamma distribution with shape parameter $k = 1/\eta/\theta$. My problem is that this was, the coefficients are wildly off (~1.6 for $x$ and ~6.4 for intercept). Is it just my input data, or did I miss something?

EDIT 2

As pointed out in the gamma distribution the $k$, the shape parameter is kept constant (as far as I understand GLM assumes that the distributions used are from the natural exponential family, based on this answer and gamma is that with a fixed shape parameter). So we have

$$y \sim \Gamma(k,\frac{1}{k\eta})$$

Here is the corrected code that now works:

x = seq(from = 1,to = 15,by = 0.1)
N = length(x)
eta = 5*x + 10

set.seed(5671)

inverse_link = function(x){1/x}
shape = 3
scales = sapply(eta,function(x){inverse_link(x)/shape})

y = sapply(scales,function(x){rgamma(1,shape=shape,scale=x)})
mod = glm (y ~ x,family=Gamma(link="inverse"))
summary(mod)
```
$\endgroup$

2 Answers 2

2
$\begingroup$

Here's how to generate from a glm (the order of some items can be moved):

  1. Choose your family and link function.

  2. choose your predictors (IV's) for each observation you want to simulate.

  3. Choose your coefficients.

  4. Evaluate the linear predictor for each observation.

  5. Transform by the inverse of the link function to get the conditional mean for each observation.

  6. Choose any other parameters.

  7. Sample the distribution at each observation, for which you now have all the parameters.


Let's see how to simulate a simple Gamma GLM with inverse link, following those steps:

  1. Choose your family and link function. (Gamma, inverse)

  2. choose your predictors (IV's) for each observation you want to simulate. ($x$)

  3. Choose your coefficients. (choose a specific $\beta_0$ & $\beta_1$ in this case)

  4. Evaluate the linear predictor for each observation. ($\eta_i=\beta_0+\beta_1 x_i,\: i=1,...,n$)

  5. Transform by the inverse of the link function to get the conditional mean for each observation. (inverse of $\eta_i=1/\mu_i$ is $\mu_i=1/\eta_i$)

  6. Choose any other parameters. (Choose the shape parameter)

  7. Sample the distribution at each observation, for which you now have all the parameters. (e.g. y=rgamma(length(x),shape,scale=mu/shape) -- noting that scale is a vector of values here)

$\endgroup$
5
  • $\begingroup$ Your answer prompted me to reformulate the code to reflect the mathematics more closely (and also because I think my original formulation is pretty much restricted to the Gaussian family). I tried applying it to a gamma distribution though and it failed miserably :( $\endgroup$
    – fbence
    Nov 20, 2019 at 7:52
  • $\begingroup$ Note that the parameterization you're using in a Gamma glm is actually a mean-shape parameterization and the shape should be identical for every observation (the mean and hence the scale vary). You may need to review how glms work. $\endgroup$
    – Glen_b
    Nov 20, 2019 at 11:26
  • $\begingroup$ I am accepting this answer, because although it didn't actually correct the code, it explained the general process, while my original formulation was a bit restricted. $\endgroup$
    – fbence
    Nov 20, 2019 at 17:14
  • $\begingroup$ I couldn't simply correct your code since you were trying to generate a model by varying the shape parameter and fixing the scale. Given the need to model the mean and fix the shape parameter, you needed to completely reformulate your model; I didn't want to choose your new model for you. I have written a more detailed recipe for a gamma glm with inverse link. By sticking (mostly) to not explicitly writing R code it should help people using other programs than R. $\endgroup$
    – Glen_b
    Nov 21, 2019 at 1:04
  • $\begingroup$ Thank you, I agree, this way it is more useful to other people. $\endgroup$
    – fbence
    Nov 21, 2019 at 13:39
3
$\begingroup$

It matters whether the error term is included in the exp() call or not. That's your big issue for this code. Consider this:

#  first, set the random seed, so that everything is reproducible
set.seed(5671)
N  = 10000
e  = rnorm(N,0,1)
x1 = runif(N,10,30)
y1 = exp(5*x1+ 10  + e)
y2 = exp(5*x1+ 10) + e
mod1.1 = glm(y1 ~ x1,family=gaussian(link="log"))
mod1.2 = lm(log(y1) ~ x1)

mod2.1 = glm(y2 ~ x1,family=gaussian(link="log"))
mod2.2 = lm(log(y2) ~ x1)
$\endgroup$
6
  • $\begingroup$ When using y = exp(5*x1+ 10) + e, then the glm will not even converge after 25 iterations, while mod2 is just bad (the latter I understand, since the residuals will not be normally distributed). $\endgroup$
    – fbence
    Nov 19, 2019 at 20:17
  • $\begingroup$ Sorry, yeah setting the seed is quite crucial here :) Only looking at mod1.1 and mod2.1, if I check the deviance residuals, in both cases I get insane numbers, e.g for the medians respectively: -1.116e+45 and 1.393e+34 which seems to suggest something is still wrong (or I don't really understand what these numbers are). Although I have to agree, the fit seems perfect in the second case. $\endgroup$
    – fbence
    Nov 19, 2019 at 20:27
  • $\begingroup$ @fbence, you're right about the convergence. It's fitting fine WRT the coefficients, but it's having trouble with the dispersion parameter (the residual SD). I'm not sure what the issue is. It looks like we need to use a different optimizer. $\endgroup$ Nov 19, 2019 at 20:38
  • $\begingroup$ What I don't really understand, is that whatever optimizer we use, calculating the error based on the good fitted parameters (which is the case here, since we have 5 and 9.9999 for the $\beta$ and $b$) should be very small, no? $\endgroup$
    – fbence
    Nov 19, 2019 at 21:34
  • $\begingroup$ I don't follow you, @fbence. The error between the fitted & true parameter values are very small. $\endgroup$ Nov 19, 2019 at 21:37

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.