AS3-0

rm(list=ls(all=T))
options(digits=4, scipen=12)
library(dplyr); library(ggplot2)

Attaching package: 愼㸱愼㸵dplyr愼㸱愼㸶

The following objects are masked from 愼㸱愼㸵package:stats愼㸱愼㸶:

    filter, lag

The following objects are masked from 愼㸱愼㸵package:base愼㸱愼㸶:

    intersect, setdiff, setequal, union

Use suppressPackageStartupMessages() to eliminate package startup messages.

---

【A】 Definitions

機率、勝率(Odd)、Logit

• Odd = \(p/(1-p)\)

• Logit = \(log(odd)\) = \(log(\frac{p}{1=p})\)

• \(o = p/(1-p)\) ; \(p = o/(1+o)\) ; \(logit = log(o)\)

par(cex=0.8, mfcol=c(1,2))
curve(x/(1-x), 0.02, 0.98, col='cyan',lwd=2, main='odd')
abline(v=seq(0,1,0.1), h=seq(0,50,5), col='lightgray', lty=3)

curve(log(x/(1-x)), 0.005, 0.995, lwd=2, col='purple', main="logit")
abline(v=seq(0,1,0.1), h=seq(-5,5,1), col='lightgray', lty=3)

Logistic Function & Logistic Regression

• Linear Model: \(y = f(x) = b_0 + b_1x_1 + b_2x_2 + ...\)

• General Linear Model(GLM): \(y = Link(f(x))\)

• Logistic Regression: \(logit(y) = log(\frac{p}{1-p}) = f(x) \text{ where } p = prob[y=1]\)

• Logistic Function: \(Logistic(F_x) = \frac{1}{1+Exp(-F_x)} = \frac{Exp(F_x)}{1+Exp(F_x)}\)

par(cex=0.8)
curve(1/(1+exp(-x)), -5, 5, col='blue', lwd=2,main="Logistic Function",
      xlab="f(x): the logit of y = 1", ylab="the probability of y = 1")
abline(v=-5:5, h=seq(0,1,0.1), col='lightgray', lty=2)

abline(v=0,h=0.5,col='pink')
points(0,0.5,pch=20,cex=1.5,col='red')

【Q】What are the definiion of ligit & logistic function? What is the relationship between them?

logit : 機率的一種算法 => log(odd) logistic function : 把logit轉成y=1的機率。

---

【B】glm(, family=binomial)

glm()的功能：在 \(\{x\}\) 的空間之中，找出區隔 \(y\) 的(類別)界線

Q = read.csv("data/quality.csv")
glm1 = glm(PoorCare~OfficeVisits+Narcotics, Q, family=binomial)
summary(glm1)

Call:
glm(formula = PoorCare ~ OfficeVisits + Narcotics, family = binomial, 
    data = Q)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.377  -0.627  -0.510  -0.154   2.119  

Coefficients:
             Estimate Std. Error z value    Pr(>|z|)    
(Intercept)   -2.5402     0.4500   -5.64 0.000000017 ***
OfficeVisits   0.0627     0.0240    2.62     0.00892 ** 
Narcotics      0.1099     0.0326    3.37     0.00076 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 147.88  on 130  degrees of freedom
Residual deviance: 116.45  on 128  degrees of freedom
AIC: 122.4

Number of Fisher Scoring iterations: 5

b = coef(glm1); b   # extract the regression coef

 (Intercept) OfficeVisits    Narcotics 
    -2.54021      0.06273      0.10990 

Given OfficeVisits=3, Narcotics=4, what are the predicted logit, odd and probability?

logit = sum(b * c(1, 3, 4))
odd = exp(logit)
prob = odd/(1+odd)
c(logit=logit, odd=odd, prob=prob)

  logit     odd    prob 
-1.9124  0.1477  0.1287 

【Q】What if OfficeVisits=2, Narcotics=3?

logit_q = sum(b * c(1, 2, 3))
odd_q = exp(logit_q)
prob_q = odd_q / (1+odd_q) # => (p / 1-p) / (1/(1-p)) = p*(1-p) / 1(1-p) = p
prob_q2 = 1 / (1+exp(-1*logit_q)) #Another way to calculate probabilitys from logit.
c(logit=logit_q, odd=odd_q, prob=prob_q)

  logit     odd    prob 
