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