AS3-2 Predicting Parole Violators

rm(list=ls(all=T))
options(digits=4, scipen=12)
library(magrittr)

---

Introduction

議題：用申請假釋者的屬性，預測他會不會違反假釋規定

學習重點：

• 設定隨機種子set.seed()，依比例分割資料
• 從邏輯式回歸的係數推算自變數的邊際效果
• 勝率和機率、勝率比和機率差 之間的換算
• 臨界機率對混淆矩陣(期望報酬)的影響payoff = matrix(c(0,-100,-10,-50),2,2); payoff
• AUC的實質意義payoff = matrix(c(100, -80, -20, 100),2,2); payoff
• 如何(從報酬矩陣)決定臨界機率
• 什麼是抽樣偏差,如何避免、如何修正

---

1 資料處理 Loading the Dataset

【1.1 讀進資料】How many parolees are contained in the dataset?

P = read.csv("data/parole.csv")
nrow(P)

[1] 675

【1.2 底線機率】How many of the parolees in the dataset violated the terms of their parole?

nrow(subset(P,P$violator == 1))

[1] 78

---

2 整理資料框 Creating Our Prediction Model

【2.1 類別變數】Which variables in this dataset are unordered factors with at least three levels?

#state,crime
P$state %<>% as.factor
P$crime %<>% as.factor

【2.2 資料框摘要】How does the output of summary() change for a factor variable as compared to a numerical variable?

summary(P)

      male            race           age       state  
 Min.   :0.000   Min.   :1.00   Min.   :18.4   1:143  
 1st Qu.:1.000   1st Qu.:1.00   1st Qu.:25.4   2:120  
 Median :1.000   Median :1.00   Median :33.7   3: 82  
 Mean   :0.807   Mean   :1.42   Mean   :34.5   4:330  
 3rd Qu.:1.000   3rd Qu.:2.00   3rd Qu.:42.5          
 Max.   :1.000   Max.   :2.00   Max.   :67.0          
  time.served    max.sentence  multiple.offenses crime  
 Min.   :0.00   Min.   : 1.0   Min.   :0.000     1:315  
 1st Qu.:3.25   1st Qu.:12.0   1st Qu.:0.000     2:106  
 Median :4.40   Median :12.0   Median :1.000     3:153  
 Mean   :4.20   Mean   :13.1   Mean   :0.536     4:101  
 3rd Qu.:5.20   3rd Qu.:15.0   3rd Qu.:1.000            
 Max.   :6.00   Max.   :18.0   Max.   :1.000            
    violator    
 Min.   :0.000  
 1st Qu.:0.000  
 Median :0.000  
 Mean   :0.116  
 3rd Qu.:0.000  
 Max.   :1.000  

#The output becomes similar to that of the table() function applied to that variable correct

---

3 資料分割 Splitting into a Training and Testing Set

To ensure consistent training/testing set splits, run the following 5 lines of code (do not include the line numbers at the beginning):

1. set.seed(144)

2. library(caTools)

3. split = sample.split(parole$violator, SplitRatio = 0.7)

4. train = subset(parole, split == TRUE)

5. test = subset(parole, split == FALSE)

【3.1 指定隨機種子、依比例分割資料】Roughly what proportion of parolees have been allocated to the training and testing sets?

set.seed(144)
library(caTools)
split = sample.split(P$violator,SplitRatio = 0.7)
TR = subset(P,split == TRUE)
TS = subset(P,split == FALSE)

#70% to the training set, 30% to the testing set

【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?

#The exact same training/testing set split as the first execution of [1]-[5]
#A different training/testing set split from the first execution of [1]-[5]
#A different training/testing set split from the first execution of [1]-[5]

---

4 建立模型 Building a Logistic Regression Model

If you tested other training/testing set splits in the previous section, please re-run the original 5 lines of code to obtain the original split.

Using glm (and remembering the parameter family=“binomial”), train a logistic regression model on the training set. Your dependent variable is “violator”, and you should use all of the other variables as independent variables.

What variables are significant in this model? Significant variables should have a least one star, or should have a probability less than 0.05 (the column Pr(>|z|) in the summary output). Select all that apply.

【4.1 顯著性】What variables are significant in this model?

glm1 = glm(violator~.,data=TR,family = binomial)
summary(glm1)

Call:
glm(formula = violator ~ ., family = binomial, data = TR)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.704  -0.424  -0.272  -0.169   2.837  

Coefficients:
                   Estimate Std. Error z value    Pr(>|z|)    
