I believe this is a statistics problem relating to my misunderstanding GLM's. But there's a chance it's a programming problem. If that turns out to be the case, I'll move over to Stack Overflow.
I've simulated data in a way that I believe would be perfectly captured by a gaussian GLM with the inverse link. But when I give the data to
stats::glm, I'm getting inaccurate estimates (or no convergence).
Here's what I mean:
set.seed(27599) n <- 1e4 x1 <- runif(n = n) x2 <- runif(n = n) linpred <- 0.5*x1 + -0.5*x2 mu <- 1/linpred y <- rnorm(n = n, mean = mu, sd = 1) d <- data.frame(x1, x2, y) fitted_glm <- glm(formula = y ~ x1 + x2, data = d, family = stats::gaussian(link = 'inverse'), mustart = rep(mean(y), n)) coef(fitted_glm) > (Intercept) x1 x2 > 19247667 -29105253 -16315984
I would expect (hope) that the intercept would be near
0 and the effects would be near
One thing I tried is applying the inverse directly to
y, rather than to the linear predictor, as in:
mu <- mean_lp y <- 1/rnorm(n = n, mean = mu, sd = 1) coef(fitted_glm) > (Intercept) x1 x2 > -1.902150e+12 2.745787e+12 -3.221891e+08
But, as you can see, that also gave nonsensical parameter estimates. Where have I gone wrong? I used a similar process to simulate and model a normal-identity and normal-log GLM and it worked as expected.