點我看教學影片 (ctrl + click) 或是按右鍵開啟新分頁觀看影片
首先將所需要的packages套件下載下來
packages = c(
"dplyr","ggplot2","googleVis","devtools","magrittr","slam","irlba","plotly",
"arules","arulesViz","Matrix","recommenderlab")
existing = as.character(installed.packages()[,1])
for(pkg in packages[!(packages %in% existing)]) install.packages(pkg)
接著載入接下來所有會用到的package套件
rm(list=ls(all=TRUE))
LOAD = TRUE
library(dplyr)
library(ggplot2)
library(googleVis)
library(Matrix)
library(slam)
library(irlba)
library(plotly)
library(arules)
library(arulesViz)
library(recommenderlab)
sed:是當前使用情况 gc trigger:是会触发垃圾回收的值 max used是上次gc()操作或者是此次啟動R後,使用最大值。 (Mb)是Ncells和Vcells的大小轉换為Mb单位時的值。 Ncells即cons cells Vcells即vector cells Load data frame and rename
load("tf0.rdata")
A = A0; X = X0; Z = Z0; rm(A0,X0,Z0); gc()
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2222872 118.8 3886542 207.6 2637877 140.9
## Vcells 8346883 63.7 11198037 85.5 8396695 64.1
Z = subset(Z, cust %in% A$cust)
首先帶大家來看一下我們資料的內容 Z為總交易項目 我們用head(Z)來觀看Z data frame裡的前六筆資料 1. date 是顧客的交易日期 2. cust 是顧客ID 3. age 是年齡代碼 4. area 是顧客所在區域 5. cat 是產品種類代碼 6. prod 是產品條碼 7. qty 是產品數量 8. cost 是已經乘以數量的產品成本 9. price 是乘以數量的金額 10.tid 是交易代碼
如果顧客在同一天買了不同的產品,那在這個Z data frame 中會把他們分開記載。
head(Z)
## date cust age area cat prod qty cost price tid
## 1 2001-01-01 141833 F F 130207 4.710105e+12 2 44 52 58565
## 2 2001-01-01 1376753 E E 110217 4.710266e+12 1 150 129 59070
## 3 2001-01-01 1603071 E G 100201 4.712019e+12 1 35 39 59211
## 4 2001-01-01 1738667 E F 530105 4.710169e+12 1 94 119 59411
## 5 2001-01-01 2141497 A B 320407 4.710431e+12 1 100 159 59936
## 6 2001-01-01 1868685 J E 110109 4.710044e+12 1 144 190 59605
再來看看X資料 這裡是將顧客同一天的消費項目都給予相同的tid,然後以tid來合併顧客在同一天的交易項目 我們觀察一下前十筆資料 1. tid 是交易代碼 2. items 是購買幾種不同種類的產品 3. pieces 是購買數目 4. total 是顧客的總消費金額 5. gross 是毛利,也就是total-總成本
head(X)
## tid date cust age area items pieces total gross
## 1 1 2000-11-01 38317 J E 2 3 76 -8
## 2 2 2000-11-01 45902 H E 4 9 555 95
## 3 3 2000-11-01 45957 G E 1 1 133 -47
## 4 4 2000-11-01 46855 D E 3 5 185 14
## 5 5 2000-11-01 58698 D E 6 6 306 52
## 6 6 2000-11-01 71284 E E 14 14 623 88
最後來看A資料 這邊是以cust,也就是每一位顧客來合併次顧客的所有交易紀錄 1. r是最近一次購買距今天數 2. s是第一次購買距今天數 2. f是購買次數 3. m平均客單價(也就是rev除以購買次數) 4. rev是顧客的花費
head(A)
## cust r s f m rev raw age area
## 1 1069 19 108 4 486.0 1944 15 K E
## 2 1113 54 109 4 557.5 2230 241 K F
## 3 1250 19 25 2 791.5 1583 354 D D
## 4 1359 87 87 1 364.0 364 104 K G
## 5 1823 36 119 3 869.0 2607 498 K D
## 6 2189 57 89 2 7028.0 14056 3299 K B
接著觀察Z這筆交易資料中的顧客人數,和產品人數
n_distinct(Z$cust) # 32241 個顧客
## [1] 32241
n_distinct(Z$prod) # 23787 個產品
## [1] 23787
製作顧客產品矩陣其實很快、也很容易 透過cpm這個函數我們製作一個縱軸為顧客,橫軸為產品的顧客產品矩陣 cpm中會有很多值為零,因為同一個顧客不可購買每一種產品 mean(cpm>0)為0.00096799,代表當中有許多值為零
library(Matrix)
library(slam)
cpm = xtabs(~ cust + prod, Z, sparse=T) # customer product matrix
dim(cpm) # 32241 23787
## [1] 32241 23787
mean(cpm > 0) # 0.00096799
## [1] 0.0009674258
我們觀察cpm數值的分佈,發現有40%的產品的總被購買次數不到六次
colSums(cpm) %>% quantile(seq(0,1,0.1))
## 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
## 1 1 2 4 6 8 13 20 35 76 8475
mean(colSums(cpm) > 10)
## [1] 0.4483541
為了不要讓極端值影響我們的分析,因此刪去購買次數小於6的產品 為了方便顧客產品矩陣和顧客資料框的合併,我們選擇先保留沒有購買產品的顧客
cpm = cpm[, colSums(cpm) >= 6] # remove the least frequent products
# cpm = cpm[rowSums(cpm) > 0, ] # remove non-buying customers
cpm = cpm[, order(-colSums(cpm))] # order product by frequency
dim(cpm) # 32241 23787>14621
## [1] 32241 14621
這裡我們用max(cpm) 發現一個顧客最高購買同一個產品49次 我們在剛剛刪掉總被購買次數小於零的產品後,發現cpm中為0的值減少了 0.00096799 -> 0.0015248 接著觀察每個顧客購買同一產品次數的比例 由此可看出,9成以上的顧客購買同一產品的次數為1次
max(cpm) # 49
## [1] 49
mean(cpm > 0) # 0.0015248
## [1] 0.001524785
table(cpm@x) %>% prop.table %>% round(4) %>% head(10)
##
## 1 2 3 4 5 6 7 8 9 10
## 0.9256 0.0579 0.0108 0.0032 0.0012 0.0006 0.0003 0.0002 0.0001 0.0001
<牛刀小試> 請你用一個指令列出被購買最多次的10個產品,和它們被購買的次數。
cpm[,1:10] %>% col_sums()
## 4714981010038 4711271000014 4719090900065 4711080010112 4710114128038
## 8475 6119 2444 2249 2178
## 4710265849066 4713985863121 4710088410139 4710583996008 4710908131589
## 2017 1976 1869 1840 1679
■ 在什麼前提之下,我們可以把購買這十個產品的次數當作變數,用來預測顧客在下一期會不會來購買呢?
■ 我們如何把這十個變數,併入顧客資料框呢?
■ 我們可不可以(在什麼前提之下我們可以)直接用cbind()
新變數併入顧客資料框呢?
■ 我們期中競賽的資料,符合直接用cbind()
併入新變數的條件嗎? 我們要如何確認這一件事呢?
我們已經依照被購買次數的大小對顧客產品矩陣的欄位做排序 接下來我們以產品的被購買次數來製作變數,排cpm
在最前邊的(N個)欄位就是變數!
在這邊我們選擇前400個最常被購買的產品當作變數對將顧客進行分群,共計200群 其實就是將購買行為相似的顧客分在同一群 最後我們再觀察分群的結果
nop= 400 # no. product = no. variables
k = 200 # no. cluster
set.seed(111); kg = kmeans(cpm[,1:nop], k)$cluster
table(kg) %>% as.vector %>% sort
## [1] 1 1 1 1 1 1 1 1 1 1 1
## [12] 1 1 1 1 1 1 1 1 1 1 1
## [23] 2 2 2 2 2 2 3 3 3 3 3
## [34] 3 3 4 4 4 4 4 4 4 5 5
## [45] 6 6 6 6 7 7 7 8 8 8 9
## [56] 9 9 9 10 10 10 10 11 11 11 11
## [67] 11 12 13 13 15 15 15 16 16 18 19
## [78] 20 20 20 20 21 22 22 22 24 24 25
## [89] 25 27 28 28 32 32 35 36 39 40 41
## [100] 42 44 45 46 47 47 48 49 49 50 51
## [111] 52 53 56 58 58 61 63 66 67 68 69
## [122] 69 72 81 85 85 86 87 90 94 96 97
## [133] 97 100 100 101 110 111 113 114 116 118 123
## [144] 123 126 130 134 136 141 141 142 143 162 165
## [155] 172 175 178 179 182 182 184 187 195 210 222
## [166] 225 228 228 237 239 242 253 254 258 258 266
## [177] 268 272 287 293 301 311 325 329 350 351 363
## [188] 396 407 407 410 418 432 448 473 523 561 1156
## [199] 1266 11215
將分群結果併入顧客資料框A
df = A %>% inner_join(data.frame(
cust = as.integer(rownames(cpm)),
kg) )
## Joining, by = "cust"
head(df) # 32241
## cust r s f m rev raw age area kg
## 1 1069 19 108 4 486.0 1944 15 K E 196
## 2 1113 54 109 4 557.5 2230 241 K F 81
## 3 1250 19 25 2 791.5 1583 354 D D 169
## 4 1359 87 87 1 364.0 364 104 K G 49
## 5 1823 36 119 3 869.0 2607 498 K D 109
## 6 2189 57 89 2 7028.0 14056 3299 K B 158
計算各群組的平均屬性
df = data.frame(
aggregate(. ~ kg, df[,c(2:7,10)], mean), # averages
size = as.vector(table(kg)), # no. customers in the group
dummy = 2001 # dummy column for googleViz
)
head(df)
## kg r s f m rev raw size
## 1 1 10.772727 113.54545 13.318182 872.7878 10697.909 1733.8182 22
## 2 2 43.934959 95.96748 3.455285 1252.9938 3819.837 533.0407 123
## 3 3 5.384615 111.15385 14.307692 1325.2611 12617.769 1951.8462 13
## 4 4 8.000000 119.00000 24.000000 955.6250 22935.000 3658.0000 1
## 5 5 6.666667 116.66667 20.333333 579.9135 10741.444 1554.2222 9
## 6 6 2.500000 114.00000 23.500000 411.9918 9775.000 1324.5000 2
## dummy
## 1 2001
## 2 2001
## 3 2001
## 4 2001
## 5 2001
## 6 2001
這部分的用法可以參考影片
互動式泡泡圖請先點選左上角的查看網站資訊 >> 網站設定 >> flash >> 允許 >> 然後重新整理即可
op = options(gvis.plot.tag='chart')
plot( gvisMotionChart(
subset(df[,c(1,4,5,6,8,2,3,7,9)],
size >= 20 & size <= 1000), # range of group size
"kg", "dummy", options=list(width=800, height=600) ) )
在這邊我們定義了Sig 這個function,我們可以透過在這個function裡面輸入集群的號碼來得知該集群的代表性商品 name # 產品名稱 share # 該產品有多少%被賣給這個集群的顧客 conf # 這個集群的人有多少%的機率購買這個產品 base # 一般顧客購買這個產品的機率 lift # 該集群的人購買這個產品的機率相較於一般顧客購買該產品機率的倍數
# use global variables: cpm, kg
Sig = function(gx, P=1000, H=10) {
print(sprintf("Group %d: No. Customers = %d", gx, sum(kg==gx)))
bx = cpm[,1:P]
data.frame(n = col_sums(bx[kg==gx,])) %>% # frequency
mutate(
share = round(100*n/col_sums(bx),2), # %prod sold to this cluster
conf = round(100*n/sum(kg==gx),2), # %buy this product, given cluster
base = round(100*col_sums(bx)/nrow(bx),2), # %buy this product, all cust
lift = round(conf/base,1), # conf/base
name = colnames(bx) # name of prod
) %>% arrange(desc(lift)) %>% head(H)
}
Sig(130)
## [1] "Group 130: No. Customers = 97"
## n share conf base lift name
## 1 29 6.65 29.90 1.35 22.1 4719090106009
## 2 19 6.62 19.59 0.89 22.0 4710632001172
## 3 12 5.38 12.37 0.69 17.9 4710093080044
## 4 389 4.59 401.03 26.29 15.3 4714981010038
## 5 20 4.59 20.62 1.35 15.3 4719090105002
## 6 14 4.53 14.43 0.96 15.0 4719090106016
## 7 8 4.40 8.25 0.56 14.7 9310022862601
## 8 10 4.31 10.31 0.72 14.3 4710189820851
## 9 9 4.25 9.28 0.66 14.1 4710088434258
## 10 15 4.19 15.46 1.11 13.9 4710189851282
巨大尺度縮減,就是將顧客所有的產品資訊,壓縮成400個特徵相量,這400個特徵向量夾帶原來所有產品的資訊, 換句話說,觀察這400個特徵向量就可以了解顧客購買所有產品的行為。
library(irlba)
if(LOAD) {
load("svd2a.rdata")
} else {
smx = cpm
smx@x = pmin(smx@x, 2) # cap at 2, similar to normalization
t0 = Sys.time()
svd = irlba(smx,
nv=400, # length of feature vector nv=400個特徵向量
maxit=800, work=800)
print(Sys.time() - t0) # 1.8795 mins
save(svd, file = "svd2a.rdata")
}
<牛刀小試2> ■ 在什麼前提之下,我們可以把顧客購買產品的特徵向量當作變數,用來預測顧客在下一期會不會來購買呢?
1.可以。 2.如果要將產品的特徵向量(X)去預測(Y),(X)必須要對預測(Y)有意義,如X與Y要有相關性、預測力或影響力。 ■ 如果可以的話,我們如何把顧客購買產品的特徵向量,併入顧客資料框呢?
1.從svd資料中,找出每位顧客對於400個產品的特徵向量 2.再以CID併入資料集A中 ■ 我們可不可以(在什麼前提之下我們可以)直接用cbind()
將特徵向量併入顧客資料框呢?
1.cbind()把矩陣橫向合併成一個大矩陣(列方式) 2.可以。前提是,row數目要相同且二個資料框cid 的次序要相同。 ■ 我們期中競賽的資料,符合直接用cbind()
併入特徵向量的條件嗎? 我們要如何確認這一件事呢?
1.符合。 2.要確認資料框row的次序要相同。
set.seed(111); kg = kmeans(svd$u, 200)$cluster
table(kg) %>% as.vector %>% sort
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [15] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [43] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [57] 1 2 2 2 2 2 3 4 4 5 7 10 14 30
## [71] 31 32 36 38 38 39 39 40 40 41 44 45 46 47
## [85] 49 54 59 62 62 69 71 77 79 79 80 82 82 84
## [99] 87 91 101 103 109 110 111 113 117 120 123 127 127 129
## [113] 132 133 134 135 136 139 141 143 143 147 147 157 159 159
## [127] 160 160 160 166 168 169 172 175 180 181 181 182 183 184
## [141] 184 188 190 190 193 194 195 196 198 198 200 201 201 202
## [155] 202 204 204 204 207 209 209 210 213 214 216 219 219 222
## [169] 225 233 234 235 236 237 237 238 239 241 248 248 248 253
## [183] 256 257 258 259 261 261 264 269 277 281 285 293 305 411
## [197] 612 896 1092 8987
計算每一群特徵向量的平均屬性
# clustster summary
df = inner_join(A, data.frame(
cust = as.integer(rownames(cpm)), kg)) %>%
group_by(kg) %>% summarise(
avg_frequency = mean(f),
avg_monetary = mean(m),
avg_revenue_contr = mean(rev),
group_size = n(),
avg_recency = mean(r),
avg_gross_profit = mean(raw)) %>%
ungroup %>%
mutate(dummy = 2001, kg = sprintf("G%03d",kg)) %>%
data.frame
## Joining, by = "cust"
head(df)
kg avg_frequency avg_monetary avg_revenue_contr group_size avg_recency
1 G001 4.088050 1205.9121 4103.918 159 39.89308 2 G002 4.517413 1234.0149 4326.552 201 31.82587 3 G003 4.323810 1087.2949 3462.862 210 34.30476 4 G004 5.684211 748.2192 3646.874 190 27.86842 5 G005 4.120000 1153.7658 3640.653 225 33.05778 6 G006 9.142857 1338.6646 9417.545 77 17.90909 avg_gross_profit dummy 1 543.9434 2001 2 712.5224 2001 3 556.2857 2001 4 535.7842 2001 5 581.5422 2001 6 1386.2597 2001
# Google Motion Chart
plot( gvisMotionChart(
subset(df, group_size >= 20 & group_size <= 1200),
"kg", "dummy", options=list(width=800, height=600) ) )
Sig(162)
## [1] "Group 162: No. Customers = 87"
## n share conf base lift name
## 1 165 38.64 189.66 1.32 143.7 4710105045313
## 2 133 35.37 152.87 1.17 130.7 4710105045320
## 3 59 13.50 67.82 1.36 49.9 4710105045450
## 4 19 9.60 21.84 0.61 35.8 4710105051147
## 5 35 8.31 40.23 1.31 30.7 4710105045443
## 6 9 5.33 10.34 0.52 19.9 4710088432650
## 7 8 4.91 9.20 0.51 18.0 4710105051185
## 8 8 4.19 9.20 0.59 15.6 4710171028326
## 9 7 3.89 8.05 0.56 14.4 4710058137059
## 10 7 3.70 8.05 0.59 13.6 4710010020061
太少被購買的產品和購買太少產品的顧客都不適合使用Collaborative Filtering這種產品推薦方法,所以我們先對顧客和產品做一次篩選
library(recommenderlab)
rx = cpm[, colSums(cpm > 0) >= 50]
### 這邊是要50次以上都是由不同顧客購買,所以在cpm>0時,不論單一顧客購買幾次,都為顯示為true,計為一次
rx = rx[rowSums(rx > 0) >= 20 & rowSums(rx > 0) <= 300, ]
### 同理,這裡是討論顧客購買至少購買20種以上300種以下的品項,所以在rx>0時,如果一個顧客購買同一種品項很多次,
### 也會只被計為一次
dim(rx) # 8846 3354 有8846名顧客,購買只少20種以上300種以下品項,3354個產品至少被50個以上不同的顧客購買。
## [1] 8846 3354
可以選擇要用
做模型。
rx = as(rx, "realRatingMatrix") # realRatingMatrix
bx = binarize(rx, minRating=1) # binaryRatingMatrix
UBCF:User Based Collaborative Filtering UBCF,是將購買行為相似的顧客作分群,若你為該群的顧客,會推薦給你群裡其他顧客所購買的產品 舉例來說,Netflix,會將觀看電影相似的客群分成一群,然後推薦群內的顧客彼此看過的影片
rUBCF = Recommender(bx[1:8800,], method = "UBCF")
pred = predict(rUBCF, bx[8801:8846,], n=4)
do.call(rbind, as(pred, "list")) %>% head(15)
## [,1] [,2] [,3] [,4]
## 2170855 "4711271000014" "4710114128038" "4714981010038" "4713985863121"
## 2171265 "4719090900065" "4710254049521" "4710036008562" "4714981010038"
## 2171340 "723125488040" "723125488064" "723125485032" "4714981010038"
## 2171425 "4710011401135" "4710011409056" "4711080010112" "4710011401142"
## 2171432 "4714981010038" "4710011406123" "4711258007371" "4710011401128"
## 2171555 "4719090900065" "4711271000014" "37000329169" "4710943109352"
## 2171883 "4711271000014" "4710583996008" "4710291112172" "4710018004704"
## 2172194 "4711271000014" "4714981010038" "4710114128038" "4710114105046"
## 2172392 "4903111345717" "4710908131589" "4710168705056" "4711271000014"
## 2172569 "4711271000014" "4714981010038" "4710128030037" "4712162000038"
## 2172583 "4714981010038" "4710085120093" "4719090900065" "4710154015206"
## 2172590 "4710011406123" "4710011401142" "4710857000028" "4710011432856"
## 2172668 "4711271000014" "4710088620156" "4719090900065" "4712425010712"
## 2172705 "4711271000014" "4714981010038" "37000445111" "37000440192"
## 2172811 "4714981010038" "37000442127" "4710719000333" "4710114128038"
IBCF:Item Based Collaborative Filtering 是將特徵相似的產品做分群,如果顧客使用這群產品中的產品, 將會推薦給顧客群裡其他的產品 舉裡來說,Netflix,將影片類型細分,像是愛情片、恐怖片等等,若你觀看愛情片,Netflix會推薦給你其他相似的愛情片
if(LOAD) {
load("recommenders.rdata")
} else{
rIBCF <- Recommender(bx[1:6000,], method = "IBCF")
}
pred = predict(rIBCF, bx[8801:8846,], n=4)
do.call(rbind, as(pred, "list")) %>% head(15)
## [,1] [,2] [,3] [,4]
## 2170855 "4719090900065" "4714981010038" "4711271000014" "4712162000038"
## 2171265 "4719090900065" "4710015103288" "4714981010038" "4711271000014"
## 2171340 "37000445111" "4710036005608" "37000442127" "723125485032"
## 2171425 "4711311617899" "4711311218836" "4710011401135" "4710011409056"
## 2171432 "4714981010038" "4710321791698" "4710857000042" "4710626111252"
## 2171555 "93432641" "93362993" "4710105045320" "4711271000014"
## 2171883 "4710670200100" "4710670200407" "4711271000014" "3228020490329"
## 2172194 "4714108700019" "4714108700064" "4909978199111" "20332433"
## 2172392 "4710706211759" "4710908131589" "4719090900058" "4710731040614"
## 2172569 "4711371850243" "84501297329" "84501293529" "4710085121007"
## 2172583 "4710085172702" "4710085120093" "34000100095" "34000231508"
## 2172590 "4710011406123" "4711271000014" "4711437000162" "4710011401142"
## 2172668 "4711371850243" "4719090900065" "4714981010038" "4711437000117"
## 2172705 "37000445111" "4710018004605" "37000442127" "37000304593"
## 2172811 "4719581980293" "4712067899287" "4719581980279" "4710908131589"
save(rIBCF, rUBCF, file="recommenders.rdata")
set.seed(4321)
scheme = evaluationScheme(
bx, method="split", train = .75, given=5)
<參考購物籃分析中的關聯規則> support 該產品被全部顧客購買的機率 confidence A被購買時B被購買的機率 Random:就是用隨機的方式來推薦你可能會喜歡的產品 Popular:推薦大家都喜歡的熱門產品給你 關聯規則:我們找出當顧客購買A產品時,購買B產品的機率,藉此找出A、B產品之間的關聯規則,然後利用這些規則來做推薦
algorithms = list(
AR53 = list(name="AR", param=list(support=0.0005, confidence=0.3)),
AR43 = list(name="AR", param=list(support=0.0004, confidence=0.3)),
RANDOM = list(name="RANDOM", param=NULL),
POPULAR = list(name="POPULAR", param=NULL),
UBCF = list(name="UBCF", param=NULL),
IBCF = list(name="IBCF", param=NULL) )
if(LOAD) {
load("results2a.rdata")
} else {
t0 = Sys.time()
results = evaluate(
scheme, algorithms,
type="topNList", # method of evaluation
n=c(5, 10, 15, 20) # no. recom. to be evaluated
)
print(Sys.time() - t0)
save(results, file="results2a.rdata")
}
## AR run fold/sample [model time/prediction time]
## 1 [4.02sec/214.6sec]
## AR run fold/sample [model time/prediction time]
## 1 [10.49sec/538.5sec]
## RANDOM run fold/sample [model time/prediction time]
## 1 [0sec/9.48sec]
## POPULAR run fold/sample [model time/prediction time]
## 1 [0sec/11.09sec]
## UBCF run fold/sample [model time/prediction time]
## 1 [0sec/75.42sec]
## IBCF run fold/sample [model time/prediction time]
## 1 [198.2sec/1.63sec]
## Time difference of 18.72 mins
# load("results.rdata")
par(mar=c(4,4,3,2),cex=0.8)
cols = c("red", "magenta", "gray", "orange", "blue", "green")
plot(results, annotate=c(1,3), legend="topleft", pch=19, lwd=2, col=cols)
abline(v=seq(0,0.006,0.001), h=seq(0,0.08,0.01), col='lightgray', lty=2)
getConfusionMatrix(results$IBCF)
## [[1]]
## TP FP FN TN precision recall TPR
## 5 1.115732 3.884268 32.97288 3311.027 0.2231465 0.03898910 0.03898910
## 10 1.699367 8.300633 32.38924 3306.611 0.1699367 0.05811506 0.05811506
## 15 2.075045 12.924955 32.01356 3301.986 0.1383363 0.07021014 0.07021014
## 20 2.385172 17.614828 31.70344 3297.297 0.1192586 0.08002344 0.08002344
## FPR
## 5 0.001171448
## 10 0.002503393
## 15 0.003898163
## 20 0.005312770