Part 2 Generalized linear models (GLM): Exercise 2

1 Logistic Regression

• Logistic regression is often used to analyse the effect of a continuous or categorical explanatory( or independent) variable on a binary response variable.
• For grouped count data in Logistic regression, it assumes a binomial probability distribution function (PDF)
• For binary response data in Logistic regression, it assumes a canonical form of the Bernoulli distribution function PDF -Book Reference: Crawley (2013)

2 Objectives:

• to analyse data with logistic regression with the problem of overdispersion

3 Binary Logistic Regression

3.1 logitbin Dataset:

• independent value X
• the binary outcome 0 or 1 in the dependent variable Y

library(readxl)
my_data = read_excel("Logittest.xlsx", sheet = 1)
logitbin = as.data.frame(my_data)
logitbin

   Object  X Y
1       1  1 0
2       2  1 0
3       3  3 1
4       4  4 0
5       5  5 1
6       6  1 0
7       7  2 0
8       8  2 0
9       9  2 0
10     10  3 0
11     11  4 0
12     12  5 0
13     13  6 0
14     14  6 0
15     15  6 1
16     16  7 0
17     17  7 0
18     18  7 1
19     19  7 1
20     20  7 0
21     21  8 1
22     22  8 1
23     23  8 1
24     24  8 1
25     25  8 1
26     26  9 1
27     27  9 1
28     28  9 1
29     29 10 1
30     30 11 1
31     31 11 1
32     32 11 1
33     33 12 1
34     34 13 1
35     35 13 1
36     36  3 0
37     37  3 0

str(logitbin)

'data.frame':   37 obs. of  3 variables:
 $ Object: num  1 2 3 4 5 6 7 8 9 10 ...
 $ X     : num  1 1 3 4 5 1 2 2 2 3 ...
 $ Y     : num  0 0 1 0 1 0 0 0 0 0 ...

3.2 family=binomial(link = logit)

logitbin_model_01 = glm(Y ~ X, family = binomial(link = logit), data = logitbin)
summary(logitbin_model_01)

Call:
glm(formula = Y ~ X, family = binomial(link = logit), data = logitbin)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.51914  -0.38083   0.07704   0.60217   2.30651  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -5.1093     1.7755  -2.878  0.00401 **
X             0.8406     0.2686   3.130  0.00175 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 51.049  on 36  degrees of freedom
Residual deviance: 24.878  on 35  degrees of freedom
AIC: 28.878

Number of Fisher Scoring iterations: 6

summary(logitbin_model_01)$dispersion

[1] 1

logitbin_model_01$deviance/logitbin_model_01$df.null

[1] 0.6910642

• The Diagnostic test for logitbin_model_01 found that
  • Residual deviance = 24.878 on 35 Degrees of freedom
  • The ratio of Residual deviance/Degrees of freedom = 24.878 /36 = 0.6910642 which is lesser than 1.2, indicating that the model is not overderdispersed . That is, the variance is not increasing with the mean.
  • Transformation to quasibinomial is not necessary.

3.3 predicted values as a logistic curve: logitbin_model_01

plot(Y ~ X, cex = 1.3, cex.axis = 1.4, cex.lab = 1.4, las = 1, data = logitbin, lty = 1, 
    lwd = 3, col = "blue")
range_x = seq(min(logitbin$X), max(logitbin$X))
lines(range_x, predict(logitbin_model_01, list(X = range_x), type = "resp"), lty = 1, 
    lwd = 2, col = "red")
grid()

3.4 Basic Diagnostic Plot: logitbin_model_01

oldpar <- par(oma = c(0, 0, 3, 0), mfrow = c(2, 2))
plot(logitbin_model_01)

4 Proportional logistic regression(PLR)

4.1 logitprop Dataset:

• two dependent variables Y0 and Y1 containing the frequency of the outcome [0,1] for each level of X.

my_data = read_excel("Logittest.xlsx", sheet = 2)
logitprop = as.data.frame(my_data)
logitprop

    X Y0 Y1
1   1  3  0
2   2  3  0
3   3  3  1
4   4  2  0
5   5  1  1
6   6  2  1
7   7  3  2
8   8  0  5
9   9  0  3
10 10  0  1
11 11  0  3
12 12  0  1
13 13  0  2

str(logitprop)

'data.frame':   13 obs. of  3 variables:
 $ X : num  1 2 3 4 5 6 7 8 9 10 ...
 $ Y0: num  3 3 3 2 1 2 3 0 0 0 ...
 $ Y1: num  0 0 1 0 1 1 2 5 3 1 ...

4.2 PLR: logitprop_model_02

logitprop_model_02 = glm(cbind(Y1, Y0) ~ X, family = binomial(link = logit), data = logitprop)
summary(logitprop_model_02)

Call:
glm(formula = cbind(Y1, Y0) ~ X, family = binomial(link = logit), 
    data = logitprop)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3073  -0.4377   0.1172   0.6303   1.3465  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -5.1093     1.7755  -2.878  0.00401 **
X             0.8406     0.2686   3.130  0.00175 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 33.2289  on 12  degrees of freedom
Residual deviance:  7.0578  on 11  degrees of freedom
AIC: 17.917

Number of Fisher Scoring iterations: 6

antilog_logitprop_model_02 <- exp(coefficients(logitprop_model_02))
antilog_logitprop_model_02

(Intercept)           X 
0.006040245 2.317778111 

antilog_logitbin_model_01 <- exp(coefficients(logitbin_model_01))
antilog_logitbin_model_01

(Intercept)           X 
0.006040245 2.317778111 

4.3 predicted values as a logistic curve: logitprop_model_02

plot(Y1/(Y1 + Y0) ~ X, cex = 1.3, cex.axis = 1.4, cex.lab = 1.4, las = 1, data = logitprop, 
    lty = 1, lwd = 3, col = "blue")
range_x = seq(min(logitprop$X), max(logitprop$X))
lines(range_x, predict(logitprop_model_02, list(X = range_x), type = "resp"), lty = 1, 
    lwd = 2, col = "red")
grid()

4.4 Basic Diagnostic Plot: logitprop_model_02

oldpar <- par(oma = c(0, 0, 3, 0), mfrow = c(2, 2))
plot(logitprop_model_02)

5 Conclusion:

• Logistic Regression of “logitbin_model_01” and proportional logistic regression of “logitprop_model_02” estimate to have the same parameters.
• Both two models shows having poor fit at intermediate values of X.
• There is no indication of overdispersion for both two models.
• Hence conclude that we have a significant positive effect of variable X on the odds of the Y1.
• The odds for Y1 increases 2.317778111 timefolds for each unit increase in X.

Reference

Crawley, Michael J. 2013. “The R Book Second Edition.” John Wiley & Sons.