library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.1
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.4.2 v purrr 1.0.1
## v tibble 3.2.1 v dplyr 1.1.2
## v tidyr 1.3.0 v stringr 1.5.0
## v readr 2.1.2 v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.3
## Warning: package 'tibble' was built under R version 4.1.3
## Warning: package 'tidyr' was built under R version 4.1.3
## Warning: package 'readr' was built under R version 4.1.3
## Warning: package 'purrr' was built under R version 4.1.3
## Warning: package 'dplyr' was built under R version 4.1.3
## Warning: package 'stringr' was built under R version 4.1.3
## Warning: package 'forcats' was built under R version 4.1.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(caret)
## Warning: package 'caret' was built under R version 4.1.2
## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.1.1
library(psych)
## Warning: package 'psych' was built under R version 4.1.2
##
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(pROC)
## Warning: package 'pROC' was built under R version 4.1.2
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
Gender <- sample(c("M","F"), size = 600, replace = T, prob = c(0.45, 0.55))
City <- sample(x = c("Lisbon", "NY", "London", "LA", "Paris"),
size = 600, replace = T, prob = c(0.24, 0.34, 0.15, 0.22, 0.05))
Marital <- sample(x = c("Single", "Married", "Divorced", "Widowed"),
replace = T, prob = c(0.3, 0.50, 0.15, 0.05), size = 600)
Happiness <- rnorm(n = 600, mean = 70, sd = 15) %>% floor()
Working <- sample(x = c(0,1), replace = T, prob = c(0.25, 0.75), size = 600)
data <- data.frame(Working, Gender, Marital, Happiness, City)
data %>% head(50) %>% knitr::kable(digits = 3)
| 1 |
M |
Single |
82 |
Lisbon |
| 0 |
M |
Married |
69 |
LA |
| 1 |
M |
Single |
108 |
NY |
| 1 |
M |
Married |
87 |
Lisbon |
| 1 |
F |
Married |
70 |
London |
| 1 |
F |
Married |
42 |
London |
| 0 |
F |
Married |
87 |
NY |
| 1 |
F |
Divorced |
79 |
LA |
| 0 |
F |
Divorced |
72 |
LA |
| 1 |
M |
Divorced |
86 |
NY |
| 1 |
F |
Married |
70 |
NY |
| 1 |
F |
Single |
64 |
LA |
| 1 |
F |
Married |
81 |
London |
| 1 |
M |
Married |
77 |
LA |
| 1 |
M |
Divorced |
70 |
Lisbon |
| 1 |
M |
Married |
67 |
NY |
| 0 |
M |
Married |
61 |
NY |
| 1 |
F |
Married |
77 |
LA |
| 1 |
F |
Married |
59 |
Lisbon |
| 1 |
F |
Married |
68 |
NY |
| 1 |
F |
Single |
104 |
Lisbon |
| 0 |
F |
Married |
68 |
LA |
| 1 |
M |
Widowed |
68 |
Lisbon |
| 1 |
F |
Single |
59 |
Lisbon |
| 1 |
F |
Widowed |
63 |
Lisbon |
| 1 |
M |
Married |
55 |
LA |
| 1 |
M |
Divorced |
57 |
NY |
| 0 |
F |
Married |
75 |
NY |
| 1 |
F |
Married |
62 |
Lisbon |
| 1 |
M |
Widowed |
101 |
Lisbon |
| 1 |
M |
Divorced |
66 |
NY |
| 1 |
M |
Married |
54 |
London |
| 0 |
M |
Married |
55 |
Lisbon |
| 0 |
F |
Married |
36 |
Lisbon |
| 1 |
M |
Single |
40 |
NY |
| 1 |
F |
Married |
68 |
NY |
| 1 |
M |
Single |
55 |
NY |
| 1 |
F |
Divorced |
61 |
LA |
| 1 |
M |
Single |
42 |
LA |
| 1 |
F |
Married |
65 |
London |
| 0 |
F |
Single |
70 |
Lisbon |
| 1 |
F |
Single |
51 |
NY |
| 1 |
F |
Single |
52 |
Lisbon |
| 0 |
M |
Married |
73 |
NY |
| 1 |
M |
Single |
64 |
Lisbon |
| 0 |
M |
Single |
51 |
London |
| 0 |
M |
Divorced |
79 |
LA |
| 1 |
F |
Single |
77 |
LA |
| 1 |
F |
Single |
79 |
London |
| 1 |
F |
Married |
66 |
Lisbon |
### Descriptive Statistics ###
data$working_status <-
ifelse(data$Working == 0, "No", "Yes")
ggplot(data, aes(x = working_status)) +
geom_bar(fill = c("red", "green"), col = "black", width = 0.4) +
stat_count(geom = "text", aes(label = ..count..), vjust = -0.5, hjust = .5) +
ggtitle("Distribution of Working Status") +
scale_y_continuous(limits = c(0, max(table(data$Working)) + 50))
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## i Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
data %>% group_by(Working) %>% summarise("Count"=n()) %>% knitr::kable()
ggplot(data, aes(x = working_status, fill = Gender)) +
geom_bar(position = "dodge") +
ggtitle("Distribution of Gender per working Status")
ggplot(data, aes(x = working_status, fill = Marital)) +
geom_bar(position = "dodge") +
ggtitle("Distribution of Marital status per working status")
ggplot(data, aes(x = working_status, fill = City)) +
geom_bar(position = "dodge") +
ggtitle("Distribution of City per working status")
col1 <- colorRampPalette(colors = c("4040555", "7788789", "99669695"))
ggplot(data, aes(x = Happiness)) +
geom_histogram(fill = col1(60), col = "black", size = 0.4) +
facet_wrap(~working_status) + ggtitle("Distribution of Happiness by Working Status")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
data %>% group_by(working_status) %>%
summarise("Count"=n(), "Mean of Happiness"=mean(Happiness)) %>%
knitr::kable(digits = 3)
| No |
146 |
68.178 |
| Yes |
454 |
68.978 |
data <- select(data, -working_status)
### End of Descriptive Statistics Part ###
### Logistic Regression Part ###
model <- glm(Working ~ . , data, family = "binomial")
summary(model)
##
## Call:
## glm(formula = Working ~ ., family = "binomial", data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0346 0.5309 0.7421 0.7995 0.8658
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.625357 0.547236 1.143 0.2531
## GenderM -0.009749 0.194445 -0.050 0.9600
## MaritalMarried 0.067327 0.288639 0.233 0.8156
## MaritalSingle 0.087818 0.309341 0.284 0.7765
## MaritalWidowed 0.305028 0.526415 0.579 0.5623
## Happiness 0.003319 0.006418 0.517 0.6050
## CityLisbon 0.705767 0.304105 2.321 0.0203 *
## CityLondon 0.084706 0.312583 0.271 0.7864
## CityNY 0.093238 0.252353 0.369 0.7118
## CityParis 0.371800 0.503504 0.738 0.4603
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 665.87 on 599 degrees of freedom
## Residual deviance: 657.75 on 590 degrees of freedom
## AIC: 677.75
##
## Number of Fisher Scoring iterations: 4
data$prob <- predict(model, type = "response")
describe(data$prob)
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 600 0.76 0.05 0.74 0.75 0.03 0.69 0.88 0.19 0.86 -0.74 0
## Several options for prediction:
### Prediction based upon 50%
data$result <- ifelse(data$prob <= 0.50, 0, 1)
table(data$Working, data$result)
##
## 1
## 0 146
## 1 454
confusionMatrix(reference = as.factor(data$Working), as.factor(data$result))
## Warning in confusionMatrix.default(reference = as.factor(data$Working), : Levels
## are not in the same order for reference and data. Refactoring data to match.
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 0 0
## 1 146 454
##
## Accuracy : 0.7567
## 95% CI : (0.7203, 0.7905)
## No Information Rate : 0.7567
## P-Value [Acc > NIR] : 0.5222
##
## Kappa : 0
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.0000
## Specificity : 1.0000
## Pos Pred Value : NaN
## Neg Pred Value : 0.7567
## Prevalence : 0.2433
## Detection Rate : 0.0000
## Detection Prevalence : 0.0000
## Balanced Accuracy : 0.5000
##
## 'Positive' Class : 0
##
### Prediction based upon mean probability of prediction value
mean_prob <- data$prob %>% mean()
data$result <- ifelse(data$prob <= mean_prob, 0, 1)
table(data$Working, data$result)
##
## 0 1
## 0 114 32
## 1 298 156
confusionMatrix(reference = as.factor(data$Working), as.factor(data$result))
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 114 298
## 1 32 156
##
## Accuracy : 0.45
## 95% CI : (0.4097, 0.4908)
## No Information Rate : 0.7567
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.0769
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.7808
## Specificity : 0.3436
## Pos Pred Value : 0.2767
## Neg Pred Value : 0.8298
## Prevalence : 0.2433
## Detection Rate : 0.1900
## Detection Prevalence : 0.6867
## Balanced Accuracy : 0.5622
##
## 'Positive' Class : 0
##
### Prediction based upon mean probability of prediction value
#### ROC Curve
plot(roc(data$Working, data$prob), col = "red", main = "ROC Curve")
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
#### AUC value
auc_val <- roc(data$Working, data$prob) %>% auc()
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
auc_val
## Area under the curve: 0.5598
data$result <- ifelse(data$prob <= auc_val, 0, 1)
table(data$Working, data$result)
##
## 1
## 0 146
## 1 454
confusionMatrix(reference = as.factor(data$Working), as.factor(data$result))
## Warning in confusionMatrix.default(reference = as.factor(data$Working), : Levels
## are not in the same order for reference and data. Refactoring data to match.
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 0 0
## 1 146 454
##
## Accuracy : 0.7567
## 95% CI : (0.7203, 0.7905)
## No Information Rate : 0.7567
## P-Value [Acc > NIR] : 0.5222
##
## Kappa : 0
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.0000
## Specificity : 1.0000
## Pos Pred Value : NaN
## Neg Pred Value : 0.7567
## Prevalence : 0.2433
## Detection Rate : 0.0000
## Detection Prevalence : 0.0000
## Balanced Accuracy : 0.5000
##
## 'Positive' Class : 0
##