rm(list=ls(all=T))
options(digits=4, scipen=12)
library(magrittr)議題：用申請假釋者的屬性，預測他會不會違反假釋規定
學習重點：
setwd("~/big data/unit3/code_data/data")The working directory was changed to C:/Users/Administrator/Documents/big data/unit3/code_data/data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.parole=read.csv("parole.csv")【1.1 讀進資料】How many parolees are contained in the dataset?
nrow(parole)[1] 675【1.2 底線機率】How many of the parolees in the dataset violated the terms of their parole?
table(parole$violator)[2] 1
78 【2.1 類別變數】Which variables in this dataset are unordered factors with at least three levels?
apply(parole,2,table)$male
0 1
130 545
$race
1 2
389 286
$age
18.4 18.5 18.7 18.8 19 19.1 19.2 19.3 19.4 19.5 19.6 19.7 19.9 20
2 1 2 2 1 3 4 1 4 2 4 2 3 3
20.2 20.3 20.4 20.5 20.6 20.7 20.8 20.9 21 21.1 21.2 21.3 21.4 21.5
4 2 1 4 3 5 2 1 2 2 4 1 4 2
21.6 21.7 21.8 21.9 22 22.1 22.2 22.3 22.4 22.5 22.6 22.8 22.9 23
2 2 2 6 4 2 1 2 3 2 3 3 1 6
23.1 23.2 23.3 23.4 23.6 23.7 23.8 24 24.2 24.3 24.4 24.5 24.6 24.7
2 6 6 3 5 3 2 1 6 2 2 3 1 1
24.8 24.9 25 25.1 25.2 25.3 25.4 25.5 25.6 25.7 25.8 25.9 26 26.3
2 3 3 3 1 4 1 1 7 4 2 3 1 2
26.4 26.5 26.6 26.8 26.9 27 27.1 27.2 27.3 27.4 27.5 27.6 27.7 27.8
2 2 1 3 4 2 2 2 1 1 2 2 2 3
27.9 28 28.1 28.2 28.3 28.4 28.5 28.7 28.8 28.9 29 29.1 29.2 29.5
1 2 3 4 2 3 2 4 5 3 1 2 2 3
29.6 29.7 29.9 30 30.1 30.2 30.3 30.4 30.7 30.8 31 31.1 31.2 31.3
2 3 4 3 2 1 1 1 2 3 6 2 3 1
31.4 31.5 31.6 31.7 31.8 32 32.1 32.2 32.3 32.4 32.5 32.6 32.7 32.8
2 2 1 1 3 2 3 3 1 4 1 1 2 4
32.9 33 33.2 33.3 33.4 33.5 33.6 33.7 33.8 33.9 34 34.1 34.2 34.3
4 2 4 1 1 2 2 4 1 3 2 1 3 2
34.4 34.5 34.6 34.7 34.8 34.9 35 35.1 35.2 35.3 35.4 35.5 35.6 35.8
1 2 1 1 1 3 3 2 1 1 3 1 3 2
35.9 36 36.1 36.2 36.3 36.4 36.5 36.6 36.7 36.8 37 37.2 37.3 37.4
2 1 2 4 4 1 5 1 3 1 1 3 2 2
37.5 37.6 37.8 38 38.1 38.2 38.3 38.4 38.5 38.6 38.7 38.8 38.9 39
3 1 2 3 2 1 2 4 2 3 2 1 2 3
39.1 39.2 39.4 39.5 39.6 39.7 39.8 39.9 40 40.1 40.3 40.4 40.6 40.8
2 4 1 1 2 4 3 1 1 3 5 2 1 1
40.9 41 41.1 41.2 41.3 41.4 41.6 41.7 41.9 42 42.1 42.3 42.4 42.5
2 1 3 3 5 2 1 3 2 3 1 3 3 3
42.6 42.8 43 43.1 43.2 43.3 43.4 43.5 43.6 43.7 43.8 44 44.1 44.2
1 1 3 1 2 2 3 4 4 1 4 3 4 1
44.3 44.4 44.5 44.6 44.7 44.8 44.9 45 45.1 45.4 45.5 45.6 45.8 45.9
2 1 1 1 2 2 3 5 3 4 1 2 1 2
46 46.1 46.2 46.3 46.4 46.5 46.6 46.7 46.8 46.9 47 47.1 47.2 47.3
2 2 1 3 2 1 1 4 1 2 2 3 2 1
47.5 47.7 47.8 48 48.2 48.4 48.5 48.7 48.8 48.9 49 49.3 49.9 50.1
3 1 3 1 3 1 3 3 2 1 2 1 2 2
50.2 50.5 50.6 50.9 51 51.1 51.2 51.3 51.4 51.7 51.8 52.1 52.5 52.6
2 1 1 1 1 4 1 1 1 1 3 1 1 1
53 53.5 53.8 53.9 54.1 54.4 54.5 54.8 54.9 55 55.7 56.4 56.5 56.8
2 3 1 1 1 2 1 1 1 2 1 1 1 2
57.5 58.5 59.4 61.4 61.6 63.4 65.1 67
2 1 2 1 1 1 1 1
$state
1 2 3 4
143 120 82 330
$time.served
0 0.1 0.2 0.3 0.5 0.7 0.8 0.9 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9
1 2 1 1 2 2 3 3 4 1 2 2 1 1 3 3 2
2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 3.2 3.3 3.4 3.5 3.6
6 1 6 1 3 1 1 6 1 11 61 18 19 12 8 6 12
3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5 5.1 5.2 5.3
17 15 12 19 19 28 13 29 27 15 22 16 22 7 30 32 26
5.4 5.5 5.6 5.7 5.8 5.9 6
13 21 16 12 13 18 26
$max.sentence
1 2 3 4 5 6 8 9 10 11 12 13 14 15 16 17 18
4 4 3 3 2 16 20 9 3 5 274 77 66 60 36 15 78
$multiple.offenses
0 1
313 362
$crime
1 2 3 4
315 106 153 101
$violator
0 1
597 78 【2.2 資料框摘要】How does the output of summary() change for a factor variable as compared to a numerical variable?
parole$state = as.factor(parole$state)
parole$crime = as.factor(parole$crime)
summary(parole$state) 1 2 3 4
143 120 82 330 summary(parole$crime) 1 2 3 4
315 106 153 101 c("Ans:The output becomes similar to that of the table() function applied to that variable")[1] "Ans:The output becomes similar to that of the table() function applied to that variable"【3.1 指定隨機種子、依比例分割資料】Roughly what proportion of parolees have been allocated to the training and testing sets?
set.seed(144)
library(caTools)
split = sample.split(parole$violator, SplitRatio = 0.7)
train = subset(parole, split == TRUE)
test = subset(parole, split == FALSE)
nrow(train)/nrow(parole)[1] 0.7007【3.2 隨機種子的功用】Now, suppose you re-ran lines [1]-[5] of Problem 3.1. What would you expect? If you instead ONLY re-ran lines [3]-[5], what would you expect? If you instead called set.seed() with a different number and then re-ran lines [3]-[5] of Problem 3.1, what would you expect?
c("Ans:The exact same training/testing set split as the first execution of [1]-[5]")[1] "Ans:The exact same training/testing set split as the first execution of [1]-[5]"c("Ans:A different training/testing set split from the first execution of [1]-[5]")[1] "Ans:A different training/testing set split from the first execution of [1]-[5]"c("Ans:A different training/testing set split from the first execution of [1]-[5]")[1] "Ans:A different training/testing set split from the first execution of [1]-[5]"【4.1 顯著性】What variables are significant in this model?
mod = glm(violator~., data=train, family="binomial")
names(which(summary(mod)$coefficients[,4]<0.05))[-1][1] "race" "state4" "multiple.offenses"【4.2 從回歸係數估計邊際效用】What can we say based on the coefficient of the multiple.offenses variable?
exp(summary(mod)$coefficients[10,1])[1] 5.013c("Ans:Our model predicts that a parolee who committed multiple offenses has 5.01 times higher odds of being a violator than a parolee who did not commit multiple offenses but is otherwise identical. ")[1] "Ans:Our model predicts that a parolee who committed multiple offenses has 5.01 times higher odds of being a violator than a parolee who did not commit multiple offenses but is otherwise identical. "【4.3 從預測值估計勝率和機率】According to the model, what are the odds this individual is a violator? What is the probability this individual is a violator?
male=1;race=1;age=50;state2=0;state3=0;state4=0;time.served=3; max.sentence=12;multiple.offenses=0;crime2=1;crime3=0;crime4=0
odds = -4.2411574 + 0.3869904*male + 0.8867192*race - 0.0001756*age + 0.4433007*state2 + 0.8349797*state3 - 3.3967878*state4 - 0.1238867*time.served + 0.0802954*max.sentence + 1.6119919*multiple.offenses + 0.6837143*crime2 - 0.2781054*crime3 - 0.0117627*crime4
exp(odds)[1] 0.18261/(1+exp(-odds)) [1] 0.1544【5.1 從測試資料預測機率】What is the maximum predicted probability of a violation?
predictions = predict(mod, newdata=test, type="response")
max(summary(predictions))[1] 0.9073【5.2 從混淆矩陣計算敏感性、明確性、正確率】What is the model’s sensitivity, specificity, accuracy?
table(test$violator, as.numeric(predictions >= 0.5))
0 1
0 167 12
1 11 12sensitivity=12/(11+12);sensitivity [1] 0.5217specificity=167/(167+12);specificity[1] 0.933accuracy=(167+12)/202;accuracy[1] 0.8861【5.3 底線機率】What is the accuracy of a simple model that predicts that every parolee is a non-violator?
table(test$violator)
0 1
179 23 accuracy=179/202;accuracy[1] 0.8861【5.4 根據報償矩陣調整臨界機率】Which of the following most likely describes their preferences and best course of action?
c("Ans:The board assigns more cost to a false negative than a false positive, and should therefore use a logistic regression cutoff less than 0.5")[1] "Ans:The board assigns more cost to a false negative than a false positive, and should therefore use a logistic regression cutoff less than 0.5"【5.5 正確率 vs 辨識率】Which of the following is the most accurate assessment of the value of the logistic regression model with a cutoff 0.5 to a parole board, based on the model’s accuracy as compared to the simple baseline model?
c("Ans:The model is likely of value to the board, and using a different logistic regression cutoff is likely to improve the model's value")[1] "Ans:The model is likely of value to the board, and using a different logistic regression cutoff is likely to improve the model's value"【5.6 計算辨識率】Using the ROCR package, what is the AUC value for the model?
library(ROCR)
pred = prediction(predictions, test$violator)
as.numeric(performance(pred, "auc")@y.values)[1] 0.8946【5.7 辨識率的定義】Describe the meaning of AUC in this context.
c("Ans:The probability the model can correctly differentiate between a randomly selected parole violator and a randomly selected parole non-violator")[1] "Ans:The probability the model can correctly differentiate between a randomly selected parole violator and a randomly selected parole non-violator"【6.1 如何避免、診斷、修正抽樣偏差】How could we improve our dataset to best address selection bias?
c("Ans:We should use a dataset tracking a group of parolees from the start of their parole until either they violated parole or they completed their term")[1] "Ans:We should use a dataset tracking a group of parolees from the start of their parole until either they violated parole or they completed their term"