AS3-0 Group-0

getwd()

[1] "/Users/katy/Desktop/rdata"

---

【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 logit & logistic function? What is the relationship between them?

＃logit：一種機率的算法，先算出odd再取log ＃logistic function：把logit轉成機率

---

【B】glm(, family=binomial)

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

#Q %>% head()

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

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

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)

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

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

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? 看觀察點最接近的等機率線為何

---

【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.625,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()
table(y = Q$PoorCare, split) %>% prop.table(2)

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

E2: Build Model

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

E3: Prediction & Evaluation

pred = predict(glm1, type='response')
mx = table(TR$PoorCare, pred > 0.5); mx
c(accuracy = sum(diag(mx))/sum(mx),
  sensitivity = mx[2,2]/sum(mx[2,]),
  specificity = mx[1,1]/sum(mx[1,]))

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)
caTools::colAUC(pred, TR$PoorCare)

---

【F】Framingham Heart Study

getwd()
source("DPP.R")

F1: Reading & Splitting

F <- read.csv("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)

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
c(accuracy = sum(diag(mx))/sum(mx),
  sensitivity = mx[2,2]/sum(mx[2,]),
  specificity = mx[1,1]/sum(mx[1,]))

F4: AUC & DPP

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

F5: Expected Result & Optimization

                                #2*2矩陣
payoff = matrix(c(0,-100,-10,-60),2,2) 
cutoff = seq(0.02, 0.64, 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】如果什麼都不做，期望報酬是多少？ 19700 【Q】如果每位病人都做呢？ 22600 【Q】以上哪一種做法期望報酬比較高？ 什麼都不做 【Q】在所有的商務情境都是這種狀況嗎？ 否 【Q】你可以模擬出「全做」比「全不做」還要好的狀況、並舉出一個會發生這種狀況的商務情境嗎？ (0,-200,-100,500) 預測顧客偏好，-200我們預測該顧客不喜歡此產品，但實際上該顧客喜歡，因此我們會損失沒有該市場的營收。-100是我們針對錯誤的目標市場去做行銷活動的損失成本。500是我們正確預測顧客喜好並對他們做行銷活動。所以在此情況之下，對所有人做行銷活動會比完全不做還好。

#
#
#

F6: Simulation

library(manipulate)
library(ROCR)
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')l
#   },
#   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，哪一種藥的期望效益是最大的呢？

#$5,$10,$15,$20,$30
#-$17415,-$17500,-$17450,-$17400,-$17655
#$20的藥期望效益最大

---

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