-2.0851  0.1243  0.1106 

#

We can plot the line of logit = 0 or prob = 0.5 on the plane of \(X\)

par(cex=0.8)
plot(Q$OfficeVisits, Q$Narcotics, col=1+Q$PoorCare,pch=20)
abline(-b[1]/b[3], -b[2]/b[3])

Furthermore, we can translate probability, logit and coefficents to intercept & slope …

\[f(x) = b_1 + b_2 x_2 + b_3 x_3 = g \Rightarrow x_3 = \frac{g - b_1}{b_3} - \frac{b_2}{b_3}x_2\]

p = seq(0.1,0.9,0.1)
logit = log(p/(1-p))
data.frame(prob = p, logit)

then mark the contours of proabilities into the scatter plot

par(cex=0.7)
plot(Q$OfficeVisits, Q$Narcotics, col=1+Q$PoorCare,
     pch=20, cex=1.3, xlab='X2', ylab='X3')
for(g in logit) {
  abline((g-b[1])/b[3], -b[2]/b[3], col=ifelse(g==0,'blue','cyan')) }

【Q】What do the blue/cyan lines means? 機率從0.1~0.9的線 【Q】Given any point in the figure above, how can you tell its (predicted) probability approximately? 因為接近0.7的線，所以該點機率大約為0.7

---

【C】The Confusion Matrix

Figure 1 - Confusion Matrix

---

【D】The Distribution of Predicted Probability (DPP)

Confusion matrix is not fixed. It changes by Threshold …

Figure 2 - Dist. Prediected Prob.

library(caTools)
DPP2 = function(pred,class,tvalue,breaks=0.01) {
  mx = table(class == tvalue, pred > 0.5) 
  tn = sum(class != tvalue & pred <= 0.5)
  fn = sum(class == tvalue & pred <= 0.5)
  fp = sum(class != tvalue & pred > 0.5)
  tp = sum(class == tvalue & pred > 0.5)
  acc = (tn + tp)/length(pred)
  sens = tp/(fn+tp)
  spec = tn/(tn+fp)
  auc = colAUC(pred,class)
  data.frame(pred,class) %>% 
    ggplot(aes(x=pred, fill=class)) +
    geom_histogram(col='gray',alpha=0.5,breaks=seq(0,1,breaks)) +
    xlim(0,1) + theme_bw() + xlab("predicted probability") + 
    ggtitle(
      sprintf("Distribution of Prob[class = \'%s\']", tvalue),
      sprintf("AUC=%.3f, Acc=%.3f, Sens=%.3f, Spec=%.3f",
              auc, acc, sens, spec) ) 
  }

N1 = 300; N2 = 100
DPP2(pred = c(rnorm(N1,0.125,0.03), rnorm(N2,0.375,0.03)),
     class = c(rep('B',N1), rep('A',N2)), 
     tvalue = 'A')

【Q】Is it possible to have AUC = ACC = SENS = SPEC = 1? Can you modify the code to make it happen?

N1 = 300; N2 = 100
DPP2(pred = c(rnorm(N1,0.125,0.03), rnorm(N2,0.625,0.03)),
     class = c(rep('B',N1), rep('A',N2)), 
     tvalue = 'A')

【Q】Is it possible to have AUC = ACC = SENS = SPEC = 0? Can you modify the code to make that happen?

N1 = 300; N2 = 100
DPP2(pred = c(rnorm(N1,0.875,0.03), rnorm(N2,0.125,0.03)),
     class = c(rep('B',N1), rep('A',N2)), 
     tvalue = 'A')

---

【E】Modeling Expert

E1: Random Split

set.seed(88)
split = sample.split(Q$PoorCare, SplitRatio = 0.75)
table(split) %>% prop.table()

split
 FALSE   TRUE 
0.2443 0.7557 

table(y = Q$PoorCare, split) %>% prop.table(2)

   split
y    FALSE   TRUE
  0 0.7500 0.7475
  1 0.2500 0.2525

