ML-tutorial
GSTAT
מבוא
בשיעורים הקודמים למדנו את השפה הבסיסית של R , איך להשתמשת בספריות חיצוניות, ואיך לבצע שאילתות ומנפולציה על נתונים בעזרת חבילות מה tidyverse.
בשיעור הזה נלמד איך לבצע ניתוחים סטטיסטים בסיסיים ומתקדמים יותר בעזרת מגוון ספריות ב R.
השיעור חולק לפי נושאים:
- סטטיסטיקה תאורית ב R
- נושאים בבניית מודלים
- אקונומטריקה קלאסית
- שיטות ב Machine Learning
מטרות השיעור
בשיעור לא נלמד על אף אחד מהנושאים לעומק.
המטרה היא להכיר לכם את הפרקטיקה של כל הנושאים שנלמד, ולהכיר לכם את השיטות והחבילות הנפוצות ביותר כיום לנושאים האלה.
הכוונה היא שעד סוף השיעור, תדעו:
- איך לבצע ניתוח שונות ב R
- איך להכין את הנתונים לבניית מודלים
- עקרונות בסיסיים בבניית מודלים
- איך להריץ רגרסיות ב R
- איך לבצע ניתוחים שונים של ML ב R
במהלך השיעור נעשה שימוש נרחב בחבילת CARET
בשביל שהניתוחים ירוצו יותר מהר נבחר תת מדגם מאוזן שהכנתי מראש עם פקודה שנלמד עליה עוד מעט…
library(tidyverse)
library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftmodel_data <- Santander_sample %>% filter(!is.na(ind_cc))
sub_sample <- createDataPartition(y = model_data$ind_cc,
p = .05, list = FALSE)
model_data <- Santander_sample[sub_sample,]
rm(sub_sample, Santander_sample)סטטיסטיקה תאורית ב R
בשיעורים הקודמים כבר למדתם איך להציג סטטיסטיקה תאורית לפי קבוצות של משתנים בעזרת הפקודה summary , table , ושאילתות בחבילת dplyr.
בחלק נשלים את התמונה בעזרת הצגת מבחנים סטטיסטיים לבדיקת שונות בין ממוצעים של משתנה מטרה בין קבוצות של משתנים מסבירים, ובדיקת קורלציות בין משתנים.
ניתוח שונות
ניתוח ANOVA נעשה בעזרת הרצת רגרסיה ליניארית ובדיקת השונות בחלקי המדגם.
הפקודה הכללית להרצת רגרסיה ליניארית היא lm
lm(ind_cc~renta, data = model_data) %>%
anova() ## Analysis of Variance Table
##
## Response: ind_cc
## Df Sum Sq Mean Sq F value Pr(>F)
## renta 1 0.13 0.128003 10.948 0.0009372 ***
## Residuals 82650 966.30 0.011691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1בדיקת התפלגויות
table, prop.test
ראינו כבר להוציא התפלגויות כלליות .
בעזרת הפקודה tableו prop.test נוכל להוסיף מבחני מובהקות לשוני באוכלוסיה
table(model_data$segmento,model_data$ind_cc) %>%
prop.test()##
## 3-sample test for equality of proportions without continuity
## correction
##
## data: .
## X-squared = 340.94, df = 2, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## prop 1 prop 2 prop 3
## 0.9696640 0.9868407 0.9969994קורילציות בין משתנים
הפקודה העיקרית היא cor כאשר בשיטה ניתן לבחור בין אחד מהשיטות הבאותץ
ניתן
- kendall
- pearson
- spearman
corr_table <- keep(model_data, is.numeric) %>%
cor(.,method = "pearson" , use ="complete.obs")## Warning in cor(., method = "pearson", use = "complete.obs"): the standard
## deviation is zerocorr_table[is.na(corr_table)] <- 0
library(d3heatmap)
d3heatmap(corr_table)rm(corr_table)התוצר של פקודת ה cor היא מטריצה אם היינו רוצים לעבוד עם ggplot על הקורילציות היינו צריכים לעבוד קשה ולהמיר את זה לטבלה
במקום זה נשתמש בחבילת d3heatmap שלוקחת מטריצה ומחזירה מפת חום אינטראקטיבית
נושאים בבניית מודלים
בניית מודלים סטטיסטיים ניתן לחלק למרכיבים הבאים:
- חלוקה למדגמים
- הכנה מוקדמת של הנתונים
- הרצת המודלים בהתאם לשיטה סטטיסטית מסוימת
- בחינת תוצאות טיב התאמה של המודל במדגמים השונים
חלוקה למדגמים
ה best practice הוא לחלק את הנתונים לקבוצת פיתוח (train) וקבוצת ביקורת (test), בד“כ ביחס של 70, 30.
זאת על מנת למזער את נזקי over fitting ועוד.
לעיתים גם נרצה להוציא גם מדגם של Out of Time על מנת לראות איך המודל חוזה קדימה.
אנו נעשה זאת בעזרת חבילת CARET והפקודה createDataPartition.
הפקודה מחזירה וקטור של אינדקסים שלפיהן נחתוך את הנתונים שלנו.
library(caret)
count(model_data, ind_cc)## # A tibble: 3 × 2
## ind_cc n
## <int> <int>
## 1 0 81674
## 2 1 978
## 3 NA 26# Only 33 missing assuming 0
model_data <- model_data %>% mutate(ind_cc = ifelse(is.na(ind_cc),0,ind_cc))
# caret::createDataPartition returns indices for rows based on a stratafied sample over the target variable y
# y = target variable
# p = percent obs
# list = format returned
inTrain <- createDataPartition(y = model_data$ind_cc, p = .7, list = FALSE)
# we use the indices to split our data into training (In sample) and testing (out of sample) tables
model_train <- model_data[inTrain,]
model_test <- model_data[-inTrain,]
rm(inTrain)טיפול בנתונים חסרים
בגדול יש שתי סיבות למה יהיו נתונים חסרים:
- הנתונים חסרים באקראי
- יש תבנית לנתונים החסרים
הטיפול בנתונים החסרים יעשה בהתאם לסיבת החוסר בנתונים
כמובן שאם הנתונים יותר מדי חסרים אז אולי לא כדאי להשתמש במשתנה הזה…
אפשר לקבל רושם כללי לכמות המשתנים החסרים בעזרת summary
- אפשר לחקור או להשלים נתונים חסרים בעזרת חבילת mice
- md.pattern - מחזיר דפוסים של חוסר בין מספר משתנים
- mice - משלים נתונים חסרים לפי שיטה נבחרת
library(mice)## Loading required package: Rcpp## mice 2.25 2015-11-09##
## Attaching package: 'mice'## The following object is masked from 'package:tidyr':
##
## completemd.pattern(model_train[,50:55])## ind_cc lead_3m_cc last_1m_cc last_6m_cc last_3m_cc
## 57669 1 1 1 1 1
## 16 1 0 1 1 1
## 190 1 1 0 0 0
## 0 16 190 190 190
## l1m_ind_ahor_fin_ult1
## 57669 1 0
## 16 1 1
## 190 0 4
## 190 776summary(model_data[50:55])## ind_cc last_1m_cc lead_3m_cc last_6m_cc
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000 Median :0.00000 Median :0.00000
## Mean :0.01183 Mean :0.01211 Mean :0.03861 Mean :0.09453
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.00000
## Max. :1.00000 Max. :1.00000 Max. :3.00000 Max. :6.00000
## NA's :263 NA's :26 NA's :263
## last_3m_cc l1m_ind_ahor_fin_ult1
## Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000
## Mean :0.04697 Mean :0.00012
## 3rd Qu.:0.00000 3rd Qu.:0.00000
## Max. :3.00000 Max. :1.00000
## NA's :263 NA's :263# imputing missing data
#############################
# This example imputes the missing data using `predictive mean method`
# but there are other options: see ?mice
# A `mids` object is returned.
# We need to use `complete` to make the imputed data
tempData <- mice(model_data[,50:53]
,m=5
,maxit=3
,meth='pmm'
,seed=500)##
## iter imp variable
## 1 1 last_1m_cc lead_3m_cc last_6m_cc
## 1 2 last_1m_cc lead_3m_cc last_6m_cc
## 1 3 last_1m_cc lead_3m_cc last_6m_cc
## 1 4 last_1m_cc lead_3m_cc last_6m_cc
## 1 5 last_1m_cc lead_3m_cc last_6m_cc
## 2 1 last_1m_cc lead_3m_cc last_6m_cc
## 2 2 last_1m_cc lead_3m_cc last_6m_cc
## 2 3 last_1m_cc lead_3m_cc last_6m_cc
## 2 4 last_1m_cc lead_3m_cc last_6m_cc
## 2 5 last_1m_cc lead_3m_cc last_6m_cc
## 3 1 last_1m_cc lead_3m_cc last_6m_cc
## 3 2 last_1m_cc lead_3m_cc last_6m_cc
## 3 3 last_1m_cc lead_3m_cc last_6m_cc
## 3 4 last_1m_cc lead_3m_cc last_6m_cc
## 3 5 last_1m_cc lead_3m_cc last_6m_cc############################
# Returning the compete data
completedData <- complete(tempData,1)
summary(completedData)## ind_cc last_1m_cc lead_3m_cc last_6m_cc
## Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000 Median :0.0000 Median :0.00000
## Mean :0.01183 Mean :0.01208 Mean :0.0386 Mean :0.09433
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.00000
## Max. :1.00000 Max. :1.00000 Max. :3.0000 Max. :6.00000####
rm(tempData, completedData)טרנספורמציות על נתונים
כמעט בכל השיטות הסטטיסטיות נרצה לבצע הכנה מוקדמת.
שתי השיטות הנפוצות ביותר הן:
- מיקוד ונרמול נתונים
- חלוקת משתנים לקטגוריות
מיקוד ונרמול
בעזרת פקודת preProcess של חבילת CARET ניתן לבצע כל מיני טרנספורמציות נפוצות על הנתונים
ניתן להשתמש בפקודה כחלק מפקודת test המריצה רגרסיות ושיטות של ML או כפקודה בפני עצמה
כאן נדגים שימוש של הפקודה לנמרמל ולמרכז את הנתונים
הסינטקס הכללי של הפקודה היא:
preProcess(data, method, parameters)
# Do the centering
trainX <- model_data %>%
select(renta, age, ind_cc) %>%
preProcess(method = c("center", "scale"))
trainX## Created from 82678 samples and 3 variables
##
## Pre-processing:
## - centered (3)
## - ignored (0)
## - scaled (3)# apply the centering
model_data_centered <- model_data %>%
select(renta,age, ind_cc) %>%
predict(trainX,.)
head(model_data_centered, 5)## renta age ind_cc
## 1 1.59818413 0.9745756 -0.1094097
## 2 -0.06913145 0.9745756 -0.1094097
## 3 0.05844565 0.0628401 -0.1094097
## 4 -0.12466046 0.1236225 -0.1094097
## 5 -0.59422756 0.4275343 -0.1094097- עוד אפשרות זה להשתמש בפקודה
scaleב R BASE
חלוקה לקטגוריות
בגדול יש שתי סוגי משתנים:
- משתנים רציפים
- משתנים קטגוריאלים
יש כל מיני סיבות למה נרצה לאגד משתנים לקטגוריות:
- חיסון מפני ערכים קיצוניים
- יותר מתאים לרגרסיה לוגיסטית
- יצירת קבוצות הומוגניות מבחינת משתנה המטרה
- ועוד
משתנים רציפים
יש מספר דרכים לחלק משתנים רציפים לקבוצות משמעותיות ביחס למשתנה מטרה
אנחנו נציג שיטה גרפית פשוטה:
model_data %>%
filter(renta<500000) %>%
ggplot(.,aes(x=renta)) +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
model_data %>%
filter(renta<500000) %>%
ggplot(.,aes(x=renta, y= ind_cc)) +
geom_smooth()## `geom_smooth()` using method = 'gam'
משני הגרפים לעיל רואים גם:
- בעזרת ההיסטוגרמה נדאג שהקבוצות לא יהיו קטנות מדי
- בעזרת ה loess נחתוך לקבוצות עם ממוצע משתנה מטרה דומה
- יש גרעין גדול ב 0
- בקושי יש ערכים מעל 300,000
בדוגמא הייתי בוחר לחתוך בערך ב - 50,000 100,000 ו 300,000
אפשר לחתוך בעזרת cut אבל מהתרגיל יש לכם כבר כלי חזק לחתוך משתנים לקטגוריות.
########################################################
# First lets make a format function maker
# Answer to question 6 format function:
grp_fmt<- function(...){
function(data_var){
data_var <- cut(data_var
, include.lowest=TRUE
, right=TRUE
, dig.lab = 7
,breaks=c(min(data_var), ... , max(data_var)))
data_var
}
}יצרנו מעטפת לפקודת cut , המאפשרת חיתוך מהיר ושמירת התיעוד של החיתוך בתצורה של פונקציה
הפונקציה מחזירה פונקציה הדומה לפורמט של SAS
ניתן להשתמש בה כמעט בכל מקום…
אבל יש לה מקום לשיפור… איך הייתם מתקנים אותה כדי שתטפל בערכים חסרים?
##############################################
# What if Im not sure where to cut?
# I can create 2 options
fmt_renta1 <- grp_fmt(50000, 150000,250000)
fmt_renta2 <- grp_fmt(50000, 100000,300000)
##############################################
# And check them one by one like this:
model_train %>%
group_by(fmt_renta1(renta)) %>%
summarise(pct = mean(ind_cc, na.rm = T), Freq = n()) %>%
ungroup() %>%
mutate(grp_pct=Freq*100/sum(Freq))## # A tibble: 4 × 4
## `fmt_renta1(renta)` pct Freq grp_pct
## <fctr> <dbl> <int> <dbl>
## 1 [0,50000] 0.008731587 15003 25.923110
## 2 (50000,150000] 0.011127419 29297 50.621166
## 3 (150000,250000] 0.016142919 9292 16.055292
## 4 (250000,2.203474e+07] 0.018911978 4283 7.400432##############################################
# I can store them in a list for future use:
fmtl_renta <- list( "fmt_renta1" = fmt_renta1
,"fmt_renta2" = fmt_renta2)
names(fmtl_renta)## [1] "fmt_renta1" "fmt_renta2"##############################################
# Or see them all at once with map or lapply:
map_df(fmtl_renta,
~ model_train %>%
group_by(.(renta)) %>%
summarise(pct = mean(ind_cc, na.rm = T),Freq = n()) %>%
ungroup() %>%
mutate(grp_pct=Freq*100/sum(Freq))
)## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character## # A tibble: 8 × 4
## `.(renta)` pct Freq grp_pct
## <chr> <dbl> <int> <dbl>
## 1 [0,50000] 0.008731587 15003 25.923110
## 2 (50000,150000] 0.011127419 29297 50.621166
## 3 (150000,250000] 0.016142919 9292 16.055292
## 4 (250000,2.203474e+07] 0.018911978 4283 7.400432
## 5 [0,50000] 0.008731587 15003 25.923110
## 6 (50000,100000] 0.009227019 17232 29.774514
## 7 (100000,300000] 0.014888337 22971 39.690713
## 8 (300000,2.203474e+07] 0.020981641 2669 4.611663הדוגמה לעיל ממחישה את הגמישות המובנית בכתיבה פונקציונאלית
משתנים קטגוריאלים
בשביל לבחון משתנה קטגוריאלי מספיק להסתכל על הממוצעים של אחוזי המשתנה הבינומי ביחס למספר התצפיות בכל קטגוריה
לרוב לא נרצה קבוצות קטנות מדי
לדוגמא במשתנה הבא רואים שיש קבוצה אחת קטנה מדי שכדאי לאחד עם אחד הקבוצות הגדולות
model_train %>%
group_by(tiprel_1mes) %>%
summarise(pct = mean(ind_cc, na.rm = T)
,Freq = n()) %>%
ungroup() %>%
mutate(grp_pct=Freq*100/sum(Freq))## # A tibble: 3 × 4
## tiprel_1mes pct Freq grp_pct
## <chr> <dbl> <int> <dbl>
## 1 A 0.0133038125 51489 88.965874730
## 2 I 0.0004699248 6384 11.030669546
## 3 P 0.0000000000 2 0.003455724נציג שתי אופציות לחתוך את המשתנה:
- בעזרת gui -> Addins -> Variable Cutting
- בעזרת הפקודה
recode
##########################################################################
# Option 1: Recode
# I made this code with the GUI:
model_train$tiprel_1mes_rec <- model_train$tiprel_1mes
model_train$tiprel_1mes_rec[model_train$tiprel_1mes == "P"] <- "A"
model_train$tiprel_1mes_rec[model_train$tiprel_1mes == ""] <- "A"
# This turns it form a char vector to a factor:
model_train$tiprel_1mes_rec <- factor(model_train$tiprel_1mes_rec)
# Lets see the results
model_train %>%
group_by(tiprel_1mes_rec) %>%
summarise(pct = mean(ind_cc, na.rm = T) ,Freq = n()) %>%
ungroup() %>%
mutate(grp_pct=Freq*100/sum(Freq))## # A tibble: 2 × 4
## tiprel_1mes_rec pct Freq grp_pct
## <fctr> <dbl> <int> <dbl>
## 1 A 0.0133032957 51491 88.96933
## 2 I 0.0004699248 6384 11.03067##########################################################################
# Option 2: Recode
# There is a dplyr function for this called `recode` and
# `recode_factor` (that turns it into an ordered factor), look at the help:
model_train %>%
group_by(recode_factor(tiprel_1mes,P="A")) %>%
summarise(pct = mean(ind_cc, na.rm = T) ,Freq = n()) %>%
ungroup() %>%
mutate(grp_pct=Freq*100/sum(Freq))## # A tibble: 2 × 4
## `recode_factor(tiprel_1mes, P = "A")` pct Freq grp_pct
## <fctr> <dbl> <int> <dbl>
## 1 A 0.0133032957 51491 88.96933
## 2 I 0.0004699248 6384 11.03067כפי שאתם רואים שימוש ב GUI זה נחמד אבל דורש יצירת משתנה חדש וכתיבת הרבה קוד….
מדדים לטיב המודל
מודלים של רגרסיה לינאירת נסתכל ב R בריבוע
לגבי מובהקות משתנים נשתמש ב
drop1העושה מבחןTYPE IIIבדומה ל SASברגרסיה בינומית אפשר להסתכל על מה שהכי מקטין את ה
devianceאת המדדים למודלים בינארים נקח מחבילת ROCR או pROC.
נראה אותם בהמשך
אקונומטריקה קלאסית
Syntax של משוואה
כאן נעבור על ה syntax שילווה אותנו ברוב הפקודות של רגרסיות ועוד.
רגרסיות ב R מוגדרות על ידי formula מהסוג:
y ~ x1 + X2
כאשר + פירושו עוד משתנה ולא חיבור מתמטי
טרנספורמציה בתוך נוסחא
ניתן לבצע כל טרנספורמציה בתוך שורת הנוסחא
y ~ exp(x1) + factor(X2)
- כמובן שאפשר גם להשתמש בפונקציה
grp_frmtשיצרנו גם בתוך שורת הרגרסיה!*
פעולה מתמטית בתוך נוסחא
אם רוצים פעולה מתמטית יש לכתוב זאת בתוך I
y ~ I(x1^2) + I(X2*x1)
אינטראקציה בין משתנים
אינטראקציה בין משתנים נסמן בסימן * כאשר אנטראקציה עם המשתנה עצמו נסמן בהעלאה בחזקה ^
y ~ x1^2 + X2*x1
פולינום על משתנה
במקום להוסיף ידנית את ערכי הפולינום אפשר לסמן פולינום על משתנה בעזרת הפקודה poly
y ~ x1 + poly(X2,4)
אבל יתכן ותרצו במקרה כזה להשתמש ב spline , יש פונקציה בחבילת splines, הפונקציה ns()
רגרסיה בלי חותך
רגרסיה בלי חותך נסמן על ידי חיסור של 1
y ~ x1 -1
עכשיו שאתם יודעים איך לכתוב נוסחא, נריץ כמה מודלים
פקודות נפוצות למודלים
כמעט לכל השיטות המחזירות מודלים יש את הפקודות הבאות, כאשר לכל שיטה יש את הנואנסים שלה:
summary- מציג סיכום תוצאות המודלpredict- מייצר ערכים חזוייםcoefficients- מחזיר את המקדמיםresid- מחזיר את השאריותdrop1- מבחן F למובהקות המשתנה - SAS TYPE IIIמי שרוצה לשדרג את העבודה עם רגרסיות - שיסתכל בחבילת
broom*
רגרסיה ליניארית
בוא נריץ רגרסיה ליניארית על renta
ה syntax הבסיסי:
lm(formula, data)
options(
contrasts = c("contr.treatment", "contr.treatment"),
na.option = na.exclude
)
reg_lm1 <- lm(renta ~ segmento + poly(age,3), data = model_train)
# regression outcomes
summary(reg_lm1)##
## Call:
## lm(formula = renta ~ segmento + poly(age, 3), data = model_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -174519 -65324 -20361 32033 21929374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 143935 3222 44.667 < 2e-16 ***
## segmento02 - PARTICULARES -30553 3281 -9.312 < 2e-16 ***
## segmento03 - UNIVERSITARIO -42334 4336 -9.763 < 2e-16 ***
## poly(age, 3)1 2439184 286018 8.528 < 2e-16 ***
## poly(age, 3)2 853683 233438 3.657 0.000255 ***
## poly(age, 3)3 -1491289 206721 -7.214 5.5e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 205700 on 57869 degrees of freedom
## Multiple R-squared: 0.008568, Adjusted R-squared: 0.008483
## F-statistic: 100 on 5 and 57869 DF, p-value: < 2.2e-16# Type III paramter test
drop1(reg_lm1, test = "F")## Single term deletions
##
## Model:
## renta ~ segmento + poly(age, 3)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 2.4480e+15 1416098
## segmento 2 4.4141e+12 2.4524e+15 1416198 52.173 < 2.2e-16 ***
## poly(age, 3) 3 6.9634e+12 2.4550e+15 1416256 54.870 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# coefficients
coefficients(reg_lm1)## (Intercept) segmento02 - PARTICULARES
## 143935.00 -30553.09
## segmento03 - UNIVERSITARIO poly(age, 3)1
## -42333.61 2439184.11
## poly(age, 3)2 poly(age, 3)3
## 853683.43 -1491289.32# residuals and predictons
resid(reg_lm1) %>% str()## Named num [1:57875] 286275 -59210 14639 -21036 -46092 ...
## - attr(*, "names")= chr [1:57875] "1" "2" "3" "4" ...predict(reg_lm1, model_train) %>% str()## Named num [1:57875] 127841 158394 108642 109731 118429 ...
## - attr(*, "names")= chr [1:57875] "1" "2" "3" "4" ...# predict values and residuals
reg_ml1_data <- model_train %>%
mutate( pred_renta = predict(reg_lm1, model_train)
,resid_renta = resid(reg_lm1)) %>%
select(renta, pred_renta, resid_renta, segmento, age)
######################################
# Lets examine some graphs...
reg_ml1_data %>%
filter(renta < 500000) %>%
ggplot(.,aes(renta, pred_renta)) +
geom_point(aes(color=segmento)) +
geom_smooth(aes(linetype = segmento), method = "lm") +
facet_wrap(~segmento, scales = "free", ncol = 1)
# Well, that looks pretty bad!
# Let see the residuals:
# With respect to segmento:
reg_ml1_data %>%
filter(renta < 500000) %>%
ggplot(.,aes(segmento, resid_renta)) +
geom_boxplot() +
coord_flip()
# With respect to renta:
reg_ml1_data %>%
filter(renta < 500000) %>%
ggplot(.,aes(renta, resid_renta)) +
geom_jitter(aes(color = segmento)) +
facet_wrap(~segmento)
טוב, זה היה דוגמא לרגרסיה גרועה! :)
אבל הבנתם את העקרון…
- לגבי מדדי החלטה ה R2 שלנו כ 0.008
רגרסיה לוגיסטית
רגרסיה לוגיסטית שייכת למשפחת הרגרסיות של glm
ה syntax הבסיסי לרגרסיה לוגיסטית היא:
glm(formula, data, family=binomial(link = “logit”))
reg_log1 <- glm(ind_cc~ segmento + sexo
,model_train
,family=binomial(link = "logit"))
drop1(reg_log1, test = "LRT")## Single term deletions
##
## Model:
## ind_cc ~ segmento + sexo
## Df Deviance AIC LRT Pr(>Chi)
## <none> 7222.3 7230.3
## segmento 2 7455.1 7459.1 232.763 < 2e-16 ***
## sexo 1 7225.7 7231.7 3.363 0.06667 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(reg_log1)##
## Call:
## glm(formula = ind_cc ~ segmento + sexo, family = binomial(link = "logit"),
## data = model_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.2563 -0.1671 -0.1555 -0.0820 3.4182
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.54435 0.10030 -35.338 <2e-16 ***
## segmento02 - PARTICULARES -0.86549 0.09693 -8.929 <2e-16 ***
## segmento03 - UNIVERSITARIO -2.29468 0.16991 -13.505 <2e-16 ***
## sexoV 0.14516 0.07959 1.824 0.0682 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 7466.6 on 57874 degrees of freedom
## Residual deviance: 7222.3 on 57871 degrees of freedom
## AIC: 7230.3
##
## Number of Fisher Scoring iterations: 8# coefficients
head(coefficients(reg_log1))## (Intercept) segmento02 - PARTICULARES
## -3.5443458 -0.8654904
## segmento03 - UNIVERSITARIO sexoV
## -2.2946789 0.1451639# residuals
head(resid(reg_log1))## 1 2 3 4 6 7
## -0.1670844 -0.2563392 -0.1670844 -0.1554597 -0.1554597 -0.1670844# predict values
pred_vals <- predict(reg_log1, model_train, type = "response")
head(pred_vals)## 1 2 3 4 6 7
## 0.01386163 0.03232104 0.01386163 0.01201115 0.01201115 0.01386163שימו לב שב predict השתמשתי ב response
מה הייתי מקבל אם לא הייתי מציין את זה?
בשביל מדד לטיב הניבוי נסתכל על מדד gini שנגזר מחישוב AUC - את החישוב נעשה בעזרת ספריית pROC
הפקודה auc לוקחת ערך בפועל וערך חזוי
library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var###### In Sample
# AUC from ROC
auc(model_train$ind_cc, pred_vals)## Area under the curve: 0.6429# Gini coefficient
auc(model_train$ind_cc, pred_vals)*2-1## [1] 0.2857619###### Out of Sample
pred_vals <- predict(reg_log1, model_test, type = "response")
# AUC from ROC
auc(model_test$ind_cc, pred_vals)## Area under the curve: 0.6573# Gini coefficient
auc(model_test$ind_cc, pred_vals)*2-1## [1] 0.3145218חבילת pROC מאפשר עוד כל מיני גרפים , חפשו בעזרה לעוד אפשרויות
- יש גם גרפים גנריים לכל רגרסיה - אבל זה בעיקר scatterplot שלוקחים הרבה זמן להדפיס….
עוד שימוש שאפשר לעשות עם כתיבה פונקציונאלית זה לבדוק מספר רגרסיות בשורה אחת:
fmrla_list <- list(
frml1 = ind_cc~ segmento + sexo
,frml2 = ind_cc~ segmento + sexo + renta
)
fmrla_list %>%
map(. ,function(x) glm(x
,model_train
,family=binomial(link = "logit"))
)## $frml1
##
## Call: glm(formula = x, family = binomial(link = "logit"), data = model_train)
##
## Coefficients:
## (Intercept) segmento02 - PARTICULARES
## -3.5443 -0.8655
## segmento03 - UNIVERSITARIO sexoV
## -2.2947 0.1452
##
## Degrees of Freedom: 57874 Total (i.e. Null); 57871 Residual
## Null Deviance: 7467
## Residual Deviance: 7222 AIC: 7230
##
## $frml2
##
## Call: glm(formula = x, family = binomial(link = "logit"), data = model_train)
##
## Coefficients:
## (Intercept) segmento02 - PARTICULARES
## -3.564e+00 -8.604e-01
## segmento03 - UNIVERSITARIO sexoV
## -2.286e+00 1.458e-01
## renta
## 1.208e-07
##
## Degrees of Freedom: 57874 Total (i.e. Null); 57870 Residual
## Null Deviance: 7467
## Residual Deviance: 7221 AIC: 7231שיטות ב Machine Learning
בחלק הזה אנחנו נסקור בקצרה 3 שיטות מובילות ב ML:
- Clustering
- Decision Trees
- PCA
את כל השיטות ניתן לעשות במגוון חבילות ובראשן CARET , שם יש פקודה אחת train שלוקחת כפרמטר את השיטה לחישוב
אנחנו ננסה להדגים עם כמה ספריות בשביל להרחיב קצת את האופקים….
מציאת אשכולות - Clustering
אפשר לחשב cluster ב R בעזרת ספריית cluster הבא מובנה עם התוכנה
על מנת לבדוק את מספר האשכולות נחשב את סכום השונות בתוך האשכולות לכמה
מספרי אשכולות ונבחר בהתאם
נשתמש בנתונים שכבר מרכזנו בפרק של טרנספורמציות…
הסינטקס הבסיסי:
kmeans(data,centers)
library(cluster)
###############################
# This is the basic function:
o_kmean<- kmeans(model_data_centered, centers=7)
############## Choose K through `Elbow` method #######
# Make a function that takes number of centers
# and returns the sum of within errors
c_func<-function(k){sum(kmeans(model_data_centered,centers=k)$withinss)}
c_func(8)## [1] 31123.61## Lets feed our function to `map_dbl`
# Iterate over h = numbers 1 to 30
h <- seq(1,30)
wss <- map_dbl(h,c_func)## Warning: Quick-TRANSfer stage steps exceeded maximum (= 4133900)
## Warning: Quick-TRANSfer stage steps exceeded maximum (= 4133900)## Warning: did not converge in 10 iterations
## Warning: did not converge in 10 iterations
## Warning: did not converge in 10 iterations
## Warning: did not converge in 10 iterations
## Warning: did not converge in 10 iterations
## Warning: did not converge in 10 iterations
## Warning: did not converge in 10 iterations# plot the results
plot(h,wss,type="b",xlab="N Clusters", ylab="Within SSE") 
######### Choose a cluster
# Based on the graph lets take k=7
o_kmean<- kmeans(model_data_centered, centers=7)
table(o_kmean$cluster)##
## 1 2 3 4 5 6 7
## 18428 18 2461 978 7419 27717 25657# Generic function to plot cluster
clusplot(model_data_centered, o_kmean$cluster, color=TRUE, shade=TRUE, lines=0)
- החישוב של
kmeansיכול להשתנות מהרצה להרצה - החבילה
pamמחשבתclusterשהאשכולות יותר עמידים בין ההרצות
עצי החלטה - Decision Trees
יש כמה ספריות שמריצים עצי החלטה.
אנחנו נדגים בעזרת הספרייה הבסיסית rpart
בנוסף לפורמולה אפשר לשחק עם הפרמטרים הספציפים.
כגון מדד החלטה לפי מה לפצל - פה בחרנו ב gini
פה הגדרנו ב method גם שהמשתנה הוא בינארי לעומת רציף או אורדינאלי
הסינטקס הבסיסי הוא:
rpart(formula, data, method, paramters)
library(rpart)
# 1. run the decision tree regression:
part_reg1 <- rpart( ind_cc ~ factor(last_3m_cc) + factor(last_6m_cc) + renta
,data = model_train
, method = "class"
,parms = list(prior = c(.98,.02), split = "gini"))
# 2. Show the outcomes
summary(part_reg1)## Call:
## rpart(formula = ind_cc ~ factor(last_3m_cc) + factor(last_6m_cc) +
## renta, data = model_train, method = "class", parms = list(prior = c(0.98,
## 0.02), split = "gini"))
## n= 57875
##
## CP nsplit rel error xerror xstd
## 1 0.02153478 0 1.0000000 1.0000000 0.03789736
## 2 0.01000000 3 0.9353957 0.9630882 0.03614006
##
## Variable importance
## factor(last_6m_cc) factor(last_3m_cc)
## 54 46
##
## Node number 1: 57875 observations, complexity param=0.02153478
## predicted class=0 expected loss=0.02 P(node) =1
## class counts: 57187 688
## probabilities: 0.980 0.020
## left son=2 (55702 obs) right son=3 (2173 obs)
## Primary splits:
## factor(last_6m_cc) splits as LRRRRRR, improve=488.60530, (190 missing)
## factor(last_3m_cc) splits as LRRR, improve=484.60000, (190 missing)
## renta < 100141.3 to the left, improve= 3.43347, (0 missing)
## Surrogate splits:
## factor(last_3m_cc) splits as LRRR, agree=0.99, adj=0.728, (0 split)
##
## Node number 2: 55702 observations
## predicted class=0 expected loss=0.006227078 P(node) =0.9569982
## class counts: 55497 205
## probabilities: 0.994 0.006
##
## Node number 3: 2173 observations, complexity param=0.02153478
## predicted class=0 expected loss=0.326514 P(node) =0.04300183
## class counts: 1690 483
## probabilities: 0.673 0.327
## left son=6 (1309 obs) right son=7 (864 obs)
## Primary splits:
## factor(last_3m_cc) splits as LLRR, improve=81.223990, (0 missing)
## factor(last_6m_cc) splits as -LLLRRR, improve=66.527150, (0 missing)
## renta < 152532.7 to the left, improve= 2.887691, (0 missing)
## Surrogate splits:
## factor(last_6m_cc) splits as -LLLRRR, agree=0.780, adj=0.447, (0 split)
## renta < 1058236 to the left, agree=0.603, adj=0.002, (0 split)
##
## Node number 6: 1309 observations
## predicted class=0 expected loss=0.2161122 P(node) =0.02461576
## class counts: 1126 183
## probabilities: 0.784 0.216
##
## Node number 7: 864 observations, complexity param=0.02153478
## predicted class=0 expected loss=0.4743228 P(node) =0.01838606
## class counts: 564 300
## probabilities: 0.526 0.474
## left son=14 (570 obs) right son=15 (294 obs)
## Primary splits:
## factor(last_3m_cc) splits as --LR, improve=18.191710, (0 missing)
## factor(last_6m_cc) splits as --LLLRR, improve=14.485360, (0 missing)
## renta < 305221.4 to the left, improve= 2.443772, (0 missing)
## Surrogate splits:
## factor(last_6m_cc) splits as --LLLLR, agree=0.823, adj=0.48, (0 split)
## renta < 478829.8 to the left, agree=0.663, adj=0.01, (0 split)
##
## Node number 14: 570 observations
## predicted class=0 expected loss=0.4045384 P(node) =0.01171303
## class counts: 407 163
## probabilities: 0.595 0.405
##
## Node number 15: 294 observations
## predicted class=1 expected loss=0.4031859 P(node) =0.00667303
## class counts: 157 137
## probabilities: 0.403 0.597part_reg1## n= 57875
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 57875 1157.5000 0 (0.980000000 0.020000000)
## 2) factor(last_6m_cc)=0 55702 344.8946 0 (0.993772922 0.006227078) *
## 3) factor(last_6m_cc)=1,2,3,4,5,6 2173 812.6054 0 (0.673486007 0.326513993)
## 6) factor(last_3m_cc)=0,1 1309 307.8815 0 (0.783887755 0.216112245) *
## 7) factor(last_3m_cc)=2,3 864 504.7238 0 (0.525677150 0.474322850)
## 14) factor(last_3m_cc)=2 570 274.2333 0 (0.595461601 0.404538399) *
## 15) factor(last_3m_cc)=3 294 155.7110 1 (0.403185916 0.596814084) *# 3. Plot the tree
par(c(1,1), xpd = NA)
plot(part_reg1)
text(part_reg1, use.n = TRUE)
# 4. Information plot:
plotcp(part_reg1)
text(part_reg1)
# 5. Create predictions and check the Confusion matrix:
rpartPred <- predict(part_reg1, model_data, type = "class")
# requires 2 factor vectors
confusionMatrix(rpartPred, model_data$ind_cc)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 81468 784
## 1 232 194
##
## Accuracy : 0.9877
## 95% CI : (0.9869, 0.9885)
## No Information Rate : 0.9882
## P-Value [Acc > NIR] : 0.8917
##
## Kappa : 0.2711
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.9972
## Specificity : 0.1984
## Pos Pred Value : 0.9905
## Neg Pred Value : 0.4554
## Prevalence : 0.9882
## Detection Rate : 0.9854
## Detection Prevalence : 0.9948
## Balanced Accuracy : 0.5978
##
## 'Positive' Class : 0
## ###################################################3
# With CARET train I wasnt able to get a result...
# But I'll leave you the code to try out:
# cvCtrl <- trainControl(method = "repeatedcv",
# repeats = 3,
# summaryFunction = twoClassSummary,
# classProbs = TRUE)
#
# set.seed(1)
#
# rpartTune <- train(fml1
# ,data = model_test
# , method = "rpart"
# ,tuneLength = 30
# ,metric = "Accuracy"
# ,trControl = cvCtrl)טרנספרומציה בעזרת Principal Components Analysis
ניתן להשתמש ב PCA כשיטה להקטנת מימדיות או למשל כשיש הרבה משתנים דומים …
אנחנו נשתמש באחד הפקודות של CARET לביצוע PCA כהכנה לשלב מידול הבא
נעשה זאת בעזרת הפקודה PreProcess
# choose only coloumns starting with `l6m`
small_data <- model_data %>% select(starts_with("l6m"))
# with CARET - preProcess
trans <- preProcess(small_data,
method=c("BoxCox", "center", "scale", "pca"))
PC <- predict(trans, small_data)## Warning in is.na(lam): is.na() applied to non-(list or vector) of type
## 'NULL'summary(PC)## PC1 PC2 PC3
## Min. :-9.6954 Min. :-10.8204 Min. :-10.6290
## 1st Qu.: 0.1722 1st Qu.: -0.5114 1st Qu.: -0.3887
## Median : 0.6315 Median : 0.3701 Median : 0.3370
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 1.0985 3rd Qu.: 0.5851 3rd Qu.: 0.4099
## Max. : 1.4682 Max. : 3.7232 Max. : 7.9477
## NA's :263 NA's :263 NA's :263
## PC4 PC5 PC6
## Min. :-13.8366 Min. :-21.5784 Min. :-33.60513
## 1st Qu.: -0.3009 1st Qu.: -0.1273 1st Qu.: 0.03613
## Median : 0.2332 Median : 0.1169 Median : 0.15040
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.2907 3rd Qu.: 0.2929 3rd Qu.: 0.29484
## Max. : 9.7139 Max. : 7.0058 Max. : 2.54533
## NA's :263 NA's :263 NA's :263
## PC7 PC8 PC9
## Min. :-26.0388 Min. :-61.67077 Min. :-217.65114
## 1st Qu.: -0.3244 1st Qu.: -0.09203 1st Qu.: -0.00005
## Median : -0.1025 Median : -0.06585 Median : 0.00031
## Mean : 0.0000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.2080 3rd Qu.: -0.02108 3rd Qu.: 0.05207
## Max. : 31.9290 Max. : 68.30237 Max. : 7.50764
## NA's :263 NA's :263 NA's :263
## PC10 PC11 PC12
## Min. :-26.34059 Min. :-5.75199 Min. :-21.76098
## 1st Qu.: -0.11193 1st Qu.:-0.04743 1st Qu.: -0.08614
## Median : 0.08484 Median : 0.08144 Median : -0.06919
## Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.13699 3rd Qu.: 0.09398 3rd Qu.: 0.12498
## Max. : 83.65843 Max. :60.16775 Max. : 20.53030
## NA's :263 NA's :263 NA's :263
## PC13 PC14 PC15
## Min. :-14.66844 Min. :-37.13095 Min. :-13.8914
## 1st Qu.: -0.07799 1st Qu.: 0.00593 1st Qu.: -0.2269
## Median : 0.10240 Median : 0.03172 Median : 0.1578
## Mean : 0.00000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.13063 3rd Qu.: 0.04764 3rd Qu.: 0.1841
## Max. : 15.09130 Max. : 7.54347 Max. : 24.8700
## NA's :263 NA's :263 NA's :263
## PC16 PC17 PC18
## Min. :-7.47574 Min. :-10.0090 Min. :-5.59104
## 1st Qu.:-0.16543 1st Qu.: -0.3184 1st Qu.:-0.31505
## Median : 0.05723 Median : -0.2348 Median : 0.08597
## Mean : 0.00000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.23509 3rd Qu.: 0.4395 3rd Qu.: 0.14893
## Max. :10.39700 Max. : 18.0026 Max. : 7.55454
## NA's :263 NA's :263 NA's :263
## PC19 PC20
## Min. :-4.6132 Min. :-3.93936
## 1st Qu.:-0.3201 1st Qu.:-0.23212
## Median :-0.2454 Median : 0.00295
## Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.4416 3rd Qu.: 0.45422
## Max. : 4.3491 Max. : 3.70967
## NA's :263 NA's :263חבילות מומלצות
מה שהוצג במערך הזה הוא רק בגדר ספתח ל ML ב R
מי שרוצה להיכנס לנושא לעומק צריך ללכלך את הידיים ולהבין את הפרמטרים של כל
פקודה ושיטה.
מספר חבילות מומלצות ל ML ב R:
- CARET - להבין את מגוון האפשרויות שלה - כעקרון יש פקודה בסיסית של
trainשעליה מזינים את השיטה שרוצים להריץ ואת הפרמטרים הרלוונטיים.
דרוש ניסיון והתנסות בשביל לשלוט ולהתאים כל שיטה…
H20 - אחד החבילות המובילות בתחום ה ML - החבילה מייצרת מחשב ורטואלי בתוך ה R ובעצם מפעילה את האלגורתמים בחישובים מקבילים כש R זה רק מעטפת
randomForest - יערות החלטה
XGBOOST - חישוב של יערות החלטה מוגבר… בערך
nnet ו neuralnet - ספריות לחישוב neural network המדמה את הרצפטורים של המוח - כביכול…
ועוד …….