Iteratively Re-weighted Least Squares (IRWLS)
S.S.B.1
setwd("C:/ssb")List of libraries used in this project
library(bookdown)
library(kableExtra)1 Algorithm for IRWLS
Algorithm
2 Algorithm for IRWLS Logistic Regression
Algorithm
datta <- read.csv("data/MichelinNY.csv", header=TRUE)
dat <- head(datta, 10)
kab <- knitr::kable(dat, caption = "Sample dataset",
booktabs = T, label = "kabletable")
kable_classic_2(kab, full_width = F)| InMichelin | Restaurant.Name | Food | Decor | Service | Price |
|---|---|---|---|---|---|
| 0 | 14 Wall Street | 19 | 20 | 19 | 50 |
| 0 | 212 | 17 | 17 | 16 | 43 |
| 0 | 26 Seats | 23 | 17 | 21 | 35 |
| 1 | 44 | 19 | 23 | 16 | 52 |
| 0 | A | 23 | 12 | 19 | 24 |
| 0 | A.O.C. | 18 | 17 | 17 | 36 |
| 1 | A.O.C. Bedford | 24 | 21 | 22 | 51 |
| 1 | Aix | 23 | 22 | 21 | 61 |
| 1 | Alain Ducasse | 27 | 27 | 27 | 179 |
| 0 | Alouette | 20 | 17 | 19 | 42 |
y <- datta$InMichelin
x1 <- datta$Decor
x2 <- datta$Food
x3 <- datta$Service
goal <- glm(y~x1+x2+x3,family=binomial)
goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 # Start with mu near y
epsilon = 0.1
mu0 = abs(y-epsilon)
eta0 = log(mu0/(1-mu0))
z0 = eta0 + (y-mu0)/(mu0*(1-mu0))
w0 = mu0*(1-mu0)
mod = lm(z0~x1+x2+x3,weights=w0)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -13.90422861 0.31978636 0.38577436 -0.03768687 ################ Iterate this #########################
eta = b0 + b1*x1 + b2*x2 + b3*x3
mu = exp(eta)/(1+exp(eta))
z = eta + (y-mu)/(mu*(1-mu))
w = mu*(1-mu)
mod = lm(z~x1+x2+x3,weights=w)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -1.243707e+01 2.604819e-01 3.425819e-01 9.356367e-04 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 ################ Iterate this #########################
eta = b0 + b1*x1 + b2*x2 + b3*x3
mu = exp(eta)/(1+exp(eta))
z = eta + (y-mu)/(mu*(1-mu))
w = mu*(1-mu)
mod = lm(z~x1+x2+x3,weights=w)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -12.623128368 0.267148400 0.348969062 -0.003730635 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 ################ Iterate this #########################
eta = b0 + b1*x1 + b2*x2 + b3*x3
mu = exp(eta)/(1+exp(eta))
z = eta + (y-mu)/(mu*(1-mu))
w = mu*(1-mu)
mod = lm(z~x1+x2+x3,weights=w)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -12.625967033 0.267238762 0.349074585 -0.003796467 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 ################ Iterate this #########################
eta = b0 + b1*x1 + b2*x2 + b3*x3
mu = exp(eta)/(1+exp(eta))
z = eta + (y-mu)/(mu*(1-mu))
w = mu*(1-mu)
mod = lm(z~x1+x2+x3,weights=w)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 #######################################################
# Try another way to start
mod = lm(y~x1+x2+x3)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -1.601393257 0.048330398 0.058303388 -0.005695744 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 ################ Iterate this #########################
eta = b0 + b1*x1 + b2*x2 + b3*x3
mu = exp(eta)/(1+exp(eta))
z = eta + (y-mu)/(mu*(1-mu))
w = mu*(1-mu)
mod = lm(z~x1+x2+x3,weights=w)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -9.16038062 0.20548859 0.24746897 -0.00936594 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 ################ Iterate this #########################
eta = b0 + b1*x1 + b2*x2 + b3*x3
mu = exp(eta)/(1+exp(eta))
z = eta + (y-mu)/(mu*(1-mu))
w = mu*(1-mu)
mod = lm(z~x1+x2+x3,weights=w)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -11.98873619 0.25686296 0.33078995 -0.00660303 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 ################ Iterate this #########################
eta = b0 + b1*x1 + b2*x2 + b3*x3
mu = exp(eta)/(1+exp(eta))
z = eta + (y-mu)/(mu*(1-mu))
w = mu*(1-mu)
mod = lm(z~x1+x2+x3,weights=w)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -12.601590571 0.266867827 0.348389976 -0.003956219 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479 ################ Iterate this #########################
eta = b0 + b1*x1 + b2*x2 + b3*x3
mu = exp(eta)/(1+exp(eta))
z = eta + (y-mu)/(mu*(1-mu))
w = mu*(1-mu)
mod = lm(z~x1+x2+x3,weights=w)
b0 = mod$coefficients[1]; b1 = mod$coefficients[2]
b2 = mod$coefficients[3]; b3 = mod$coefficients[4]
c(b0,b1,b2,b3)## (Intercept) x1 x2 x3
## -12.625930288 0.267238236 0.349073584 -0.003796788 goal$coefficients## (Intercept) x1 x2 x3
## -12.625967612 0.267238778 0.349074607 -0.003796479U of Universe, ssb@universe↩︎