TR = subset(Q, split == TRUE)
TS = subset(Q, split == FALSE)

E2: Build Model

glm1 = glm(PoorCare ~ OfficeVisits + Narcotics, TR, family=binomial)
summary(glm1)

Call:
glm(formula = PoorCare ~ OfficeVisits + Narcotics, family = binomial, 
    data = TR)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0630  -0.6316  -0.5050  -0.0969   2.1669  

Coefficients:
             Estimate Std. Error z value   Pr(>|z|)    
(Intercept)   -2.6461     0.5236   -5.05 0.00000043 ***
OfficeVisits   0.0821     0.0305    2.69     0.0072 ** 
Narcotics      0.0763     0.0321    2.38     0.0173 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 111.888  on 98  degrees of freedom
Residual deviance:  89.127  on 96  degrees of freedom
AIC: 95.13

Number of Fisher Scoring iterations: 4

E3: Prediction & Evaluation

pred = predict(glm1, type='response')
mx = table(TR$PoorCare, pred > 0.5); mx

   
    FALSE TRUE
  0    70    4
  1    15   10

c(accuracy = sum(diag(mx))/sum(mx),
  sensitivity = mx[2,2]/sum(mx[2,]),
  specificity = mx[1,1]/sum(mx[1,]))

   accuracy sensitivity specificity 
     0.8081      0.4000      0.9459 

E4: ROC & AUC

library(ROCR)
ROCRpred = prediction(pred, TR$PoorCare)
ROCRperf = performance(ROCRpred, "tpr", "fpr")
par(cex=0.8)
plot(ROCRperf, colorize=TRUE, print.cutoffs.at=seq(0,1,0.1))

as.numeric(performance(ROCRpred, "auc")@y.values)

[1] 0.7746

caTools::colAUC(pred, TR$PoorCare)

          [,1]
0 vs. 1 0.7746

---

【F】Framingham Heart Study

source("DPP.R")

F1: Reading & Splitting

F = read.csv("data/framingham.csv")
set.seed(1000)
split = sample.split(F$TenYearCHD, SplitRatio = 0.65)
TR = subset(F, split==TRUE)
TS = subset(F, split==FALSE)

F2: Logistic Regression Model

glm2 = glm(TenYearCHD ~ ., TR, family=binomial)
summary(glm2)

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

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.849  -0.601  -0.426  -0.284   2.837  

Coefficients:
                Estimate Std. Error z value        Pr(>|z|)    
(Intercept)     -7.88657    0.89073   -8.85         < 2e-16 ***
male             0.52846    0.13544    3.90 0.0000955212349 ***
age              0.06206    0.00834    7.44 0.0000000000001 ***
education       -0.05892    0.06243   -0.94          0.3453    
currentSmoker    0.09324    0.19401    0.48          0.6308    
cigsPerDay       0.01501    0.00783    1.92          0.0551 .  
BPMeds           0.31122    0.28741    1.08          0.2789    
prevalentStroke  1.16579    0.57121    2.04          0.0413 *  
prevalentHyp     0.31582    0.17176    1.84          0.0660 .  
diabetes        -0.42149    0.40799   -1.03          0.3016    
totChol          0.00384    0.00138    2.79          0.0053 ** 
sysBP            0.01134    0.00457    2.48          0.0130 *  
diaBP           -0.00474    0.00800   -0.59          0.5535    
BMI              0.01072    0.01616    0.66          0.5069    
heartRate       -0.00810    0.00531   -1.52          0.1274    
glucose          0.00893    0.00284    3.15          0.0016 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2020.7  on 2384  degrees of freedom
Residual deviance: 1792.3  on 2369  degrees of freedom
  (371 observations deleted due to missingness)
AIC: 1824

Number of Fisher Scoring iterations: 5

F3: Prediction & Evaluation

pred = predict(glm2, TS, type="response")
y = TS$TenYearCHD[!is.na(pred)]             # remove NA
pred = pred[!is.na(pred)]
mx = table(y, pred > 0.5); mx

   
y   FALSE TRUE
  0  1069    6
  1   187   11

