I have an exercise where I have to derive the both $w_i^{(m-1)}$ and $z_i^{(m-1)}$ from the iterative weighted least squared updating equation $b^{(m)} = \left( X^\top W^{(m-1)} X \right)^{-1} X^\top W^{(m-1)} z^{(m-1)}$ for the BeetleMortality data with the probit link $\phi$. Then I have to implement it in R. From what I understand, $w_i^{(m-1)}=\frac{1}{\text{Var}(Y)} (\frac{\partial \mu_i}{\partial \eta_i}) ^2 = \frac{\phi'(\eta_i)^{(m-1)}}{n\mu_i^{(m-1)}(1 - \mu_i^{(m-1)})}$ and $z_i^{(m-1)} = \eta_i^{(m-1)} (y_i -\mu_i) \frac{\partial \eta_i}{\partial \mu_i}= \eta_i^{(m-1)} + \frac{y_i + \mu_i^{(m-1)}}{\phi'(\eta_i)^{(m-1)}}$ where $\phi'$ is the normal distribution PDF. I have implemented this in R, where I am pretty sure everything in the implementation is correct. Therefore I am lead to belief that the problem is in how I have derived these values. Did I make any mistake? I am relatively new to all this so it could very well be that I made some glaring mistake which I cannot find. Any feedback is appreciated!
1 Answer
Yes, there are mistakes in your derivations. The main issues are:
- You used the derivative of the normal PDF ($\phi'(\eta)$) instead of the PDF itself ($\phi(\eta)$) when expressing the weights.
- In the expression for $z_i^{(m-1)}$, you should have $(y_i - \mu_i^{(m-1)})$ in the numerator, not $(y_i + \mu_i^{(m-1)})$. Also, you should use $\phi(\eta_i)$, not $\phi'(\eta_i)$.
In fitting a generalized linear model (GLM) with a given link function using Iteratively Reweighted Least Squares (IRLS), the update step for the coefficient vector $b$ is often given by:
$$ b^{(m)} = (X^\top W^{(m-1)} X)^{-1} X^\top W^{(m-1)} z^{(m-1)} $$
Here, $W^{(m-1)}$ is the diagonal matrix of working weights and $z^{(m-1)}$ is the working response vector at iteration $(m-1)$. For a binary/bionomial model with responses $Y_i$ (often successes out of $n$ trials) and a probit link, we have:
- The link: $\eta_i = g(\mu_i) = \Phi^{-1}(\mu_i)$ where $\Phi$ is the standard normal CDF, so $\mu_i = \Phi(\eta_i)$
- The derivative: $\frac{d \mu_i}{d \eta_i} = \phi(\eta_i)$, where $\phi$ is the standard normal PDF: $\phi(\eta) = \frac{1}{\sqrt{2\pi}}e^{-\eta^2/2}$
For a binomial model with $n$ trials and probability $\mu_i$, the variance is:
$$ \text{Var}(Y_i) = n \mu_i (1 - \mu_i) $$
The standard IRLS formula for weights in a GLM is:
$$ w_i^{(m-1)} = \frac{\left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2}{\text{Var}(Y_i)} $$
Plugging in the probit link derivatives:
$$ \frac{\partial \mu_i}{\partial \eta_i} = \phi(\eta_i) $$
Thus:
$$ w_i^{(m-1)} = \frac{\phi(\eta_i^{(m-1)})^2}{n \mu_i^{(m-1)} (1 - \mu_i^{(m-1)})} $$
You wrote:
$$ w_i^{(m-1)} = \frac{\phi'(\eta_i^{(m-1)})}{n \mu_i^{(m-1)} (1 - \mu_i^{(m-1)})} $$
This is incorrect. You should not be using $\phi'(\eta)$ (the derivative of the PDF), but $\phi(\eta)$ (the PDF itself). The correct form is:
$$ w_i^{(m-1)} = \frac{\phi(\eta_i^{(m-1)})^2}{n \mu_i^{(m-1)} (1 - \mu_i^{(m-1)})} $$
The working response in IRLS is generally given by:
$$ z_i^{(m-1)} = \eta_i^{(m-1)} + \frac{y_i - \mu_i^{(m-1)}}{\frac{\partial \mu_i}{\partial \eta_i}} $$
Since $\frac{\partial \mu_i}{\partial \eta_i} = \phi(\eta_i^{(m-1)})$, we have:
$$ z_i^{(m-1)} = \eta_i^{(m-1)} + \frac{y_i - \mu_i^{(m-1)}}{\phi(\eta_i^{(m-1)})} $$
You wrote:
$$ z_i^{(m-1)} = \eta_i^{(m-1)} + \frac{y_i + \mu_i^{(m-1)}}{\phi'(\eta_i^{(m-1)})} $$
Two errors here:
- It should be $y_i - \mu_i^{(m-1)}$, not $y_i + \mu_i^{(m-1)}$
- It should be divided by $\phi(\eta_i^{(m-1)})$ not $\phi'(\eta_i^{(m-1)})$
The corrected form is:
$$ z_i^{(m-1)} = \eta_i^{(m-1)} + \frac{y_i - \mu_i^{(m-1)}}{\phi(\eta_i^{(m-1)})} $$
