Example 1: Study Hours

df <- read.table("studytime.txt", sep="\t", header=TRUE)
summary(df)
     Passed         Hours      
 Min.   :0.00   Min.   : 50.0  
 1st Qu.:0.00   1st Qu.:116.2  
 Median :1.00   Median :189.0  
 Mean   :0.56   Mean   :178.2  
 3rd Qu.:1.00   3rd Qu.:241.5  
 Max.   :1.00   Max.   :299.0  

Create Plots

plot(Passed ~ Hours, df)

We will covert Passed to a factor variable, and will plot the data again.

df$Passed <- factor(df$Passed)
summary(df)
 Passed     Hours      
 0:44   Min.   : 50.0  
 1:56   1st Qu.:116.2  
        Median :189.0  
        Mean   :178.2  
        3rd Qu.:241.5  
        Max.   :299.0  
plot(Passed ~ Hours, df)

plot(Hours ~ Passed, df)

Create the Logistic Regression Model

# Create Logistic Regression Model 
m1 <- glm(Passed ~ Hours, data=df, family=binomial(link="logit"))

Generate Predictions

# Create new data consisting of vector of study times
newdata = data.frame(Hours = c(50, 100, 150, 200, 250, 300))
# Create predictions
logodds <- predict(m1, newdata)     # This produces a prediction of log(odds)
logodds
         1          2          3          4          5          6 
-2.2020131 -1.2189354 -0.2358576  0.7472202  1.7302980  2.7133758 
odds <- exp(logodds)                # Prediction of odds
odds
         1          2          3          4          5          6 
 0.1105803  0.2955447  0.7898932  2.1111234  5.6423351 15.0800968 
odds / (1 + odds)             # Prediction of probabilities
         1          2          3          4          5          6 
0.09956985 0.22812386 0.44130744 0.67857270 0.84945053 0.93781132 
predict(m1, newdata, type="response")
         1          2          3          4          5          6 
0.09956985 0.22812386 0.44130744 0.67857270 0.84945053 0.93781132 

Example 2: Graduate School Admissions

df <- read.table("admissions.txt", sep="\t", header=TRUE)
summary(df)
     admit             gre             gpa             rank      
 Min.   :0.0000   Min.   :220.0   Min.   :2.260   Min.   :1.000  
 1st Qu.:0.0000   1st Qu.:520.0   1st Qu.:3.130   1st Qu.:2.000  
 Median :0.0000   Median :580.0   Median :3.395   Median :2.000  
 Mean   :0.3175   Mean   :587.7   Mean   :3.390   Mean   :2.485  
 3rd Qu.:1.0000   3rd Qu.:660.0   3rd Qu.:3.670   3rd Qu.:3.000  
 Max.   :1.0000   Max.   :800.0   Max.   :4.000   Max.   :4.000  

Create Plots

df$admit = factor(df$admit)
df$rank = factor(df$rank)
summary(df)
 admit        gre             gpa        rank   
 0:273   Min.   :220.0   Min.   :2.260   1: 61  
 1:127   1st Qu.:520.0   1st Qu.:3.130   2:151  
         Median :580.0   Median :3.395   3:121  
         Mean   :587.7   Mean   :3.390   4: 67  
         3rd Qu.:660.0   3rd Qu.:3.670          
         Max.   :800.0   Max.   :4.000          
plot(gre ~ admit, df)

plot(gpa ~ admit, df)

plot(rank ~ admit, df)

Create the Logistic Regression Model

# Create Logistic Regression Model 
m2 <- glm(admit ~ gre + gpa + rank, data=df, family=binomial(link="logit"))
summary(m2)

