Blog3: Logistic Regression

## load library
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("gridExtra")) install.packages("gridExtra")

Introduction

In the situation the response variable is based “yes”/“no” responses, such as whether a particular restaurant is recommended by being included in a prestigious guide. Ideally such responses follow a binomial distribution in which case the appropriate model is a logistic regression model.

Dataset

MichelinFood <- read.table("https://gattonweb.uky.edu/sheather/book/docs/datasets/MichelinFood.txt", header=TRUE)

attach(MichelinFood)
plot(Food,proportion,ylab="Sample proportion",xlab="Zagat Food Rating")

glm() is the function that tells R to run a generalized linear model.

Inside the parentheses we give R important information about the model. To the left of the ~ is the dependent variable: success. It must be coded 0 & 1 for glm to read it as binary.

#Logistic Model
m1 <- glm(cbind(InMichelin,NotInMichelin)~Food,family=binomial)
summary(m1)

## 
## Call:
## glm(formula = cbind(InMichelin, NotInMichelin) ~ Food, family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4850  -0.7987  -0.1679   0.5913   1.5889  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -10.84154    1.86236  -5.821 5.84e-09 ***
## Food          0.50124    0.08768   5.717 1.08e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 61.427  on 13  degrees of freedom
## Residual deviance: 11.368  on 12  degrees of freedom
## AIC: 41.491
## 
## Number of Fisher Scoring iterations: 4

The estimates (coefficients of the predictors – numeracy and anxiety) are now in logits. The coefficient of numeracy is: 0.50124, so that a one unit change in numeracy produces approximately a 0.50124 unit change in the log odds.

For our example, we have a Null Deviance of about 61.427 on 13 degrees of freedom. This value indicates poor fit (a significant difference between fitted values and observed values). Including the independent variables (numeracy and anxiety) decreased the deviance by nearly 40 points on 1 degrees of freedom. The Residual Deviance is 11.368 on 12 degrees of freedom

x <- seq(15,28,0.05)
y <- 1/(1+exp(-1*(m1$coeff[1] + m1$coeff[2]*x)))
plot(Food,proportion,ylab="Probability of inclusion in the Michelin Guide",xlab="Zagat Food Rating")
lines(x,y)

Reference

Simon J. Sheather, “A Modern Approach to Regression with R”