1
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – rcorty Sep 21 '17 at 14:03
  • $\begingroup$ 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... $\endgroup$ – bbrot Sep 21 '17 at 15:17
  • $\begingroup$ Usually, I think there's no need to manually set the initial values. R is quite clever about that. $\endgroup$ – bbrot Sep 21 '17 at 15:31
  • $\begingroup$ got it -- thanks for your input. Not sure why I didn't think to use y! $\endgroup$ – rcorty Sep 21 '17 at 18:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.