Logisitic regression¶

Xiaoqi Zheng, 0328/2020¶

1. Binary classification¶

# Select some columns form mtcars.
input <- mtcars[,c("am","cyl","hp","wt")]
print(head(input))

                  am cyl  hp    wt
Mazda RX4          1   6 110 2.620
Mazda RX4 Wag      1   6 110 2.875
Datsun 710         1   4  93 2.320
Hornet 4 Drive     0   6 110 3.215
Hornet Sportabout  0   8 175 3.440
Valiant            0   6 105 3.460

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.0000  0.0000  0.4062  1.0000  1.0000

lr_model = glm(formula = am ~ cyl + hp + wt, data = input, family = binomial) #binomial

print(summary(lr_model))

Call:
glm(formula = am ~ cyl + hp + wt, family = binomial, data = input)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.17272  -0.14907  -0.01464   0.14116   1.27641  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) 19.70288    8.11637   2.428   0.0152 *
cyl          0.48760    1.07162   0.455   0.6491  
hp           0.03259    0.01886   1.728   0.0840 .
wt          -9.14947    4.15332  -2.203   0.0276 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 43.2297  on 31  degrees of freedom
Residual deviance:  9.8415  on 28  degrees of freedom
AIC: 17.841

Number of Fisher Scoring iterations: 8

test_prob = predict(object = lr_model,newdata = input,type = 'response')
test_prob

glm.pred = ifelse(test_prob > 0.5, "1", "0")
table(glm.pred,input$am)

        
glm.pred  0  1
       0 18  2
       1  1 11

## training accuracy
mean(glm.pred == input$am)

## plot ROC curve
library(ROCR)
p <- predict(lr_model, newdata=input, type="response")
pr <- prediction(p, input$am)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
plot(prf,main = "ROC curve",lwd = 2,col = "red")

library(pROC)
auc = auc(response = input$am,predictor = p)
text(0.8,0.2,labels = paste0("AUC = ",round(auc,4)))

Setting levels: control = 0, case = 1

Setting direction: controls < cases

2. Multi-group classification¶

Multinomial logistic regression can be implemented with mlogit() from mlogit package and multinom() from nnet package. We will use the latter for this example.

cmcData <- read.csv("http://archive.ics.uci.edu/ml/machine-learning-databases/cmc/cmc.data", stringsAsFactors=FALSE, header=F)
colnames(cmcData) <- c("wife_age", "wife_edu", "hus_edu", "num_child", "wife_rel", "wife_work", "hus_occu", "sil", "media_exp", "cmc")
head(cmcData)

# Prepare Training and Test Data
set.seed(100)
trainingRows <- sample(1:nrow(cmcData), 0.7*nrow(cmcData))
training <- cmcData[trainingRows, ]
test <- cmcData[-trainingRows, ]

library(nnet)
multinomModel <- multinom(cmc ~ ., data=training) # multinom Model
summary (multinomModel) # model summary

# weights:  33 (20 variable)
initial  value 1132.669270 
iter  10 value 1005.819763
iter  20 value 975.922721
final  value 974.824494 
converged

Call:
multinom(formula = cmc ~ ., data = training)

Coefficients:
  (Intercept)   wife_age  wife_edu     hus_edu num_child   wife_rel   wife_work
2  -3.2863389 -0.0475747 0.8641614  0.03270041 0.3143566 -0.4876721 -0.01258458
3   0.7014815 -0.1030796 0.3194865 -0.01320350 0.2836229 -0.2368861  0.13903133
    hus_occu       sil  media_exp
2 -0.1210229 0.3289317 -0.1245895
3  0.1346765 0.1527346 -0.6424231

Std. Errors:
  (Intercept)   wife_age  wife_edu   hus_edu  num_child  wife_rel wife_work
2   0.9333410 0.01436321 0.1335053 0.1637365 0.05029112 0.2375021 0.2005847
3   0.7513528 0.01318821 0.1018969 0.1199011 0.04445760 0.2408192 0.1795966
   hus_occu        sil media_exp
2 0.1156028 0.11741720  0.426903
3 0.1005491 0.08660681  0.316616

Residual Deviance: 1949.649 
AIC: 1989.649

predicted_scores <- predict (multinomModel, test, "probs") # predict on new data
predicted_class <- predict (multinomModel, test)

table(predicted_class, test$cmc)

               
predicted_class   1   2   3
              1 118  28  46
              2  13  38  35
              3  61  30  73

# model accuracy
mean(predicted_class == test$cmc)

	wife_age	wife_edu	hus_edu	num_child	wife_rel	wife_work	hus_occu	sil	media_exp	cmc
	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>
1	24	2	3	3	1	1	2	3	0	1
2	45	1	3	10	1	1	3	4	0	1
3	43	2	3	7	1	1	3	4	0	1
4	42	3	2	9	1	1	3	3	0	1
5	36	3	3	8	1	1	3	2	0	1
6	19	4	4	0	1	1	3	3	0	1