getting inaccurate estimates from glm with family = gaussian(link = 'inverse')

Question

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 0.5 and -0.5.

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.

Answer 1

Note that if the range of your linear predictor goes across 0, its reciprocal is going to shoot off toward $\pm\infty$ there; you would expect some difficulties with convergence in that case.

See how you go if you make it 1+0.5*x1-0.5*x2 instead. (With the reciprocal link, normally your means should be all-positive. That's not to say it can't work when it does cross 0 -bbrot is right that it converges when you remove the start values- but it's making things hard for yourself and it may not always behave so well.)

Answer 2

This is a numerical problem due to unfortunate initial values for $\mu$. Remove the mustart-option from the glm-call and everything works as expected.