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. Now these could very well be wrong, but I like to think not. Now, myhave implemented this in R implementation is as follows:
library(glmx)
library(msme)
data("BeetleMortality")
data <- as.data.frame(BeetleMortality)
# initializing variables & coefficients
x <- data$dose
X <- cbind(1, x) # design matrix with intercept
y <- data$died
n <- data$n
b <- rep(0, ncol(X))
print(X)
# IRWLS algorithm
for (iter in 1:5) {
eta <- X %*% b
mu <- pnorm(eta) # fitted probabilities CDF
# computing weights
phi_eta <- dnorm(eta) # ftandard normal PDF at eta
w <- (phi_eta^2) / (mu * (1 - mu))
# computing working response
z <- eta + (y - mu) / (phi_eta)
# updating coefficients
W <- diag(as.vector(w))
b_new <- solve(t(X) %*% W %*% X) %*% (t(X) %*% W %*% z)
b <- b_new # update coefficients
}
cat("Estimated coefficients (IRWLS):\n")
print(b)
glm_fit <- glm(cbind(y, n - y) ~ dose, data = BeetleMortality, family = binomial(link =
probit))
cat("Estimated coefficients (glm):\n")
print(coef(glm_fit))
The code runs, but after only one iterationwhere I get NaN values foram pretty sure everything in the coefficientsimplementation is correct. I have tried re-deriving $w_i^{(m-1)}$ and $z_i^{(m-1)}$, butTherefore I always getam lead to belief that the same. As you can see,problem is in how I have a GLM at the end of the code, and runs perfectly finederived 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 appreaciatedappreciated!