(Intercept)       -4.241157   1.293885   -3.28       0.001 ** 
male               0.386990   0.437961    0.88       0.377    
race               0.886719   0.395066    2.24       0.025 *  
age               -0.000176   0.016085   -0.01       0.991    
state2             0.443301   0.481662    0.92       0.357    
state3             0.834980   0.556270    1.50       0.133    
state4            -3.396788   0.611586   -5.55 0.000000028 ***
time.served       -0.123887   0.120423   -1.03       0.304    
max.sentence       0.080295   0.055375    1.45       0.147    
multiple.offenses  1.611992   0.385305    4.18 0.000028683 ***
crime2             0.683714   0.500355    1.37       0.172    
crime3            -0.278105   0.432836   -0.64       0.521    
crime4            -0.011763   0.571304   -0.02       0.984    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 340.04  on 472  degrees of freedom
Residual deviance: 251.48  on 460  degrees of freedom
AIC: 277.5

Number of Fisher Scoring iterations: 6

#race,state4,multiple.offeses

【4.2 從回歸係數估計邊際效用】What can we say based on the coefficient of the multiple.offenses variable?

The following two properties might be useful to you when answering this question:

1. If we have a coefficient c for a variable, then that means the log odds (or Logit) are increased by c for a unit increase in the variable.

2. If we have a coefficient c for a variable, then that means the odds are multiplied by e^c for a unit increase in the variable.

exp(1.611992)  #exp(multiple.offenses的estimate) 

[1] 5.013

#將logit轉回odd
# 5.013
#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.

Consider a parolee who is male, of white race, aged 50 years at prison release, from the state of Maryland, served 3 months, had a maximum sentence of 12 months, did not commit multiple offenses, and committed a larceny. Answer the following questions based on the model’s predictions for this individual. (HINT: You should use the coefficients of your model, the Logistic Response Function, and the Odds equation to solve this problem.)

【4.3 從預測值估計勝率和機率】According to the model, what are the odds this individual is a violator? What is the probability this individual is a violator?

#log(odds) or logit =   -4.241157 + 0.386990*male + 0.886719* race + (-0.000176)*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
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
logit = -4.241157 + 0.386990*male + 0.886719* race + (-0.000176)*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
odd = exp(logit);odd

[1] 0.1826

p = 1/1+odd ; p

[1] 1.183

---

5 驗證模型 Evaluating the Model on the Testing Set

Use the predict() function to obtain the model’s predicted probabilities for parolees in the testing set, remembering to pass type=“response”.

【5.1 從測試資料預測機率】What is the maximum predicted probability of a violation?

pred = predict(glm1,TS,type="response")
max(pred)

[1] 0.9073

Use the predict() function to obtain the model’s predicted probabilities for parolees in the testing set, remembering to pass type=“response”.

【5.2 從混淆矩陣計算敏感性、明確性、正確率】What is the model’s sensitivity, specificity, accuracy?

A = table(actual = TS$violator , pred = pred >= 0.5);A

      pred
actual FALSE TRUE
     0   167   12
     1    11   12

sen  = 12/(11+12);sen

[1] 0.5217

spec = 167/(167+12);spec

[1] 0.933

acc = sum(diag(A))/sum(A);acc

[1] 0.8861

【5.3 底線機率】What is the accuracy of a simple model that predicts that every parolee is a non-violator?

table(TS$violator) #在test里non-violator有179

  0   1 
179  23 

#如果全部都視為non
ACC = 179/(179+23) ; ACC

[1] 0.8861

【5.4 根據報償矩陣調整臨界機率】 Consider a parole board using the model to predict whether parolees will be violators or not. The job of a parole board is to make sure that a prisoner is ready to be released into free society, and therefore parole boards tend to be particularily concerned about releasing prisoners who will violate their parole

Which of the following most likely describes their preferences and best course of action?

#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?

#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)
preds = prediction(pred, TS$violator)
as.numeric(performance(preds, "auc")@y.values)

[1] 0.8946

【5.7 辨識率的定義】Describe the meaning of AUC in this context.

#The probability the model can correctly differentiate between a randomly selected parole violator and a randomly selected parole non-violator.
#

---

6 抽樣偏差 Identifying Bias in Observational Data

Our goal has been to predict the outcome of a parole decision, and we used a publicly available dataset of parole releases for predictions. In this final problem, we’ll evaluate a potential source of bias associated with our analysis. It is always important to evaluate a dataset for possible sources of bias.

The dataset contains all individuals released from parole in 2004, either due to completing their parole term or violating the terms of their parole. However, it does not contain parolees who neither violated their parole nor completed their term in 2004, causing non-violators to be underrepresented. This is called “selection bias” or “selecting on the dependent variable,” because only a subset of all relevant parolees were included in our analysis, based on our dependent variable in this analysis (parole violation).

【6.1 如何避免、診斷、修正抽樣偏差】How could we improve our dataset to best address selection bias?

#因為TEST資料集只限於某個年份，但忽略了他在之後假釋期間的後續動向，沒有完整記錄
#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.

---