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

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,
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.

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.)
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.
• I see that this starting value caused problems in this case. In general, though, mean(y) seems like a reasonable starting value to me. Would you agree? – rcorty Sep 21 '17 at 14:03
• In some cases, mean(y) may be a good initial value, but here $y$ has a huge variance (var(y) = 590846.6, that's because $y$ shoots off to infinity, as Glen_b mentions). In other words, mean(y) is pretty far away from the true value of $\mu$ for most observations. It's much better to simply use $y$ as initial values. Try it out, it works... – bbrot Sep 21 '17 at 15:17
• got it -- thanks for your input. Not sure why I didn't think to use y! – rcorty Sep 21 '17 at 18:37
• @rcorty Setting mustart=y is what glm does by default, so (as bbrot says) there was never a need for you to set mustart yourself. – Gordon Smyth May 17 at 6:51