Starting coefficient vector for GLM

Question

I would like to know how R chooses its starting coefficient vector for a GLM when its start argument is left blank and defaults to NULL. For my personal implementation of a GLM, I have simply initialized $ \boldsymbol \beta_0 $ to be all $1$s. However, while this generally is fine, it can cause the iterative algorithm to diverge.

Basically, I am just looking for a simple algorithm/formula that takes into consideration the data points and family of the GLM to choose the original coefficient vector, $ \boldsymbol \beta_0 $.

Answer 1

R's glm does not (by default) start with an initial value for $\beta$, it starts with an initial value for $\mu$. The initial value for $\mu$ depends on the family; it is close to $y$ but chosen to be in the domain of the likely link function. For example, for binomial, with $y=r/n$, $$\mu=\frac{r+1/2}{n-r+1/2}$$ and for Poisson, $\mu=y+0.1$, and for Gamma, $\mu=y$.

The initialising value for $\mu$ is used to compute the working response and working weights, and these are used to compute the first value of $\beta$ (after the first iteration)

You can specify an initial beta, and for some link/variance combinations you have to (eg, binomial(log), where the obvious $\beta=0$ doesn't work but $\beta^T=(-1,0,0,\dots,0)$ does)

Answer 2

Well, after much searching and going through papers on the theory behind GLM, I found this algorithm for the initial values, which numerically agrees with R using maxit = 1 to force R to output its initial coefficient estimates.