Understanding glm and link functions: how to generate data?

Question

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)
```

Answer 1

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)

Answer 2

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)