From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida  

 

you are kindly asked to include the complete citation if you used this material in a publication

Code 5.1 GLM logistic regression in R

==================================================


# Data
x <- c(13,10,15,9,18,22,29,13,17,11,27,21,16,14,18,8)
y <- c(1,1,1,0,0,1,0,1,1,0,0,0,1,1,0,0)

# Fit
mu <- (y + 0.5)/2                                                      # initialize mu
eta <- log(mu/(- mu))                                            # initialize eta with the Bernoulli link

for (i in 1:8) {
    w <- mu*(- mu)                                                 # variance function
    z <- eta + (y - mu)/(mu*(1-mu))                          # working response
    mod <- lm(z ~ x, weights=w)                              # weighted regression
    eta <- mod$fit                                                      # linear predictor
    mu <- 1/(1+exp(-eta))                                          # fitted value
    cat(i, coef(mod), "\n")                                          # displayed iteration log

}

# Output
summary(mod)

==================================================

Output on screen:

Call:

lm(formula = z ~ x, weights = w)

 

Weighted Residuals:

       Min               1Q      Median              3Q           Max

-1.38600     -0.93301     0.08861      0.92411      1.25563

 

Coefficients:

 

                        Estimate         Std. Error    t value        Pr(>|t|)

(Intercept)        1.28606            1.65744      0.776          0.451

x                     -0.07915            0.09701     -0.816          0.428

Residual standard error: 1.067 on 14 degrees of freedom Multiple R-squared: 0.0454, Adjusted R-squared: -0.02279 F-statistic: 0.6658 on 1 and 14 DF, p-value: 0.4282