c(accuracy = sum(diag(mx))/sum(mx),
  sensitivity = mx[2,2]/sum(mx[2,]),
  specificity = mx[1,1]/sum(mx[1,]))

   accuracy sensitivity specificity 
    0.84839     0.05556     0.99442 

F4: AUC & DPP

par(cex=0.7)
auc = DPP(pred, y, 1, b=seq(0,1,0.02))  # 0.74211

F5: Expected Result & Optimization

Figure 3 - Startegic Optimization

payoff = matrix(c(0,-100,-10,-60),2,2) 
cutoff = seq(0.02, 0.7, 0.01)
result = sapply(cutoff, function(p) sum(table(y,pred>p)*payoff) )
i = which.max(result)
par(cex=0.7)
plot(cutoff, result, type='l', col='cyan', lwd=2, main=sprintf(
  "Optomal Expected Result: $%d @ %.2f",result[i],cutoff[i]))
abline(v=seq(0,1,0.05),h=seq(-23000,-17000,500),col='lightgray',lty=3)

abline(v=cutoff[i],col='red')

【Q】如果什麼都不做，期望報酬是多少？ 趨近於-20500元
【Q】如果每位病人都做呢？ -22580元
【Q】以上哪一種做法期望報酬比較高？ 什麼都不做
【Q】在所有的商務情境都是這種狀況嗎？ 當誤判的成本極高時，常會出現這樣的狀況。
【Q】你可以模擬出「全做」比「全不做」還要好的狀況、並舉出一個會發生這種狀況的商務情境嗎？ 如果做錯了不會有任何損失，但做對了會有很棒的利潤，則可以去做。
當你不做卻發生意外(False Negative)的時候產生的成本遠高於做時候，就可以去做。（飛安定期檢查)

payoff = matrix(c(100,0,0,100),2,2) 
cutoff = seq(0.02, 0.7, 0.01)
result = sapply(cutoff, function(p) sum(table(y,pred>p)*payoff) )
i = which.max(result)
par(cex=0.7)
plot(cutoff, result, type='l', col='cyan', lwd=2, main=sprintf(
  "Optomal Expected Result: $%d @ %.2f",result[i],cutoff[i]))
abline(v=seq(0,1,0.05),h=seq(-23000,-17000,500),col='lightgray',lty=3)

abline(v=cutoff[i],col='red')

F6: Simulation

library(manipulate)
p0 = par(mfrow=c(2,1),cex=0.8)
manipulate({
  Y0 = -22000; Y1 = -12000
  mx = matrix(c(true_neg, false_neg, false_pos, true_pos),2,2) 
  cx = seq(0.02, 0.64, 0.01)
  rx = sapply(cx, function(p) sum(table(y, pred>p)*mx) )
  i = which.max(rx)
  plot(cx, rx, type='l',col='cyan',lwd=2,main=sprintf(
    "Optomal Expected Result: $%d @ %.2f, T:%d",rx[i],cx[i],sum(pred>cx[i])),
    ylim=c(Y0,Y1))
  abline(v=cx[i],col='red')
  abline(v=seq(0,1,0.1),h=seq(Y0,Y1,2000),col='lightgray',lty=3)
  DPP(pred, y, 1, b=seq(0,1,0.02))
  abline(v=cx[i],col='red')
  },
  true_neg  = slider(-100,100,0,step=5),
  false_neg = slider(-100,100,-100,step=5),
  false_pos = slider(-100,100,-10,step=5),
  true_pos  = slider(-100,100,-60,step=5)
  ) 
par(p0)

【Q】有五種成本分別為 $5, $10, $15, $20, $30 的藥，它們分別可以將風險成本從 $100 降低到 $70, $60, $50, $40, $25，哪一種藥的期望效益是最大的呢？

#1 $5$70  Expected Result:-17415
#2 $10$60 Expected Result:-17500
#3 $15$50 Expected Result:-17450
#4 $20$40 Expected Result:-17400
#5 $30$25 Expected Result:-17655
#ans: -17400
#哪一種藥的期望效益是最大的呢？($20/$40)這組最大。"

---

【G】分析流程：資料、模型、預測、決策

Figure 4 - 資料、模型、預測、決策