pacman::p_load(caTools, ggplot2, dplyr)
D = read.csv("data/quality.csv") # Read in dataset
set.seed(88)
split = sample.split(D$PoorCare, SplitRatio = 0.75) # split vector
TR = subset(D, split == TRUE)
TS = subset(D, split == FALSE)
glm1 = glm(PoorCare ~ OfficeVisits + Narcotics, TR, family=binomial)
summary(glm1)
因為這個資料集很小,我們使用全部的資料來做模擬 (通常我們是使用測試資料集)
pred = predict(glm1, D, type="response")
y = D$PoorCare
data.frame(pred, y) %>%
ggplot(aes(x=pred, fill=factor(y))) +
geom_histogram(bins=20, col='white', position="stack", alpha=0.5) +
ggtitle("Distribution of Predicted Probability (DPP,FULL)") +
xlab("predicted probability")
報酬矩陣 Payoff Matrix
payoff = matrix(c(0,-100,-10,-50),2,2)
rownames(payoff) = c("FALSE","TRUE")
colnames(payoff) = c("NoAct","Act")
payoff
NoAct Act
FALSE 0 -10
TRUE -100 -50
期望報酬 Expected Payoff
cutoff = seq(0, 1, 0.01)
result = sapply(cutoff, function(p) {
cm = table(factor(y==1, c(F,T)), factor(pred>p, c(F,T)))
sum(cm * payoff) # sum of confusion * payoff matrix
})
i = which.max(result)
par(cex=0.7, mar=c(4,4,3,1))
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.1),h=seq(-6000,0,100),col='lightgray',lty=3)
points(cutoff[i], result[i], pch=20, col='red', cex=2)
使用manipulate
套件做策略模擬
library(manipulate)
manipulate({
payoff = matrix(c(TN,FN,FP,TP),2,2)
cutoff = seq(0, 1, 0.01)
result = sapply(cutoff, function(p) {
cm = table(factor(y==1, c(F,T)), factor(pred>p, c(F,T)))
sum(cm * payoff) # sum of confusion * payoff matrix
})
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.1),h=seq(-10000,0,100),col='lightgray',lty=3)
points(cutoff[i], result[i], pch=20, col='red', cex=2)
},
TN = slider(-100,0, 0,step=5),
FN = slider(-100,0,-100,step=5),
FP = slider(-100,0, -10,step=5),
TP = slider(-100,0, -50,step=5)
)
🗿 練習:
執行Sim12.R
,先依預設的報酬矩陣回答下列問題:
【A】 最佳臨界機率是? 它所對應的期望報酬是多少?
【B】 什麼都不做時,臨界機率和期望報酬各是多少?
【C】 每位保戶都做時,臨界機率和期望報酬各是多少?
【D】 以上哪一種做法期的望報酬比較高?
【E】 在所有的商務情境都是這種狀況嗎?
藉由調整報酬矩陣:
【F】 模擬出「全不做」比「全做」還要好的狀況
【G】 並舉出一個會發生這種狀況的商務情境
有五種成本分別為$5, $10, $15, $20, $35
的介入方法,它們分別可以將風險成本從$100
降低到$70, $60, $50, $40, $30
…
【H】 它們的最佳期望報酬分別是多少?
【I】 哪一種介入方法的最佳期望報酬是最大的呢?
🌞 參考模擬程式
pacman::p_load(plotly)
PMX = list(
c(0,-100, -5,-75), c(0,-100,-10,-70), c(0,-100,-15,-65),
c(0,-100,-20,-60), c(0,-100,-35,-65)) %>%
lapply(matrix, nrow=2, ncol=2)
do.call(rbind, lapply(1:length(PMX), function(i) {
r = sapply(cutoff, function(p) {
cm = table(factor(y==1,c(F,T)), factor(pred>p,c(F,T)))
sum(cm * PMX[[i]]) })
data.frame(drug=paste0("drug_",i), cutoff=cutoff, exp.payoff=r)
})) %>%
ggplot(aes(x=cutoff, y=exp.payoff, col=drug)) +
geom_line() + theme_bw() -> p
ggplotly(p)