the predictor is equal to (in the code case we don't have the intercept): $\eta_i = \sum_{j=1}^{2}\beta_jx_{ij}=\beta_1x_{i1}+\beta_{i2}$$\eta_i = \sum_{j=1}^{2}\beta_jx_{ij}=\beta_1x_{i1}+\beta_{i2}x_{i2}$
As stated in the first link above $W$ is a diagonal matrix, where each element of the diagonal is the second partial derivative in respect of the vector of parameters $\beta$ of fitted values of the Logistic Regression
the residual $z =\frac{y_i - E[y_i]}{h'(\eta_i)}$ where $h'(\eta_n)$ is the first partial derivative of the fitted values in respect of the vector of the same parameters, and it is equal to $h'(\eta) = \frac{1}{1+e^\eta}*(1-\frac{1}{1+e^\eta})$
#LOGISTIC REGRESSION Estimation (IRLS)
#LOGIT
set.seed(5)
p <- 2 ##per##for p > 3 questothe algoritmoestimates èdo nonnot consistenteconverge
n <- 20
x <- matrix(rnorm(n * p), n, p)
betas <- runif(p, -2, 2)
hc <- function(x) 1 /(1 + exp(-x)) # inverse canonical link
p.true <- hc(x %*% betas)
y <- rbinom(n, 1, p.true)
tol=1e-8
#IRLS using the 'lm' function:
b.init = rep(1,p)
b.old <- b.init
change <- Inf
IRLS_canoni_ = list()
while(change > tol) {
eta <- x %*% b.old # linear predictor
y.hat <- hc(eta)
h.prime_eta <- y.hat * (1 - y.hat) #first derivative
z <- (y - y.hat) / h.prime_eta
b.new <- b.old + lm(z ~ x - 1, weights = h.prime_eta)$coef # WLS regression
change <- sqrt(sum((b.new - b.old)^2))
b.old <- b.new
IRLS_canoni_$eta = cbind(IRLS_canoni_$eta,eta)
IRLS_canoni_$y.hat = cbind(IRLS_canoni_$y.hat,y.hat)
IRLS_canoni_$h.prime_eta = cbind(IRLS_canoni_$h.prime_eta, h.prime_eta)
IRLS_canoni_$z = cbind(IRLS_canoni_$z, z)
IRLS_canoni_$b.old = cbind(IRLS_canoni_$b.old, b.old)
print(b.old)
Sys.sleep(.1)
}
b.old
my_IRLS_canonical(x, y, rep(1,p), hc)
glm(y ~ x - 1, family=binomial())$coef #model with no intercept
glm1 = glm(y ~ x, family=binomial())
##Trying to obtain same results with matrix notation (IRLS):
deriv2 = function(x) exp(x)/(1+exp(x))^2 #second derivative
b.init = rep(1,p)
b.old1 <- b.init
change1 <- Inf
IRLS_matrix = list()
while(change1 > tol) {
eta1 <- x %*% b.old1 # linear predictor
y.hat1 <- hc(eta1)
h.prime_eta1 <- y.hat1 * (1 - y.hat1)
z1 <- (y - y.hat1) / h.prime_eta1
Wdiag = deriv2(eta)
W = matrix(0,n,n)
diag(W) = Wdiag
H = -(t(x)%*%(W)%*%x) #not using it
b.new1 = b.old1 + ((solve(t(x) %*% W %*% x)) %*% (t(x)%*%W%*%z1))
change1 = sqrt(sum((b.new1 - b.old1)^2))
b.old1 = b.new1
IRLS_matrix$eta = cbind(IRLS_matrix$eta, eta1)
IRLS_matrix$y.hat = cbind(IRLS_matrix$y.hat, y.hat1)
IRLS_matrix$h.prime_eta = cbind(IRLS_matrix$h.prime_eta, h.prime_eta1)
IRLS_matrix$z = cbind(IRLS_matrix$z, z1)
IRLS_matrix$b.old = cbind(IRLS_matrix$b.old, b.old1)
print(b.new1)
Sys.sleep(.1)
}
b.new1
glm(y ~ x - 1, family=binomial())$coef #model with no intercept
IRLS_canoni_$eta[,1] == IRLS_matrix$eta[,1]
IRLS_canoni_$y.hat[,1] == IRLS_matrix$y.hat[,1]
IRLS_canoni_$h.prime_eta[,1] == IRLS_matrix$h.prime_eta[,1]
IRLS_canoni_$z[,1] == IRLS_matrix$z[,1]
IRLS_canoni_$b.old[,1] == IRLS_matrix$b.old[,1]