Statistical Learning - Group 16 - Final Presentation
Angelina Khatiwada & Rijin Baby
18/06/2021
Data Set: in-vehicle coupon recommendation
This data was collected via a survey on Amazon Mechanical Turk. The survey describes different driving scenarios including the destination, current time, weather, passenger, etc., and then ask the person whether he will accept the coupon if he is the driver.
- Multivariate dataset
- Number of Instances: 12684
- Number of Attributes: 23
Key Attribute Information:
- destination: Home, No Urgent Place, Work
- passanger: Alone Friend(s), Kid(s), Partner
- coupon: Bar, Carry out & Take away, Coffee House, Restaurant(<20), Restaurant(20-50)
- expiration: 1d, 2h (the coupon expires in 1 day or in 2 hours)
- Bar: never, less1, 1~3, gt8, nan4~8 (feature meaning: how many times do you go to a bar every month?)
- CoffeeHouse: never, less1, 4~8, 1~3, gt8, nan (feature meaning: how many times do you go to a coffeehouse every month?)
- CarryAway:n4~8, 1~3, gt8, less1, never (feature meaning: how many times do you get take-away food every month?)
- RestaurantLessThan20: 4~8, 1~3, less1, gt8, never (feature meaning: how many times do you go to a restaurant with an average expense per person of less than $20 every month?)
- toCoupon_GEQ15min:0,1 (feature meaning: driving distance to the restaurant/bar for using the coupon is greater than 15 minutes)
- direction_same:0, 1 (feature meaning: whether the restaurant/bar is in the same direction as your current destination)
- Y:1, 0 (whether the coupon is accepted)
Data Preparation
1) Load Libraries & Read Dataset
2) Missing & unique value check
## Columns with missing value## [1] "car :100%" "Bar :0.84%"
## [3] "CoffeeHouse :1.71%" "CarryAway :1.19%"
## [5] "RestaurantLessThan20 :1.02%" "Restaurant20To50 :1.49%"## Column with unique value## toCoupon_GEQ5min
## 203) Creating New Variables - occupation class - Expiration weightage (numeric, scaled)
| occupation_class | Actual_occupation |
|---|---|
| Craft and related trades workers | Installation Maintenance & Repair |
| Craft and related trades workers | Transportation & Material Moving |
| Craft and related trades workers | Food Preparation & Serving Related |
| Craft and related trades workers | Building & Grounds Cleaning & Maintenance |
| Others | Office & Administrative Support |
| Others | Production Occupations |
| Others | Farming Fishing & Forestry |
| Professionals | Architecture & Engineering |
| Professionals | Education&Training&Library |
| Professionals | Healthcare Practitioners & Technical |
| Professionals | Management |
| Professionals | Arts Design Entertainment Sports & Media |
| Professionals | Computer & Mathematical |
| Professionals | Legal |
| Professionals | Business & Financial |
| Retired | Retired |
| Service and sales | Sales & Related |
| Service and sales | Personal Care & Service |
| Service and sales | Protective Service |
| Student | Student |
| Technicians & prof | Healthcare Support |
| Technicians & prof | Life Physical Social Science |
| Technicians & prof | Community & Social Services |
| Technicians & prof | Construction & Extraction |
| Unemployed | Unemployed |
4) Missing imputation knn approach
library(VIM)
cleaned_data <- kNN(
coupon_data,
variable = c("Bar","CoffeeHouse","CarryAway","RestaurantLessThan20","Restaurant20To50")
, k = 5)
cleaned_data <- cleaned_data[,1:ncol(coupon_data)]
# coupon_data_final %>% map(table)
# colMeans(is.na(cleaned_data))*100Plotting data columns
Modeling
- Droped the variable ‘direction_opp’ that is perfectly correlated with the varible ‘direction_same’.
- ‘occupation’ will be replaced by ‘occupation_class’ in the model.
- ‘expiration’ replaced by numeric variable.
Dataset Summary
## 'data.frame': 12684 obs. of 22 variables:
## $ destination : chr "No Urgent Place" "No Urgent Place" "No Urgent Place" "No Urgent Place" ...
## $ passanger : chr "Alone" "Friends" "Friends" "Friends" ...
## $ weather : chr "Sunny" "Sunny" "Sunny" "Sunny" ...
## $ temperature : chr "55" "80" "80" "80" ...
## $ coupon : chr "Restaurant(<20)" "Coffee House" "Carry out & Take away" "Coffee House" ...
## $ gender : chr "Female" "Female" "Female" "Female" ...
## $ age : chr "21" "21" "21" "21" ...
## $ maritalStatus : chr "Unmarried partner" "Unmarried partner" "Unmarried partner" "Unmarried partner" ...
## $ has_children : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
## $ education : chr "Some college - no degree" "Some college - no degree" "Some college - no degree" "Some college - no degree" ...
## $ income : chr "$37500 - $49999" "$37500 - $49999" "$37500 - $49999" "$37500 - $49999" ...
## $ Bar : chr "never" "never" "never" "never" ...
## $ CoffeeHouse : chr "never" "never" "never" "never" ...
## $ CarryAway : chr "1~3" "less1" "less1" "less1" ...
## $ RestaurantLessThan20: chr "4~8" "4~8" "4~8" "4~8" ...
## $ Restaurant20To50 : chr "1~3" "1~3" "1~3" "1~3" ...
## $ toCoupon_GEQ15min : Factor w/ 2 levels "0","1": 1 1 2 2 2 2 2 2 2 2 ...
## $ toCoupon_GEQ25min : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ direction_same : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ Y : Factor w/ 2 levels "0","1": 2 1 2 1 1 2 2 2 2 1 ...
## $ occupation_class : chr "Unemployed" "Unemployed" "Unemployed" "Unemployed" ...
## $ expiration_weightage: num 1.334 0.111 0.111 0.111 1.334 ...
Logistic Regression
The multiple logistic regression is used to predict the probability of class membership based on multiple predictor variables.
Run model with all parameters on dataset without splitting
- The most significant variables (highest z-score and lowest p-value): destination, weather, coupon, CoffeHouse, direction, expiration_weightage
- Insignificant variables (p-value > 0.05): has_children
The coefficient estimate of the variable destinationNo Urgent Place is b = 0.72: means that an No UrgentPlace is associated with increase in the probability of accepting coupon more than baseline dummy destinationWork.
Being Male, for example, increases the probability of accepting coupon more than being Female. Being age50plus decreases the probability of accepting coupon more than being 21-26.
If we have at least one significant dummy, we can not exclude the variable.
set.seed(123)
full_log_model <- glm(Y ~., data = cleaned_data, family=binomial(link='logit'))
summary(full_log_model)##
## Call:
## glm(formula = Y ~ ., family = binomial(link = "logit"), data = cleaned_data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5150 -1.0280 0.5651 0.9452 2.4271
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.842614 0.243444 -7.569 3.76e-14 ***
## destinationNo Urgent Place 0.728894 0.079446 9.175 < 2e-16 ***
## destinationWork -0.097692 0.057683 -1.694 0.090341 .
## passangerFriends 0.198204 0.072183 2.746 0.006036 **
## passangerKids -0.373699 0.093329 -4.004 6.23e-05 ***
## passangerPartner 0.244527 0.087174 2.805 0.005031 **
## weatherSnowy 0.146718 0.131864 1.113 0.265859
## weatherSunny 0.572450 0.079709 7.182 6.88e-13 ***
## temperature55 0.350954 0.095173 3.688 0.000226 ***
## temperature80 0.186184 0.085180 2.186 0.028833 *
## couponCarry out & Take away 1.718958 0.072760 23.625 < 2e-16 ***
## couponCoffee House 0.540722 0.064185 8.424 < 2e-16 ***
## couponRestaurant(<20) 1.522600 0.070056 21.734 < 2e-16 ***
## couponRestaurant(20-50) 0.418126 0.079787 5.241 1.60e-07 ***
## genderMale 0.219835 0.043048 5.107 3.28e-07 ***
## age26 0.004076 0.069442 0.059 0.953196
## age31 -0.158493 0.075461 -2.100 0.035700 *
## age36 -0.075832 0.086472 -0.877 0.380513
## age41 0.025218 0.096172 0.262 0.793151
## age46 0.013301 0.105949 0.126 0.900096
## age50plus -0.213604 0.086328 -2.474 0.013348 *
## agebelow21 0.003143 0.119334 0.026 0.978987
## maritalStatusMarried partner 0.154955 0.112572 1.376 0.168667
## maritalStatusSingle 0.305867 0.117405 2.605 0.009181 **
## maritalStatusUnmarried partner 0.110064 0.120425 0.914 0.360736
## maritalStatusWidowed 0.261376 0.230633 1.133 0.257088
## has_children1 0.092140 0.058092 1.586 0.112716
## educationBachelors degree -0.157473 0.077792 -2.024 0.042943 *
## educationGraduate degree -0.323460 0.090673 -3.567 0.000361 ***
## educationHigh School Graduate 0.174644 0.104765 1.667 0.095513 .
## educationSome college - no degree 0.048642 0.078520 0.619 0.535596
## educationSome High School 0.654488 0.285999 2.288 0.022113 *
## income$12500 - $24999 -0.093274 0.084546 -1.103 0.269926
## income$25000 - $37499 0.045049 0.079759 0.565 0.572199
## income$37500 - $49999 -0.081985 0.081278 -1.009 0.313119
## income$50000 - $62499 0.162556 0.081009 2.007 0.044788 *
## income$62500 - $74999 -0.322062 0.096822 -3.326 0.000880 ***
## income$75000 - $87499 -0.232408 0.095601 -2.431 0.015056 *
## income$87500 - $99999 -0.275381 0.096396 -2.857 0.004280 **
## incomeLess than $12500 -0.108887 0.099748 -1.092 0.275000
## Bar4~8 -0.112211 0.086349 -1.300 0.193770
## Bargt8 -0.416026 0.143365 -2.902 0.003710 **
## Barless1 -0.198242 0.063643 -3.115 0.001840 **
## Barnever -0.181372 0.061517 -2.948 0.003195 **
## CoffeeHouse4~8 -0.032526 0.070143 -0.464 0.642854
## CoffeeHousegt8 -0.388980 0.085206 -4.565 4.99e-06 ***
## CoffeeHouseless1 -0.447871 0.057515 -7.787 6.86e-15 ***
## CoffeeHousenever -0.943201 0.063176 -14.930 < 2e-16 ***
## CarryAway4~8 -0.016472 0.050861 -0.324 0.746042
## CarryAwaygt8 -0.071545 0.074833 -0.956 0.339039
## CarryAwayless1 -0.141711 0.063818 -2.221 0.026382 *
## CarryAwaynever 0.018829 0.188179 0.100 0.920298
## RestaurantLessThan204~8 0.060047 0.052078 1.153 0.248901
## RestaurantLessThan20gt8 0.182197 0.085667 2.127 0.033437 *
## RestaurantLessThan20less1 0.092708 0.062651 1.480 0.138940
## RestaurantLessThan20never 0.289231 0.164846 1.755 0.079336 .
## Restaurant20To504~8 0.071296 0.099458 0.717 0.473466
## Restaurant20To50gt8 0.075333 0.178821 0.421 0.673551
## Restaurant20To50less1 -0.127060 0.050898 -2.496 0.012547 *
## Restaurant20To50never -0.300305 0.068515 -4.383 1.17e-05 ***
## toCoupon_GEQ15min1 -0.101592 0.045957 -2.211 0.027064 *
## toCoupon_GEQ25min1 0.183414 0.082056 2.235 0.025402 *
## direction_same1 0.536336 0.069056 7.767 8.06e-15 ***
## occupation_classOthers 0.244455 0.122058 2.003 0.045202 *
## occupation_classProfessionals 0.076542 0.098637 0.776 0.437750
## occupation_classRetired -0.073831 0.143914 -0.513 0.607939
## occupation_classService and sales 0.140400 0.109172 1.286 0.198429
## occupation_classStudent 0.054337 0.115699 0.470 0.638611
## occupation_classTechnicians & prof 0.333932 0.122577 2.724 0.006445 **
## occupation_classUnemployed 0.037047 0.106977 0.346 0.729112
## expiration_weightage 0.684226 0.035094 19.497 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 17345 on 12683 degrees of freedom
## Residual deviance: 14991 on 12613 degrees of freedom
## AIC: 15133
##
## Number of Fisher Scoring iterations: 3An odds ratio measures the association between a x and y. It represents the ratio of the odds that an event will occur given the presence of the predictor x, compared to the odds of the event occurring in the absence of that predictor. The associated beta coefficient (b) in the logistic regression function corresponds to the log of the odds ratio for that predictor. If the odds ratio is 0.5, then the odds that the event = 1 are 0.5 higher when the predictor x is present (x = 1) versus x is absent (x = 0).
Chi-square test to check the overall effect of variables on the dependent variable
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 12683 17345
## destination 2 217.51 12681 17128 < 2.2e-16 ***
## passanger 3 77.17 12678 17051 < 2.2e-16 ***
## weather 2 104.77 12676 16946 < 2.2e-16 ***
## temperature 2 6.24 12674 16940 0.044132 *
## coupon 4 897.31 12670 16042 < 2.2e-16 ***
## gender 1 25.80 12669 16017 3.784e-07 ***
## age 7 52.77 12662 15964 4.109e-09 ***
## maritalStatus 4 10.99 12658 15953 0.026729 *
## has_children 1 1.39 12657 15952 0.239020
## education 5 35.14 12652 15916 1.414e-06 ***
## income 8 39.17 12644 15877 4.579e-06 ***
## Bar 4 61.01 12640 15816 1.777e-12 ***
## CoffeeHouse 4 296.51 12636 15520 < 2.2e-16 ***
## CarryAway 4 5.34 12632 15514 0.254013
## RestaurantLessThan20 4 10.39 12628 15504 0.034317 *
## Restaurant20To50 4 24.51 12624 15479 6.315e-05 ***
## toCoupon_GEQ15min 1 38.38 12623 15441 5.833e-10 ***
## toCoupon_GEQ25min 1 3.11 12622 15438 0.077887 .
## direction_same 1 33.88 12621 15404 5.868e-09 ***
## occupation_class 7 19.09 12614 15385 0.007896 **
## expiration_weightage 1 393.52 12613 14991 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Insignificant coefficients (with p-value > 0.05) according to chi2 check are:
- has_children
- CarryAway
- toCoupon_GEQ25min
Multicollinearity check
car::vif(full_log_model)## GVIF Df GVIF^(1/(2*Df))
## destination 4.226976 2 1.433862
## passanger 3.711907 3 1.244322
## weather 3.718154 2 1.388615
## temperature 4.020295 2 1.416004
## coupon 1.663346 4 1.065670
## gender 1.189669 1 1.090719
## age 5.084769 7 1.123176
## maritalStatus 3.363088 4 1.163703
## has_children 2.110278 1 1.452680
## education 2.228544 5 1.083433
## income 3.173868 8 1.074854
## Bar 2.410418 4 1.116250
## CoffeeHouse 2.160243 4 1.101065
## CarryAway 1.996979 4 1.090302
## RestaurantLessThan20 2.686835 4 1.131501
## Restaurant20To50 2.525245 4 1.122763
## toCoupon_GEQ15min 1.335446 1 1.155615
## toCoupon_GEQ25min 1.861952 1 1.364534
## direction_same 2.114196 1 1.454028
## occupation_class 4.523701 7 1.113835
## expiration_weightage 1.176458 1 1.084646There is no problem of mulricollinearity.
Keeping insignificant variables in the model may contribute to overfitting. Eliminated manually, using stepwise regression and penalized regression methods.
Stepwise Selection
Both directions
## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## Y ~ destination + passanger + weather + temperature + coupon +
## gender + age + maritalStatus + has_children + education +
## income + Bar + CoffeeHouse + CarryAway + RestaurantLessThan20 +
## Restaurant20To50 + toCoupon_GEQ15min + toCoupon_GEQ25min +
## direction_same + occupation_class + expiration_weightage
##
## Final Model:
## Y ~ destination + passanger + weather + temperature + coupon +
## gender + age + maritalStatus + has_children + education +
## income + Bar + CoffeeHouse + Restaurant20To50 + toCoupon_GEQ15min +
## toCoupon_GEQ25min + direction_same + occupation_class + expiration_weightage
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 12613 14991.42 15133.42
## 2 - CarryAway 4 5.724530 12617 14997.14 15131.14
## 3 - RestaurantLessThan20 4 7.139874 12621 15004.28 15130.28Forward direction
## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## Y ~ destination + passanger + weather + temperature + coupon +
## gender + age + maritalStatus + has_children + education +
## income + Bar + CoffeeHouse + CarryAway + RestaurantLessThan20 +
## Restaurant20To50 + toCoupon_GEQ15min + toCoupon_GEQ25min +
## direction_same + occupation_class + expiration_weightage
##
## Final Model:
## Y ~ destination + passanger + weather + temperature + coupon +
## gender + age + maritalStatus + has_children + education +
## income + Bar + CoffeeHouse + CarryAway + RestaurantLessThan20 +
## Restaurant20To50 + toCoupon_GEQ15min + toCoupon_GEQ25min +
## direction_same + occupation_class + expiration_weightage
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 12613 14991.42 15133.42With forward direction we get the same result as we had with a complete model, forward begins with a Null model (intercept only model)
Backward direction
## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## Y ~ destination + passanger + weather + temperature + coupon +
## gender + age + maritalStatus + has_children + education +
## income + Bar + CoffeeHouse + CarryAway + RestaurantLessThan20 +
## Restaurant20To50 + toCoupon_GEQ15min + toCoupon_GEQ25min +
## direction_same + occupation_class + expiration_weightage
##
## Final Model:
## Y ~ destination + passanger + weather + temperature + coupon +
## gender + age + maritalStatus + has_children + education +
## income + Bar + CoffeeHouse + Restaurant20To50 + toCoupon_GEQ15min +
## toCoupon_GEQ25min + direction_same + occupation_class + expiration_weightage
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 12613 14991.42 15133.42
## 2 - CarryAway 4 5.724530 12617 14997.14 15131.14
## 3 - RestaurantLessThan20 4 7.139874 12621 15004.28 15130.28The backward procedure eliminated exactly the same variables as the “both” procedure.
Different criteria can be assigned to the stepAIC() function for stepwise selection. The default is AIC, which is performed by assigning the argument k to 2.
We also tried running bestglm, but there is a hard-coded constraint in bestglm (15 predictors means there are 2^15 = 32768 candidate models). In our case we have around 70 variables counting dummies.
Running 10 fold Cross-Validation on Complete and Stepwise model using training data
Initially splitting data into train (0.67) and test sets to avoid overfitting
1) Model on training set with all variables
## Generalized Linear Model
##
## 8499 samples
## 21 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 7650, 7649, 7649, 7650, 7649, 7648, ...
## Resampling results:
##
## Accuracy Kappa
## 0.6862035 0.35110172) Model on training set by dropping columns suggested by stepwise model (CarryAway, RestaurantLessThan20)
## Generalized Linear Model
##
## 8499 samples
## 19 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 7649, 7649, 7649, 7649, 7649, 7650, ...
## Resampling results:
##
## Accuracy Kappa
## 0.6887878 0.3564379Full model AIC: 10139
Stepwise model AIC: 10134
According to the bias-variance trade-off, all things equal, simpler model should be always preferred because it is less likely to overfit the training data (less variance).
Logistic model Predictions on Test set
Logistic Complete Model Accuracy
log_reg_prob1 <- predict(log_reg_full, test, type = 'prob')
log_reg_pred1 <- ifelse(log_reg_prob1[2] > 0.5, 1, 0)
mean(log_reg_pred1 == test$Y)## [1] 0.6824373Stepwise Model Accuracy
## [1] 0.6862605Stepwise model Confusion matrix
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1039 546
## 1 767 1833
##
## Accuracy : 0.6863
## 95% CI : (0.672, 0.7003)
## No Information Rate : 0.5685
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.351
##
## Mcnemar's Test P-Value : 1.268e-09
##
## Sensitivity : 0.7705
## Specificity : 0.5753
## Pos Pred Value : 0.7050
## Neg Pred Value : 0.6555
## Prevalence : 0.5685
## Detection Rate : 0.4380
## Detection Prevalence : 0.6213
## Balanced Accuracy : 0.6729
##
## 'Positive' Class : 1
## High sensitivity: fewer False Negative errors, which are severe types of error for us.
Low specificity: Many False Positive.
Sensitivity = TP/TP+FN
The ability of a test to correctly identify people with a coupon
Speficity = TN/TP+FN
We can change the threshold from 0.5 to 0.45 in order to decrease FN. Focus on sensitivity which is more important for business opportunities.
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 884 415
## 1 922 1964
##
## Accuracy : 0.6805
## 95% CI : (0.6662, 0.6946)
## No Information Rate : 0.5685
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.3261
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.8256
## Specificity : 0.4895
## Pos Pred Value : 0.6805
## Neg Pred Value : 0.6805
## Prevalence : 0.5685
## Detection Rate : 0.4693
## Detection Prevalence : 0.6896
## Balanced Accuracy : 0.6575
##
## 'Positive' Class : 1
## ROC curve
test_roc = roc(test$Y ~ log_reg_prob1$"1", plot = TRUE, print.auc = TRUE)
Penalised models: Lasso and Elastic Net
With dummy variables we have 70 variables, some of the coefficients are already shrinked. We tried to apply Lasso and Elastic Net to find a reduced set of variables resulting to an optimal performing mode (penalty for having too many variables - regularization).
Lasso regression: the coefficients of some less contributive variables are zero. Only the most significant variables are kept in the final model.
Elastic net regression: the combination of ridge and lasso regression. It shrinks some coefficients toward zero and set some coefficients to exactly zero.
set.seed(123)
x=model.matrix(Y~., data=train)[,-20]
y=as.numeric(train$Y)
x.test <- model.matrix(Y ~., test)[,-20]
fit.lasso= glmnet(x, y, family = "binomial", alpha = 1) #lambda = NULL
plot(fit.lasso,xvar="lambda",label=TRUE)
The log of the optimal value of lambda is around -6, which is the one that minimizes the prediction error. The purpose of regularization is to balance accuracy and simplicity.
lasso_model <- glmnet(x, y, family = "binomial", alpha = 1, lambda = cv.lasso$lambda.min)
# Display regression coefficients
print('The value of the optimal lambda:')## [1] "The value of the optimal lambda:"cv.lasso$lambda.min## [1] 0.001394676coef(lasso_model)## 71 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) -1.395440651
## (Intercept) .
## destinationNo Urgent Place 0.687677416
## destinationWork -0.029055342
## passangerFriends 0.161064848
## passangerKids -0.346999778
## passangerPartner 0.205737684
## weatherSnowy .
## weatherSunny 0.494608065
## temperature55 0.231844081
## temperature80 0.107032017
## couponCarry out & Take away 1.658855726
## couponCoffee House 0.463145232
## couponRestaurant(<20) 1.446256175
## couponRestaurant(20-50) 0.311719740
## genderMale 0.207668395
## age26 0.015276996
## age31 -0.157707910
## age36 -0.069876006
## age41 .
## age50plus -0.172175907
## agebelow21 .
## maritalStatusMarried partner .
## maritalStatusSingle 0.059809522
## maritalStatusUnmarried partner -0.096114244
## maritalStatusWidowed -0.053111199
## has_children1 0.007920614
## educationBachelors degree -0.107266369
## educationGraduate degree -0.288333043
## educationHigh School Graduate 0.215209216
## educationSome college - no degree 0.067963413
## educationSome High School 0.527705570
## income$12500 - $24999 -0.007796724
## income$25000 - $37499 0.087179080
## income$37500 - $49999 .
## income$50000 - $62499 0.182097445
## income$62500 - $74999 -0.273247786
## income$75000 - $87499 -0.151325805
## income$87500 - $99999 -0.180826363
## incomeLess than $12500 -0.076945301
## Bar4~8 .
## Bargt8 -0.201588609
## Barless1 -0.110958718
## Barnever -0.134647931
## CoffeeHouse4~8 0.007775047
## CoffeeHousegt8 -0.304124375
## CoffeeHouseless1 -0.412770446
## CoffeeHousenever -0.903520517
## CarryAway4~8 -0.025069505
## CarryAwaygt8 -0.136041211
## CarryAwayless1 -0.107593147
## CarryAwaynever -0.134109570
## RestaurantLessThan204~8 0.063105560
## RestaurantLessThan20gt8 0.090030765
## RestaurantLessThan20less1 0.017064766
## RestaurantLessThan20never 0.281946770
## Restaurant20To504~8 0.117941942
## Restaurant20To50gt8 0.103178598
## Restaurant20To50less1 -0.170603732
## Restaurant20To50never -0.270385703
## toCoupon_GEQ15min1 -0.074680603
## toCoupon_GEQ25min1 0.041470774
## direction_same1 0.478435414
## occupation_classOthers 0.126544850
## occupation_classProfessionals 0.028604913
## occupation_classRetired -0.068766314
## occupation_classService and sales .
## occupation_classStudent -0.019337819
## occupation_classTechnicians & prof 0.156277125
## occupation_classUnemployed -0.010457848
## expiration_weightage 0.637153176#Make predictions on the test data
probabilities <- lasso_model %>% predict(x.test)
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
#Model accuracy
observed.classes <- test$Y
mean(predicted.classes == observed.classes)## [1] 0.6551971Results do not show any improvement in the the model performance on the test data.
ELASTIC NET
For elastic net regression, alpha is between 0 and 1.
We automatically select the best tuning parameters alpha and lambda
# Build the model
set.seed(123)
elastic <- train(
Y ~., data = train, method = "glmnet",
family = "binomial",
trControl = trainControl("cv", number = 10),
tuneLength = 5
)
# Best tuning parameter
elastic$bestTune## alpha lambda
## 13 0.55 0.001672172Training model with the best parameters and predicting on the test set
elastic_model <- glmnet(x, y, family = "binomial", alpha = elastic$bestTune$alpha, lambda = elastic$bestTune$lambda) #lambda = NULL
# Display regression coefficients
#coef(lasso_model)
#Make predictions on the test data
probabilities <- elastic_model %>% predict(x.test)
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
#Model accuracy
observed.classes <- test$Y
mean(predicted.classes == observed.classes)## [1] 0.6568698Elastic Net performs a bit better than Lasso.
Linear Discriminant Analysis
lda.fit=lda(Y~.,data = train)
lda.fit## Call:
## lda(Y ~ ., data = train)
##
## Prior probabilities of groups:
## 0 1
## 0.4315802 0.5684198
##
## Group means:
## destinationNo Urgent Place destinationWork passangerFriends passangerKids
## 0 0.4209378 0.2870774 0.1984733 0.08969466
## 1 0.5549576 0.2196233 0.3044918 0.07120679
## passangerPartner weatherSnowy weatherSunny temperature55 temperature80
## 0 0.08069793 0.13331516 0.7489095 0.3290622 0.4738277
## 1 0.09314842 0.09232043 0.8298489 0.2862761 0.5435728
## couponCarry out & Take away couponCoffee House couponRestaurant(<20)
## 0 0.1090513 0.3710469 0.1496728
## 1 0.2411509 0.2777893 0.2769613
## couponRestaurant(20-50) genderMale age26 age31 age36 age41
## 0 0.15130862 0.4612868 0.1870229 0.1723010 0.11068702 0.08560523
## 1 0.09169944 0.5056924 0.2119644 0.1506934 0.09728835 0.08818050
## age46 age50plus agebelow21 maritalStatusMarried partner
## 0 0.05261723 0.1589422 0.03489640 0.4201200
## 1 0.05237011 0.1275098 0.04595322 0.3877044
## maritalStatusSingle maritalStatusUnmarried partner maritalStatusWidowed
## 0 0.3454198 0.1761178 0.015267176
## 1 0.3930863 0.1718071 0.008693852
## has_children1 educationBachelors degree educationGraduate degree
## 0 0.4375682 0.3546892 0.1559433
## 1 0.3968123 0.3361623 0.1324777
## educationHigh School Graduate educationSome college - no degree
## 0 0.06733915 0.3241549
## 1 0.07555372 0.3595529
## educationSome High School income$12500 - $24999 income$25000 - $37499
## 0 0.005179935 0.1420393 0.1488550
## 1 0.009107845 0.1440696 0.1629062
## income$37500 - $49999 income$50000 - $62499 income$62500 - $74999
## 0 0.1442203 0.1251363 0.07279171
## 1 0.1446905 0.1397226 0.06085697
## income$75000 - $87499 income$87500 - $99999 incomeLess than $12500 Bar4~8
## 0 0.08015267 0.07360960 0.08042530 0.07142857
## 1 0.05733803 0.06706686 0.08528255 0.09356241
## Bargt8 Barless1 Barnever CoffeeHouse4~8 CoffeeHousegt8 CoffeeHouseless1
## 0 0.02671756 0.2824427 0.4449291 0.1185932 0.08315158 0.2876227
## 1 0.02960050 0.2829642 0.3792176 0.1575243 0.09066446 0.2651625
## CoffeeHousenever CarryAway4~8 CarryAwaygt8 CarryAwayless1 CarryAwaynever
## 0 0.3009815 0.3369684 0.1286805 0.1682116 0.012540894
## 1 0.1900228 0.3498241 0.1270958 0.1330987 0.009107845
## RestaurantLessThan204~8 RestaurantLessThan20gt8 RestaurantLessThan20less1
## 0 0.271265 0.09351145 0.1782988
## 1 0.294349 0.10908715 0.1535914
## RestaurantLessThan20never Restaurant20To504~8 Restaurant20To50gt8
## 0 0.01826609 0.04280262 0.01635769
## 1 0.01676671 0.06603188 0.02649555
## Restaurant20To50less1 Restaurant20To50never toCoupon_GEQ15min1
## 0 0.5016358 0.1949291 0.6090513
## 1 0.4713310 0.1598013 0.5261851
## toCoupon_GEQ25min1 direction_same1 occupation_classOthers
## 0 0.16248637 0.2003817 0.05806979
## 1 0.09045746 0.2171393 0.06582488
## occupation_classProfessionals occupation_classRetired
## 0 0.3977644 0.04770992
## 1 0.3914303 0.03311944
## occupation_classService and sales occupation_classStudent
## 0 0.1164122 0.1115049
## 1 0.1128131 0.1314428
## occupation_classTechnicians & prof occupation_classUnemployed
## 0 0.05916031 0.1567612
## 1 0.06375492 0.1457255
## expiration_weightage
## 0 0.7111097
## 1 0.8655550
##
## Coefficients of linear discriminants:
## LD1
## destinationNo Urgent Place 0.813177729
## destinationWork -0.070080891
## passangerFriends 0.208370081
## passangerKids -0.403107630
## passangerPartner 0.321473973
## weatherSnowy 0.187316197
## weatherSunny 0.607018816
## temperature55 0.400603495
## temperature80 0.245265120
## couponCarry out & Take away 1.994767404
## couponCoffee House 0.637269307
## couponRestaurant(<20) 1.735688877
## couponRestaurant(20-50) 0.494646962
## genderMale 0.236254826
## age26 0.008815704
## age31 -0.214391396
## age36 -0.121976162
## age41 -0.022040646
## age46 0.030075786
## age50plus -0.219693586
## agebelow21 -0.002374899
## maritalStatusMarried partner 0.064875722
## maritalStatusSingle 0.168788119
## maritalStatusUnmarried partner -0.038418038
## maritalStatusWidowed -0.033118049
## has_children1 0.085937730
## educationBachelors degree -0.162979134
## educationGraduate degree -0.361950647
## educationHigh School Graduate 0.285132939
## educationSome college - no degree 0.069865151
## educationSome High School 0.537923590
## income$12500 - $24999 -0.077022823
## income$25000 - $37499 0.050644130
## income$37500 - $49999 -0.072436049
## income$50000 - $62499 0.171546333
## income$62500 - $74999 -0.384698313
## income$75000 - $87499 -0.231434583
## income$87500 - $99999 -0.297018677
## incomeLess than $12500 -0.155226888
## Bar4~8 -0.088871251
## Bargt8 -0.370006357
## Barless1 -0.190356486
## Barnever -0.206910019
## CoffeeHouse4~8 0.000917293
## CoffeeHousegt8 -0.378030794
## CoffeeHouseless1 -0.483216073
## CoffeeHousenever -1.047960002
## CarryAway4~8 -0.064923378
## CarryAwaygt8 -0.196022289
## CarryAwayless1 -0.151099339
## CarryAwaynever -0.265174018
## RestaurantLessThan204~8 0.108687206
## RestaurantLessThan20gt8 0.169837647
## RestaurantLessThan20less1 0.077156231
## RestaurantLessThan20never 0.395936973
## Restaurant20To504~8 0.103955408
## Restaurant20To50gt8 0.163389659
## Restaurant20To50less1 -0.188367799
## Restaurant20To50never -0.326505292
## toCoupon_GEQ15min1 -0.117267708
## toCoupon_GEQ25min1 0.154603029
## direction_same1 0.635661820
## occupation_classOthers 0.264532648
## occupation_classProfessionals 0.115526192
## occupation_classRetired -0.021122386
## occupation_classService and sales 0.104580884
## occupation_classStudent 0.019922008
## occupation_classTechnicians & prof 0.281574402
## occupation_classUnemployed 0.039960347
## expiration_weightage 0.729621403plot(lda.fit)
lda.pred=predict(lda.fit,test)$class
table(lda.pred,test$Y)##
## lda.pred 0 1
## 0 1037 554
## 1 769 1825mean(lda.pred==test$Y)## [1] 0.683871LDA predicts more False Positives as well, accuracy is similar to Logistic Regression (Robust result).
Improving the accuracy: - Remove potential outliers (Robust Logistic Regression, remove outliers within categorical variables and combinations of categorical variables) - Random forest and boosting models to capture non-linearities and interactions between the variables - Another option to apply Generalised Additive Models.
RANDOM FOREST
- Non-linear classifier, made up of an ensemble of decision trees that are independent of each other and each will predict the outcome variable using their own set of rules
- Random Forest Classifier is more of Accuracy focused algorithm and is best till it’s used with proper fit, else it gets overfit quickly
- Random Forest Classifier performs better with more categorical data than numeric and logistic regression is a little confusing when comes to categorical data
- RFC can also tell us how much each feature contributes to class prediction using Feature Importance and tree graphs for better interpretation.
Dataset Summary
## 'data.frame': 12684 obs. of 22 variables:
## $ destination : Factor w/ 3 levels "Home","No Urgent Place",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ passanger : Factor w/ 4 levels "Alone","Friends",..: 1 2 2 2 2 2 2 3 3 3 ...
## $ weather : Factor w/ 3 levels "Rainy","Snowy",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ temperature : Factor w/ 3 levels "30","55","80": 2 3 3 3 3 3 2 3 3 3 ...
## $ coupon : Factor w/ 5 levels "Bar","Carry out & Take away",..: 4 3 2 3 3 4 2 4 2 1 ...
## $ gender : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
## $ age : Factor w/ 8 levels "21","26","31",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ maritalStatus : Factor w/ 5 levels "Divorced","Married partner",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ has_children : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
## $ education : Factor w/ 6 levels "Associates degree",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ income : Factor w/ 9 levels "$100000 or More",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Bar : Factor w/ 5 levels "1~3","4~8","gt8",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ CoffeeHouse : Factor w/ 5 levels "1~3","4~8","gt8",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ CarryAway : Factor w/ 5 levels "1~3","4~8","gt8",..: 1 4 4 4 4 4 1 4 4 4 ...
## $ RestaurantLessThan20: Factor w/ 5 levels "1~3","4~8","gt8",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Restaurant20To50 : Factor w/ 5 levels "1~3","4~8","gt8",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ toCoupon_GEQ15min : Factor w/ 2 levels "0","1": 1 1 2 2 2 2 2 2 2 2 ...
## $ toCoupon_GEQ25min : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ direction_same : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ Y : Factor w/ 2 levels "0","1": 2 1 2 1 1 2 2 2 2 1 ...
## $ occupation_class : Factor w/ 8 levels "Craft and related trades workers",..: 8 8 8 8 8 8 8 8 8 8 ...
## $ expiration_weightage: num [1:12684, 1] 1.334 0.111 0.111 0.111 1.334 ...
## ..- attr(*, "scaled:scale")= num 18Run model on complete dataset
##
## Call:
## randomForest(formula = Y ~ ., data = train)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 4
##
## OOB estimate of error rate: 25.19%
## Confusion matrix:
## 0 1 class.error
## 0 2452 1216 0.3315158
## 1 925 3906 0.1914717
Tune parameters
# Perform training with parameters
rf_classifier = randomForest(Y ~., data = train, ntree=200, mtry=5, importance=TRUE)
rf_classifier##
## Call:
## randomForest(formula = Y ~ ., data = train, ntree = 200, mtry = 5, importance = TRUE)
## Type of random forest: classification
## Number of trees: 200
## No. of variables tried at each split: 5
##
## OOB estimate of error rate: 25.17%
## Confusion matrix:
## 0 1 class.error
## 0 2477 1191 0.3247001
## 1 948 3883 0.1962327 # plot(rf_classifier)
varImpPlot(rf_classifier)
Most important variables: Coupon, coffee house, age, income, occupation, bar, restaurantlessthan20, expiration
MeanDecreaseAccuracy: gives a rough estimate of the loss in prediction performance when that particular variable is omitted from the training set. Caveat: if two variables are somewhat redundant, then omitting one of them may not lead to massive gains in prediction performance, but would make the second variable more important.
MeanDecreaseGini: GINI is a measure of node impurity. Think of it like this, if you use this feature to split the data, how pure will the nodes be? Highest purity means that each node contains only elements of a single class. Assessing the decrease in GINI when that feature is omitted leads to an understanding of how important that feature is to split the data correctly.
Model Evaluation
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1212 410
## 1 594 1969
##
## Accuracy : 0.7601
## 95% CI : (0.7469, 0.773)
## No Information Rate : 0.5685
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.505
##
## Mcnemar's Test P-Value : 7.676e-09
##
## Sensitivity : 0.8277
## Specificity : 0.6711
## Pos Pred Value : 0.7682
## Neg Pred Value : 0.7472
## Prevalence : 0.5685
## Detection Rate : 0.4705
## Detection Prevalence : 0.6124
## Balanced Accuracy : 0.7494
##
## 'Positive' Class : 1
## # classification perfomance
# auc(test$Y,as.numeric(as.character(prediction_for_table)))
test_roc = roc(test$Y , as.numeric(as.character(prediction_for_table)), plot = TRUE, print.auc = TRUE)
BOOSTING
- An ensemble modeling technique which attempts to build a strong classifier from the number of weak classifiers. It is done building a model by using weak models in series. Firstly, a model is built from the training data. Then the second model is built which tries to correct the errors present in the first model. This procedure is continued and models are added until either the complete training data set is predicted correctly or the maximum number of models are added
- Trees are grown sequentially & each tree uses information from previous tree
Train Model
# gaussian - regression, bernoulli - classification
rf_boost <- gbm(Y~.,data=train_boost,distribution="bernoulli",n.trees=5000,shrinkage=0.01,interaction.depth=1, cv.folds=3) #, verbose=F
best.iter = gbm.perf(rf_boost, method="cv") #Check the best iteration number
# plot - Summary of the model results, with the importance plot of predictors
par(mar = c(5, 8, 1, 1))
summary(rf_boost, cBars = 10,las = 2)
## var rel.inf
## coupon coupon 30.1105941
## CoffeeHouse CoffeeHouse 11.3105901
## expiration_weightage expiration_weightage 9.7430308
## income income 6.9154155
## destination destination 5.6196868
## weather weather 5.5998788
## passanger passanger 5.2079665
## occupation_class occupation_class 4.0144568
## age age 3.7807999
## education education 3.1227566
## Restaurant20To50 Restaurant20To50 1.9246713
## toCoupon_GEQ25min toCoupon_GEQ25min 1.7884993
## Bar Bar 1.7506133
## direction_same direction_same 1.7388944
## maritalStatus maritalStatus 1.5690136
## RestaurantLessThan20 RestaurantLessThan20 1.2921685
## gender gender 1.2840792
## CarryAway CarryAway 1.1991101
## temperature temperature 0.9763819
## toCoupon_GEQ15min toCoupon_GEQ15min 0.9205985
## has_children has_children 0.1307938## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 1010 546
## 1 796 1833
##
## Accuracy : 0.6793
## 95% CI : (0.6649, 0.6935)
## No Information Rate : 0.5685
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.3353
##
## Mcnemar's Test P-Value : 1.068e-11
##
## Sensitivity : 0.7705
## Specificity : 0.5592
## Pos Pred Value : 0.6972
## Neg Pred Value : 0.6491
## Prevalence : 0.5685
## Detection Rate : 0.4380
## Detection Prevalence : 0.6282
## Balanced Accuracy : 0.6649
##
## 'Positive' Class : 1
## Tune model
# Using the caret package the get the model preformance in the best iteration
set.seed(123)
fitControl = trainControl(method="cv", number=3, returnResamp = "all")
train_boost$Y <- as.factor(train_boost$Y)
test_boost$Y <- as.factor(test_boost$Y)
model2 = train(Y~., data=train_boost, method="gbm",distribution="bernoulli", trControl=fitControl,
verbose=F, tuneGrid=data.frame(.n.trees=best.iter, .shrinkage=0.01, .interaction.depth=1, .n.minobsinnode=1)) #
# model2
# confusionMatrix(model2)
mPred = predict(model2, test_boost[,-y_col]) #, na.action = na.pass
# postResample(mPred, test_boost$Y)
# confusionMatrix(mPred, test_boost$Y)
confusionMatrix(
as.factor(mPred),
as.factor(test_boost$Y),
positive = "1"
)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 996 549
## 1 810 1830
##
## Accuracy : 0.6753
## 95% CI : (0.6608, 0.6894)
## No Information Rate : 0.5685
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.3264
##
## Mcnemar's Test P-Value : 1.753e-12
##
## Sensitivity : 0.7692
## Specificity : 0.5515
## Pos Pred Value : 0.6932
## Neg Pred Value : 0.6447
## Prevalence : 0.5685
## Detection Rate : 0.4373
## Detection Prevalence : 0.6308
## Balanced Accuracy : 0.6604
##
## 'Positive' Class : 1
## - Performance of boosting can be an effect of noise in the data
- Larger noise leads to lack of bias
- Boosting works on weak classifiers that have high bias and low variation works iteratively on weak learners, and more weightage is given to misclassified learners in next iteration
Unsupervised Learning
CLUSTERING
- Clustering is used to find structure in unlabeled data
- Interested to see if the drivers could be segmented easily based on the data present
Dataset summary
## 'data.frame': 12684 obs. of 23 variables:
## $ destination : Factor w/ 3 levels "Home","No Urgent Place",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ passanger : Factor w/ 4 levels "Alone","Friends",..: 1 2 2 2 2 2 2 3 3 3 ...
## $ weather : Factor w/ 3 levels "Rainy","Snowy",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ temperature : Factor w/ 3 levels "30","55","80": 2 3 3 3 3 3 2 3 3 3 ...
## $ coupon : Factor w/ 5 levels "Bar","Carry out & Take away",..: 4 3 2 3 3 4 2 4 2 1 ...
## $ gender : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
## $ age : Factor w/ 8 levels "21","26","31",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ maritalStatus : Factor w/ 5 levels "Divorced","Married partner",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ has_children : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
## $ education : Factor w/ 6 levels "Associates degree",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ income : Factor w/ 9 levels "$100000 or More",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Bar : Factor w/ 5 levels "1~3","4~8","gt8",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ CoffeeHouse : Factor w/ 5 levels "1~3","4~8","gt8",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ CarryAway : Factor w/ 5 levels "1~3","4~8","gt8",..: 1 4 4 4 2 4 1 4 4 4 ...
## $ RestaurantLessThan20: Factor w/ 5 levels "1~3","4~8","gt8",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Restaurant20To50 : Factor w/ 5 levels "1~3","4~8","gt8",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ toCoupon_GEQ15min : Factor w/ 2 levels "0","1": 1 1 2 2 2 2 2 2 2 2 ...
## $ toCoupon_GEQ25min : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ direction_same : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ Y : Factor w/ 2 levels "0","1": 2 1 2 1 1 2 2 2 2 1 ...
## $ occupation_class : Factor w/ 8 levels "Craft and related trades workers",..: 8 8 8 8 8 8 8 8 8 8 ...
## $ expiration_weightage: num 1.334 0.111 0.111 0.111 1.334 ...
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...Gower distance
(for mixed type data - categorical and numerical)
Cons: Requires an NxN distance matrix to be calculated, which is computationally intensive. One possible solution for this is to sample data, cluster the smaller sample, then treat the clustered sample as training data for k Nearest Neighbors and “classify” the rest of the data.
We create a sample of 1000 observations.
library(cluster)
library(Rtsne)
gower_distance <- function(datafr) {
gower_dist <- daisy(datafr,
metric = "gower",
type = list(logratio = 3))
return(gower_dist)
}Clustering for entire driving scenario (with all the variables) Gower distance
set.seed(123)
#removing Y and ID
gower_dist <- gower_distance(sample_clust[,-c(20,23)])## Warning in daisy(datafr, metric = "gower", type = list(logratio = 3)): binary
## variable(s) 21 treated as interval scaledsummary(gower_dist)## 499500 dissimilarities, summarized :
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.5238 0.6190 0.6140 0.7143 0.9524
## Metric : mixed ; Types = N, N, N, N, N, N, N, N, N, N, N, N, N, N, N, N, N, N, N, N, I
## Number of objects : 1000Gower distance = 0 minimum. 1 - maximum distance.
We have more dissimilar observaions in the dataset.
# Output most similar pair
gower_mat <- as.matrix(gower_dist)
sample_clust[
which(gower_mat == min(gower_mat[gower_mat != min(gower_mat)]),
arr.ind = TRUE)[1, ], ]## destination passanger weather temperature coupon gender
## 2777 Work Alone Sunny 55 Restaurant(<20) Female
## 2757 No Urgent Place Alone Sunny 55 Restaurant(<20) Female
## age maritalStatus has_children education income Bar
## 2777 50plus Married partner 1 Bachelors degree $50000 - $62499 never
## 2757 50plus Married partner 1 Bachelors degree $50000 - $62499 never
## CoffeeHouse CarryAway RestaurantLessThan20 Restaurant20To50
## 2777 never less1 less1 less1
## 2757 never less1 less1 less1
## toCoupon_GEQ15min toCoupon_GEQ25min direction_same Y occupation_class
## 2777 0 0 0 1 Professionals
## 2757 0 0 0 1 Professionals
## expiration_weightage ID
## 2777 1.33374 2777
## 2757 1.33374 2757# Output most dissimilar pair
sample_clust[
which(gower_mat == max(gower_mat[gower_mat != max(gower_mat)]),
arr.ind = TRUE)[1, ], ]## destination passanger weather temperature coupon gender age
## 466 Home Alone Sunny 55 Bar Male 31
## 2463 No Urgent Place Partner Sunny 80 Coffee House Female 21
## maritalStatus has_children education income
## 466 Single 1 High School Graduate $50000 - $62499
## 2463 Unmarried partner 0 Some college - no degree $12500 - $24999
## Bar CoffeeHouse CarryAway RestaurantLessThan20 Restaurant20To50
## 466 gt8 gt8 gt8 gt8 gt8
## 2463 1~3 less1 1~3 1~3 less1
## toCoupon_GEQ15min toCoupon_GEQ25min direction_same Y occupation_class
## 466 0 0 1 1 Professionals
## 2463 1 0 0 1 Unemployed
## expiration_weightage ID
## 466 1.333740 466
## 2463 0.111145 2463
K-medoids Clustering
Partitioning around medoids with custom distance matrix
The k-medoids problem is a clustering problem similar to k-means.In contrast to the k-means algorithm, k-medoids chooses actual data points as centers (medoids or exemplars), and thereby allows for greater interpretability of the cluster centers than in k-means.
k-medoids can be used with arbitrary dissimilarity measures, whereas k-means generally requires Euclidean distance for efficient solutions. Because k-medoids minimizes a sum of pairwise dissimilarities instead of a sum of squared Euclidean distances, it is more robust to noise and outliers than k-means.
Selecting the number of clusters k using silhouette width, an internal validation metric which is an aggregated measure of how similar an observation is to its own cluster compared its closest neighboring cluster. The metric can range from -1 to 1, where higher values are better.
library(Rtsne)
# Calculate silhouette width for many k using PAM
sil_width <- c(NA)
for(i in 2:10){
pam_fit <- pam(gower_dist,
diss = TRUE,
k = i)
sil_width[i] <- pam_fit$silinfo$avg.width
}
# Plot sihouette width (higher is better)
plot(1:10, sil_width,
xlab = "Number of clusters",
ylab = "Silhouette Width")
lines(1:10, sil_width)
Interpret the clusters by running summary on each cluster:
In cluster 1 (430): higher frequency for Alone, Male, 21, Single, no_children, Some college - no degree/Bachelors degree, lower income, Student
In cluster 2 (570): higher frequency for No Urgent Place, Friends/Alone, Female, 50plus, Married_Partner, has_children, Bachelors degree/Graduate, higher income,Restaurant20To50, toCoupon_GEQ15min, Professionals
Other characteristics are overlapping for both the clusters, so we can explore if only the personal characteriscts and not the driving scenario lead to the cluster separation.
Distribution of Y is the same in the clusters (58% accepted in the first group and 55% in the second).
set.seed(123)
pam_fit <- pam(gower_dist, diss = TRUE, k = 2)
pam_results <- sample_clust %>%
dplyr::select(-ID) %>%
mutate(cluster = pam_fit$clustering) %>%
group_by(cluster) %>%
do(the_summary = summary(.))
#pam_results$the_summary sample_clust[pam_fit$medoids, ]## destination passanger weather temperature coupon gender
## 12584 Home Alone Sunny 80 Restaurant(20-50) Male
## 9753 No Urgent Place Friends Sunny 80 Coffee House Female
## age maritalStatus has_children education
## 12584 21 Single 0 Some college - no degree
## 9753 50plus Married partner 1 Bachelors degree
## income Bar CoffeeHouse CarryAway RestaurantLessThan20
## 12584 $12500 - $24999 never never 1~3 1~3
## 9753 $100000 or More never less1 4~8 1~3
## Restaurant20To50 toCoupon_GEQ15min toCoupon_GEQ25min direction_same Y
## 12584 less1 0 0 0 0
## 9753 1~3 1 0 0 1
## occupation_class expiration_weightage ID
## 12584 Student 1.33374 12584
## 9753 Professionals 1.33374 9753One way to visualize many variables in a lower dimensional space is with t-distributed stochastic neighborhood embedding, or t-SNE. A dimensionality reduction technique that tries to preserve local structure so as to make clusters visible in a 2D or 3D visualization.
tsne_obj <- Rtsne(gower_dist, is_distance = TRUE)
tsne_data <- tsne_obj$Y %>%
data.frame() %>%
setNames(c("X", "Y")) %>%
mutate(cluster = factor(pam_fit$clustering),
name = sample_clust$ID)
ggplot(aes(x = X, y = Y), data = tsne_data) +
geom_point(aes(color = cluster))
Clusters are overlapping. Blue points are closer to blue cluster medoid for the metrics, however, it’s impossible to show it in a 2D space due to the fact that only 18% of variance of the data explained after dim reduction.
Agglomerative clustering
aggl.clust.c <- hclust(gower_dist, method = "complete")
plot(aggl.clust.c,
main = "Agglomerative, complete linkages")
#
# aggl.clust.c <- hclust(gower_dist, method = "average")
# plot(aggl.clust.c,
# main = "Agglomerative, average linkages")
#
# aggl.clust.c <- hclust(gower_dist, method = "single")
# plot(aggl.clust.c,
# main = "Agglomerative, single linkages")
aggl.clust.w <- hclust(gower_dist, method = "ward.D2")
plot(aggl.clust.w,
main = "Agglomerative, ward linkages")
# hierarchical clustering using Ward linkage
dendro <- as.dendrogram(aggl.clust.w)
dendro.col <- dendro %>%
set("branches_k_color", k = 2, value = c( "gold3", "darkcyan")) %>%
set("branches_lwd", 0.6) %>%
set("labels_colors",
value = c("darkslategray")) %>%
set("labels_cex", 0.5)
ggd1 <- as.ggdend(dendro.col)
ggplot(ggd1, theme = theme_minimal()) +
labs(x = "Num. observations", y = "Height", title = "Dendrogram, k = 2")
groups <- cutree(aggl.clust.w, k=2) # cut tree into 5 clusters
clusplot(sample_clust, groups, color=TRUE, shade=TRUE,
labels=2, lines=0, main= 'Driving scenario clusters')
Clustering evaluation - Rand index
Rand Index looks at similarity between any two clustering methods. The Rand Index gives a value between 0 and 1,where 1 means the two clustering outcomes match identically
Another strategy - groups based on personal characteristics
Dataset Information
## 'data.frame': 1000 obs. of 9 variables:
## $ gender : Factor w/ 2 levels "Female","Male": 1 2 1 1 1 2 2 1 2 2 ...
## $ age : Factor w/ 8 levels "21","26","31",..: 1 1 2 3 2 1 5 6 2 1 ...
## $ maritalStatus : Factor w/ 5 levels "Divorced","Married partner",..: 4 3 4 2 4 4 3 1 3 3 ...
## $ has_children : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 1 ...
## $ education : Factor w/ 6 levels "Associates degree",..: 5 5 2 5 5 5 4 2 4 5 ...
## $ income : Factor w/ 9 levels "$100000 or More",..: 2 1 4 1 5 3 2 3 2 9 ...
## $ Y : Factor w/ 2 levels "0","1": 2 2 2 1 2 1 2 2 2 2 ...
## $ occupation_class: Factor w/ 8 levels "Craft and related trades workers",..: 8 6 7 3 2 5 1 2 5 6 ...
## $ ID : int 2463 2511 10419 8718 12483 2986 1842 9334 3371 11638 ...
Group by characteristics of coupon only
## 'data.frame': 1000 obs. of 6 variables:
## $ coupon : Factor w/ 5 levels "Bar","Carry out & Take away",..: 3 3 3 5 1 2 3 4 3 1 ...
## $ toCoupon_GEQ15min : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 1 1 1 2 ...
## $ toCoupon_GEQ25min : Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 1 1 1 ...
## $ Y : Factor w/ 2 levels "0","1": 2 2 2 1 2 1 2 2 2 2 ...
## $ expiration_weightage: num 0.111 0.111 1.334 0.111 0.111 ...
## $ ID : int 2463 2511 10419 8718 12483 2986 1842 9334 3371 11638 ...
pam_results$the_summary ## [[1]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0: 0 0:59 0: 0
## Carry out & Take away: 0 1:66 1: 7 1:66
## Coffee House :66
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :1
## 1st Qu.:0.1111 1st Qu.:1
## Median :1.3337 Median :1
## Mean :0.7595 Mean :1
## 3rd Qu.:1.3337 3rd Qu.:1
## Max. :1.3337 Max. :1
##
## [[2]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0:85 0:85 0: 0
## Carry out & Take away: 0 1: 0 1: 0 1:85
## Coffee House :85
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :2
## 1st Qu.:0.1111 1st Qu.:2
## Median :0.1111 Median :2
## Mean :0.6865 Mean :2
## 3rd Qu.:1.3337 3rd Qu.:2
## Max. :1.3337 Max. :2
##
## [[3]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0: 0 0:43 0:47
## Carry out & Take away: 0 1:47 1: 4 1: 0
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) :47
##
## expiration_weightage cluster
## Min. :0.1111 Min. :3
## 1st Qu.:0.1111 1st Qu.:3
## Median :1.3337 Median :3
## Mean :0.9956 Mean :3
## 3rd Qu.:1.3337 3rd Qu.:3
## Max. :1.3337 Max. :3
##
## [[4]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar :30 0: 0 0:20 0: 0
## Carry out & Take away: 0 1:30 1:10 1:30
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :4
## 1st Qu.:0.4168 1st Qu.:4
## Median :1.3337 Median :4
## Mean :1.0077 Mean :4
## 3rd Qu.:1.3337 3rd Qu.:4
## Max. :1.3337 Max. :4
##
## [[5]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0: 0 0:25 0:28
## Carry out & Take away:28 1:28 1: 3 1: 0
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :5
## 1st Qu.:0.1111 1st Qu.:5
## Median :0.1111 Median :5
## Mean :0.5041 Mean :5
## 3rd Qu.:1.3337 3rd Qu.:5
## Max. :1.3337 Max. :5
##
## [[6]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0:95 0:95 0: 0
## Carry out & Take away: 0 1: 0 1: 0 1:95
## Coffee House : 0
## Restaurant(<20) :94
## Restaurant(20-50) : 1
##
## expiration_weightage cluster
## Min. :0.1111 Min. :6
## 1st Qu.:0.1111 1st Qu.:6
## Median :1.3337 Median :6
## Mean :0.9090 Mean :6
## 3rd Qu.:1.3337 3rd Qu.:6
## Max. :1.3337 Max. :6
##
## [[7]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0: 0 0:76 0: 0
## Carry out & Take away:101 1:101 1:25 1:101
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :7
## 1st Qu.:0.1111 1st Qu.:7
## Median :1.3337 Median :7
## Mean :0.9464 Mean :7
## 3rd Qu.:1.3337 3rd Qu.:7
## Max. :1.3337 Max. :7
##
## [[8]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0: 0 0:85 0:97
## Carry out & Take away: 0 1:97 1:12 1: 0
## Coffee House :97
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :8
## 1st Qu.:0.1111 1st Qu.:8
## Median :0.1111 Median :8
## Mean :0.5397 Mean :8
## 3rd Qu.:1.3337 3rd Qu.:8
## Max. :1.3337 Max. :8
##
## [[9]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar :27 0:28 0:28 0: 0
## Carry out & Take away: 0 1: 0 1: 0 1:28
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) : 1
##
## expiration_weightage cluster
## Min. :0.1111 Min. :9
## 1st Qu.:1.3337 1st Qu.:9
## Median :1.3337 Median :9
## Mean :1.2027 Mean :9
## 3rd Qu.:1.3337 3rd Qu.:9
## Max. :1.3337 Max. :9
##
## [[10]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar :23 0:30 0:30 0:30
## Carry out & Take away: 7 1: 0 1: 0 1: 0
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :10
## 1st Qu.:0.1111 1st Qu.:10
## Median :1.3337 Median :10
## Mean :0.8855 Mean :10
## 3rd Qu.:1.3337 3rd Qu.:10
## Max. :1.3337 Max. :10
##
## [[11]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0:61 0:61 0: 0
## Carry out & Take away:60 1: 0 1: 0 1:61
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) : 1
##
## expiration_weightage cluster
## Min. :0.1111 Min. :11
## 1st Qu.:0.1111 1st Qu.:11
## Median :1.3337 Median :11
## Mean :0.8527 Mean :11
## 3rd Qu.:1.3337 3rd Qu.:11
## Max. :1.3337 Max. :11
##
## [[12]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0: 0 0:52 0: 0
## Carry out & Take away: 0 1:58 1: 6 1:58
## Coffee House : 0
## Restaurant(<20) :58
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :12
## 1st Qu.:0.1111 1st Qu.:12
## Median :0.1111 Median :12
## Mean :0.6170 Mean :12
## 3rd Qu.:1.3337 3rd Qu.:12
## Max. :1.3337 Max. :12
##
## [[13]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0:30 0:30 0:30
## Carry out & Take away: 1 1: 0 1: 0 1: 0
## Coffee House : 0
## Restaurant(<20) :29
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :13
## 1st Qu.:0.1111 1st Qu.:13
## Median :0.1111 Median :13
## Mean :0.5187 Mean :13
## 3rd Qu.:1.3337 3rd Qu.:13
## Max. :1.3337 Max. :13
##
## [[14]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0:42 0:42 0:30
## Carry out & Take away: 8 1: 0 1: 0 1:12
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) :34
##
## expiration_weightage cluster
## Min. :0.1111 Min. :14
## 1st Qu.:0.1111 1st Qu.:14
## Median :1.3337 Median :14
## Mean :0.8389 Mean :14
## 3rd Qu.:1.3337 3rd Qu.:14
## Max. :1.3337 Max. :14
##
## [[15]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar :74 0: 0 0:45 0:74
## Carry out & Take away: 0 1:74 1:29 1: 0
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :15
## 1st Qu.:0.1111 1st Qu.:15
## Median :1.3337 Median :15
## Mean :0.9868 Mean :15
## 3rd Qu.:1.3337 3rd Qu.:15
## Max. :1.3337 Max. :15
##
## [[16]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0: 2 0:24 0: 0
## Carry out & Take away: 0 1:26 1: 4 1:28
## Coffee House : 0
## Restaurant(<20) : 0
## Restaurant(20-50) :28
##
## expiration_weightage cluster
## Min. :0.1111 Min. :16
## 1st Qu.:1.3337 1st Qu.:16
## Median :1.3337 Median :16
## Mean :1.1154 Mean :16
## 3rd Qu.:1.3337 3rd Qu.:16
## Max. :1.3337 Max. :16
##
## [[17]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0: 0 0:17 0:34
## Carry out & Take away: 0 1:34 1:17 1: 0
## Coffee House : 0
## Restaurant(<20) :34
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :17
## 1st Qu.:0.1111 1st Qu.:17
## Median :0.1111 Median :17
## Mean :0.3269 Mean :17
## 3rd Qu.:0.1111 3rd Qu.:17
## Max. :1.3337 Max. :17
##
## [[18]]
## coupon toCoupon_GEQ15min toCoupon_GEQ25min Y
## Bar : 0 0:66 0:66 0:66
## Carry out & Take away: 3 1: 0 1: 0 1: 0
## Coffee House :63
## Restaurant(<20) : 0
## Restaurant(20-50) : 0
##
## expiration_weightage cluster
## Min. :0.1111 Min. :18
## 1st Qu.:0.1111 1st Qu.:18
## Median :1.3337 Median :18
## Mean :0.8336 Mean :18
## 3rd Qu.:1.3337 3rd Qu.:18
## Max. :1.3337 Max. :18