Skip to main content
1 of 3
enfield
  • 101
  • 1

Help with R implementation

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, my 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 iteration I get NaN values for the coefficients. I have tried re-deriving $w_i^{(m-1)}$ and $z_i^{(m-1)}$, but I always get the same. As you can see, I have a GLM at the end of the code, and runs perfectly fine. 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 appreaciated!

enfield
  • 101
  • 1