Call:
glm(formula = admit ~ gre + gpa + rank, family = binomial(link = "logit"), 
    data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6268  -0.8662  -0.6388   1.1490   2.0790  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.989979   1.139951  -3.500 0.000465 ***
gre          0.002264   0.001094   2.070 0.038465 *  
gpa          0.804038   0.331819   2.423 0.015388 *  
rank2       -0.675443   0.316490  -2.134 0.032829 *  
rank3       -1.340204   0.345306  -3.881 0.000104 ***
rank4       -1.551464   0.417832  -3.713 0.000205 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 458.52  on 394  degrees of freedom
AIC: 470.52

Number of Fisher Scoring iterations: 4

Generate Predictions

nd = data.frame( gre = c(720, 720, 720, 720, 550, 550, 550, 550),
                 gpa = c(3.2, 3.2, 3.2, 3.2, 3.8, 3.8, 3.8, 3.8),
                 rank = c('1', '2', '3', '4', '1', '2', '3', '4') )
predict(m2, nd, type="response")
        1         2         3         4         5         6         7         8 
0.5531306 0.3864841 0.2447380 0.2078168 0.5770800 0.4098356 0.2631993 0.2243201 
logodds <- predict(m2, nd)
odds <- exp(logodds)
prob <- odds / (1 + odds)
prob
        1         2         3         4         5         6         7         8 
0.5531306 0.3864841 0.2447380 0.2078168 0.5770800 0.4098356 0.2631993 0.2243201 
LS0tDQp0aXRsZTogIkxlc3NvbiAyNSAtIExvZ2lzdGljIFJlZ3Jlc3Npb24iDQphdXRob3I6ICJSb2JiaWUgQmVhbmUiDQpvdXRwdXQ6DQogIGh0bWxfbm90ZWJvb2s6DQogICAgdGhlbWU6IGZsYXRseQ0KICAgIHRvYzogeWVzDQogICAgdG9jX2RlcHRoOiAyDQotLS0NCg0KIyBFeGFtcGxlIDE6IFN0dWR5IEhvdXJzDQoNCg0KYGBge3J9DQpkZiA8LSByZWFkLnRhYmxlKCJzdHVkeXRpbWUudHh0Iiwgc2VwPSJcdCIsIGhlYWRlcj1UUlVFKQ0Kc3VtbWFyeShkZikNCmBgYA0KDQojIyBDcmVhdGUgUGxvdHMNCg0KYGBge3J9DQpwbG90KFBhc3NlZCB+IEhvdXJzLCBkZikNCmBgYA0KDQpXZSB3aWxsIGNvdmVydCBgUGFzc2VkYCB0byBhIGZhY3RvciB2YXJpYWJsZSwgYW5kIHdpbGwgcGxvdCB0aGUgZGF0YSBhZ2Fpbi4gDQoNCmBgYHtyfQ0KZGYkUGFzc2VkIDwtIGZhY3RvcihkZiRQYXNzZWQpDQpzdW1tYXJ5KGRmKQ0KYGBgDQoNCmBgYHtyfQ0KcGxvdChQYXNzZWQgfiBIb3VycywgZGYpDQpgYGANCg0KDQpgYGB7cn0NCnBsb3QoSG91cnMgfiBQYXNzZWQsIGRmKQ0KYGBgDQoNCiMjIENyZWF0ZSB0aGUgTG9naXN0aWMgUmVncmVzc2lvbiBNb2RlbCANCg0KYGBge3J9DQojIENyZWF0ZSBMb2dpc3RpYyBSZWdyZXNzaW9uIE1vZGVsIA0KbTEgPC0gZ2xtKFBhc3NlZCB+IEhvdXJzLCBkYXRhPWRmLCBmYW1pbHk9Ymlub21pYWwobGluaz0ibG9naXQiKSkNCnN1bW1hcnkobTEpDQpgYGANCg0KIyMgR2VuZXJhdGUgUHJlZGljdGlvbnMNCg0KYGBge3J9DQojIENyZWF0ZSBuZXcgZGF0YSBjb25zaXN0aW5nIG9mIHZlY3RvciBvZiBzdHVkeSB0aW1lcw0KbmV3ZGF0YSA9IGRhdGEuZnJhbWUoSG91cnMgPSBjKDUwLCAxMDAsIDE1MCwgMjAwLCAyNTAsIDMwMCkpDQoNCiMgQ3JlYXRlIHByZWRpY3Rpb25zDQpsb2dvZGRzIDwtIHByZWRpY3QobTEsIG5ld2RhdGEpICAgICAjIFRoaXMgcHJvZHVjZXMgYSBwcmVkaWN0aW9uIG9mIGxvZyhvZGRzKQ0KbG9nb2Rkcw0KYGBgDQoNCmBgYHtyfQ0Kb2RkcyA8LSBleHAobG9nb2RkcykgICAgICAgICAgICAgICAgIyBQcmVkaWN0aW9uIG9mIG9kZHMNCm9kZHMNCmBgYA0KDQpgYGB7cn0NCm9kZHMgLyAoMSArIG9kZHMpICAgICAgICAgICAgICMgUHJlZGljdGlvbiBvZiBwcm9iYWJpbGl0aWVzDQpgYGANCg0KDQpgYGB7cn0NCnByZWRpY3QobTEsIG5ld2RhdGEsIHR5cGU9InJlc3BvbnNlIikNCmBgYA0KDQojIEV4YW1wbGUgMjogR3JhZHVhdGUgU2Nob29sIEFkbWlzc2lvbnMNCg0KDQpgYGB7cn0NCmRmIDwtIHJlYWQudGFibGUoImFkbWlzc2lvbnMudHh0Iiwgc2VwPSJcdCIsIGhlYWRlcj1UUlVFKQ0Kc3VtbWFyeShkZikNCmBgYA0KDQojIyBDcmVhdGUgUGxvdHMNCg0KYGBge3J9DQpkZiRhZG1pdCA9IGZhY3RvcihkZiRhZG1pdCkNCmRmJHJhbmsgPSBmYWN0b3IoZGYkcmFuaykNCnN1bW1hcnkoZGYpDQpgYGANCg0KDQpgYGB7cn0NCnBsb3QoZ3JlIH4gYWRtaXQsIGRmKQ0KYGBgDQoNCg0KYGBge3J9DQpwbG90KGdwYSB+IGFkbWl0LCBkZikNCmBgYA0KDQoNCmBgYHtyfQ0KcGxvdChyYW5rIH4gYWRtaXQsIGRmKQ0KYGBgDQoNCg0KIyMgQ3JlYXRlIHRoZSBMb2dpc3RpYyBSZWdyZXNzaW9uIE1vZGVsIA0KDQpgYGB7cn0NCiMgQ3JlYXRlIExvZ2lzdGljIFJlZ3Jlc3Npb24gTW9kZWwgDQptMiA8LSBnbG0oYWRtaXQgfiBncmUgKyBncGEgKyByYW5rLCBkYXRhPWRmLCBmYW1pbHk9Ymlub21pYWwobGluaz0ibG9naXQiKSkNCnN1bW1hcnkobTIpDQpgYGANCg0KIyMgR2VuZXJhdGUgUHJlZGljdGlvbnMNCg0KYGBge3J9DQpuZCA9IGRhdGEuZnJhbWUoIGdyZSA9IGMoNzIwLCA3MjAsIDcyMCwgNzIwLCA1NTAsIDU1MCwgNTUwLCA1NTApLA0KICAgICAgICAgICAgICAgICBncGEgPSBjKDMuMiwgMy4yLCAzLjIsIDMuMiwgMy44LCAzLjgsIDMuOCwgMy44KSwNCiAgICAgICAgICAgICAgICAgcmFuayA9IGMoJzEnLCAnMicsICczJywgJzQnLCAnMScsICcyJywgJzMnLCAnNCcpICkNCg0KcHJlZGljdChtMiwgbmQsIHR5cGU9InJlc3BvbnNlIikNCmBgYA0KDQoNCmBgYHtyfQ0KbG9nb2RkcyA8LSBwcmVkaWN0KG0yLCBuZCkNCm9kZHMgPC0gZXhwKGxvZ29kZHMpDQpwcm9iIDwtIG9kZHMgLyAoMSArIG9kZHMpDQoNCnByb2INCmBgYA0KDQoNCg0KDQo=