Repository navigation
Expand file tree
/
Copy pathHR Predictive Analytics Tutorial 1-4.r
More file actions
473 lines (325 loc) · 14.4 KB
/
Copy pathHR Predictive Analytics Tutorial 1-4.r
File metadata and controls
473 lines (325 loc) · 14.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
#Full description of code and analysis can be found at: https://www.littlemissdata.com/blog/predictive-analytics-tutorial-part-1
# If you are using local R install:
#You can download the file from https://raw.githubusercontent.com/lgellis/Predictive_Analytics_Tutorial
#and use the following line of code to import:
#hr = read.table("/filepath/HR_comma_sep.csv", head=TRUE, sep=",")
#OR
#You can import the file directly from the github repo with these line of code:
#install.packages("data.table")
#library(data.table)
#hr= fread('https://raw.githubusercontent.com/lgellis/Predictive_Analytics_Tutorial/master/HR_comma_sep.csv')
#If you are using DSX, import the file place your cursor in cell and "insert to code" as per the instructions here:
#https://www.littlemissdata.com/blog/predictive-analytics-tutorial-part-1
head(df.data.2)
#rename data frame
hr <- df.data.2
#view columns
names(hr)
#quick view
head(hr)
#dim shows number of rows and columns
dim(hr)
#str shows the structure of each column including the type of data and sample values
str(hr)
#show summary statistics for each column. For numeric: min, max, median, mean, 2st and 3rd quartile. For factors show the counts of each.
summary(hr)
#install packages - do this one time
install.packages("data.table")
install.packages("corrplot")
install.packages("ggplot2")
install.packages("gcookbook")
install.packages("caret")
install.packages("hexbin")
install.packages("leaps")
install.packages("plyr")
install.packages("plotly")
install.packages("waffle")
install.packages("dummies")
install.packages("caTools")
install.packages("wesanderson")
install.packages("visreg")
install.packages("car")
install.packages("leaps")
install.packages("MASS")
# Load the relevant libraries - do this every time
library(data.table)
library(corrplot)
library(ggplot2)
library (gcookbook)
library(caret)
library(hexbin)
library(leaps)
library(plyr)
library(plotly)
library(waffle)
library(dummies)
library(caTools)
library(wesanderson)
library(visreg)
library(car)
library(rpart)
library(leaps)
library(MASS)
#attach allows us to reference columns by their name
attach(hr)
#Check Correlations of numeric columns
corMatrix <-cor(hr[1:8], use="complete.obs", method="pearson")
#round to two decimals
round(corMatrix, 2)
corrplot(corMatrix, method="circle")
corrplot(corMatrix, method="square")
corrplot(corMatrix, method="number")
corrplot(corMatrix, method="shade")
corrplot(corMatrix, type = "upper")
corrplot.mixed(corMatrix)
#rename
colNames <- c("satLevel", "lastEval", "numProj", "avgHrs", "timeCpny", "wrkAcdnt", "left", "fiveYrPrmo", "job", "salary")
setnames(hr, colNames)
names(hr)
attach(hr)
corrplot(corMatrix, order = "hclust", addrect = 2, col = heat.colors(100), tl.col = "black")
# Change background color to lightblue
corrplot(corMatrix, type = "upper", order = "hclust", col = c("black", "white"), bg = "lightblue", tl.col = "black")
# Run a histogram for all numeric variables to understand distribution
hist(avgHrs/4, main="Distribution of Average Hours per Week", xlab="Avg Hours", breaks=7, col="lightblue")
hist(satLevel, main="Distribution of Satisfaction Level", xlab="Satisfaction Level", breaks=7, col="lightblue")
hist(lastEval, main="Distribution of Last Evaluations", xlab="Last Eval", breaks=7, col="lightblue")
hist(numProj, main="Distribution of Number of Projects", xlab="Number of Projects", breaks=7, col="lightblue")
# side by side pretty bar
ggplot(data=hr, aes(x=numProj, y=avgHrs, fill=salary)) +
geom_bar(position="dodge", color="black", stat="identity") +
scale_fill_brewer(palette = "Pastel1")
#density plot
qplot(avgHrs/4, data=hr, geom="density", fill=salary, alpha=I(.5),
main="Avg Weekly Hours by Salary Category", xlab="Average Weekly Hours",
ylab="Density")
attach(hr)
#create factor variables with a better format
hr$leftFactor <- factor(left,levels=c(0,1),
labels=c("Did Not Leave Company","Left Company"))
hr$promoFactor <- factor(fiveYrPrmo,levels=c(0,1),
labels=c("Did Not Get Promoted","Did Get Promoted"))
hr$wrkAcdntFactor <- factor(wrkAcdnt,levels=c(0,1),
labels=c("No Accident","Accident"))
attach(hr)
#density plot
qplot(avgHrs/4, data=hr, geom="density", fill=leftFactor, alpha=I(.5),
main="Avg Weekly Hours by Retention", xlab="Average Weekly Hours",
ylab="Density")
#boxplot
boxplot(avgHrs~job,data=hr, main="HR Data",
xlab="Job Title", ylab="Avg Hours", col="lightblue")
#violin plot
hrBox <-ggplot(hr, aes(y=avgHrs, x=job))
hrBox + geom_violin(trim=FALSE, fill="lightblue")
#many dimension charts
qplot(avgHrs/4, timeCpny, data=hr, shape=leftFactor, color=salary, facets=numProj~promoFactor, size=I(3),
xlab="avg hrs", ylab="time at company")
hrScat <-ggplot(hr, aes(x=avgHrs/4, y=satLevel))
hrScat + geom_point()
#make the points more transparent so that it's less intense
hrScat + geom_point(alpha=.01)
hrScat + stat_bin2d
hrScat + stat_binhex()
hrScat + stat_binhex() + scale_fill_gradient(low="lightblue", high="red")
LEvSL <-ggplot(hr, aes(x=satLevel, y=lastEval))
LEvSL + geom_point()
LEvSL + stat_binhex() + scale_fill_gradient(low="lightblue", high="red")
# Example of a Bagplot
#install.packages("aplpack")
library(aplpack)
bagplot(avgHrs/4, satLevel, xlab="Avg Weekly Hrs", ylab="SatLevel",
main="Bagplot Avg Hours by Sat Level")
attach(hr)
typeof(hr)
hr <- as.data.frame(hr)
typeof(hr)
#create dummy variables for job and salary All values need to be numeric for inclusion in the model
hr2 <- cbind (hr, dummy(job), dummy(salary))
#Ensure dummy variables are created
names(hr2)
head(hr2)
#Log Transforms for: satLevel, lastEval
hr2$satLevelLog <- log(satLevel)
hr2$lastEvalLog <- log(lastEval)
#SQRT Transforms for: satLevel, lastEval
hr2$satLevelSqrt <- sqrt(satLevel)
hr2$lastEvalSqrt <- sqrt(lastEval)
#Scale Transforms for: satLevel, lastEval
hr2$satLevelScale <- scale(satLevel)
hr2$lastEvalScale <- scale(lastEval)
#Peek at new columns
head(hr2)
#Ensure all columns present
names(hr2)
#Create the binary variable for highly rated but low satisfaction users
hr2$greatEvalLowSat <- ifelse(lastEval>0.8 & satLevel <0.2, 1, 0)
#Visualize effects that it has on average hours
attach(hr2)
x <- ggplot(hr2, aes(factor(greatEvalLowSat), avgHrs))
x <- x + geom_boxplot (aes(fill=factor(greatEvalLowSat)), outlier.color="black", outlier.size=1)
x <- x + coord_flip()
x + scale_fill_manual(values=wes_palette(n=2, name="GrandBudapest2"))
#Visualize distribution by job
ggplot(hr, aes(avgHrs, job, z= greatEvalLowSat)) + geom_tile(aes(fill = greatEvalLowSat)) + theme_bw() + scale_fill_gradient(low="#98afc7", high="#646d7e")
#Get distribution table
summary <- count(avgHrs, c('greatEvalLowSat'))
summary
#Starting code credit to: https://stackoverflow.com/questions/17200114/how-to-split-data-into-training-testing-sets-using-sample-function
set.seed(101)
sample = sample.split(hr2, SplitRatio = .75)
train = subset(hr2, sample == TRUE)
test = subset(hr2, sample == FALSE)
#See if the test/train sets have beens split appropriately
dim(train)
dim(test)
#Ensure data characteristics are roughly the same
summary(train)
summary(test)
#1) Throw everything in!
model.lr1 <- lm(avgHrs ~ ., train)
summary(model.lr1)
#adj rsquared - higher the better
summary(model.lr1)
#Adj R-squared - higher is better
#AIC, BIC - lower the better
aic <- AIC(model.lr1)
bic <-BIC(model.lr1)
pred.model.lr1 <- predict(model.lr1, newdata = test) # validation predictions
# mean Prediction Error and Std Error - lower is better
meanPred <- mean((test$avgHrs - pred.model.lr1)^2) # mean prediction error
stdError <- sd((test$avgHrs - pred.model.lr1)^2)/length(test$avgHrs) # std error
str(summary(model.lr1))
### Create a data frame to store the model metrics
# Note stringsAsFactors = FALSE
modelMetrics <- data.frame( "Model" = character(), "adjRsq" = integer(), "AIC"= integer(), "BIC"= integer(), "Mean Prediction Error"= integer(), "Standard Error"= integer(), stringsAsFactors=FALSE)
# Append a row
modelMetrics[nrow(modelMetrics) + 1, ] <- c( "model.lr1", 0.279, aic, bic, meanPred, stdError)
modelMetrics
#2) Backward selection
step <- stepAIC(model.lr1, direction="backward")
step$anova # display results
model.bkwd.lr1 <- lm(avgHrs ~ satLevel + lastEval + numProj + timeCpny + left + satLevelLog +
lastEvalLog + satLevelSqrt + lastEvalSqrt, train)
summary(model.bkwd.lr1)
#Adj R-squared - higher is better
#AIC, BIC - lower the better
aic <- AIC(model.bkwd.lr1)
bic <-BIC(model.bkwd.lr1)
pred.model.bkwd.lr1 <- predict(model.bkwd.lr1, newdata = test) # validation predictions
# mean Prediction Error and Std Error - lower is better
meanPred <- mean((test$avgHrs - pred.model.bkwd.lr1 )^2) # mean prediction error
stdError <- sd((test$avgHrs - pred.model.bkwd.lr1 )^2)/length(test$avgHrs) # std error
# Append a row to our modelMetrics
modelMetrics[nrow(modelMetrics) + 1, ] <- c( "model.bkwd.lr1", 0.2795, aic, bic, meanPred, stdError)
modelMetrics
#3) Create a Model with the highest Correlated Values
#Pull a correlation plot on all of the variables to see what are the most correlated to avgHrs
#remove the non-numeric
corData <- train[ -c(9:13) ]
head(corData)
corMatrix <-cor(corData, use="complete.obs", method="pearson")
#round to two decimals
corrplot(corMatrix, type = "upper", order = "hclust", col = c("black", "white"), bg = "lightblue", tl.col = "black")
#some of the highest correlated are: numProj, lastEvalLog, timeCmpy, left, greatEvalLowSat
# Create a Model with the highest Correlated Values
model.lr2 <- lm(avgHrs ~ numProj + lastEvalLog + timeCpny + left + greatEvalLowSat, train)
summary(model.lr2)
#Adj R-squared - higher is better
#AIC, BIC - lower the better
aic <- AIC(model.lr2)
bic <-BIC(model.lr2)
pred.model.lr2 <- predict(model.lr2, newdata = test) # validation predictions
# mean Prediction Error and Std Error - lower is better
meanPred <- mean((test$avgHrs - pred.model.lr2 )^2) # mean prediction error
stdError <- sd((test$avgHrs - pred.model.lr2 )^2)/length(test$avgHrs) # std error
# Append a row to our modelMetrics
modelMetrics[nrow(modelMetrics) + 1, ] <- c( "model.lr2", 0.2353, aic, bic, meanPred, stdError)
modelMetrics
#4) Best Subsets
model.bestSub=regsubsets(avgHrs ~ ., train, nvmax =25)
summary(model.bestSub)
reg.summary =summary(model.bestSub)
which.min (reg.summary$bic )
which.max (reg.summary$adjr2 )#just for fun
#Plot the variable bic values by number of variables
plot(reg.summary$bic ,xlab=" Number of Variables ",ylab=" BIC",type="l")
points (6, reg.summary$bic [6], col =" red",cex =2, pch =20)
coef(model.bestSub, 6)
bestSubModel = lm(avgHrs ~ satLevel + numProj + timeCpny + NAsupport + NAtechnical + NAhigh, data=train)
bestSubModel
summary(bestSubModel)
#Adj R-squared - higher is better
#AIC, BIC - lower the better
aic <- AIC(bestSubModel)
bic <-BIC(bestSubModel)
pred.bestSubModel <- predict(bestSubModel, newdata = test) # validation predictions
# mean Prediction Error and Std Error - lower is better
meanPred <- mean((test$avgHrs - pred.bestSubModel )^2) # mean prediction error
stdError <- sd((test$avgHrs - pred.bestSubModel )^2)/length(test$avgHrs) # std error
# Append a row to our modelMetrics
modelMetrics[nrow(modelMetrics) + 1, ] <- c( "bestSubModel", 0.1789, aic, bic, meanPred, stdError)
modelMetrics
#standard model plots
par(mfrow = c(2, 2)) # Split the plotting panel into a 2 x 2 grid
plot(model.bkwd.lr1)
par(mfrow = c(1, 1))
#qq plot
qqnorm(resid(model.bkwd.lr1))
qqPlot(model.bkwd.lr1)
#residual plot
#Cooks Distance
cutoff <- 4/((nrow(train)-length(model.bkwd.lr1$coefficients)-2))
plot(model.bkwd.lr1, which=4, cook.levels=cutoff)
#visualize the model
visreg2d(model.bkwd.lr1, "numProj", "lastEvalLog", plot.type="persp" )
visreg2d(model.bkwd.lr1, "numProj", "lastEvalLog", plot.type="image" )
visreg(model.bkwd.lr1, "numProj")
#import in the data set
currentEmployees <- fread('https://raw.githubusercontent.com/lgellis/Predictive_Analytics_Tutorial/master/currentEmployees.csv')
currentEmployees
typeof(currentEmployees)
#convert to dataframe
currentEmployees <-as.data.frame(currentEmployees)
#Add the existing data transformations
colNames <- c("satLevel", "lastEval", "numProj", "timeCpny", "wrkAcdnt", "left", "fiveYrPrmo", "job", "salary")
setnames(currentEmployees, colNames)
attach(currentEmployees)
currentEmployees
#create factor variables with a better format
currentEmployees$leftFactor <- factor(left,levels=c(0,1),
labels=c("Did Not Leave Company","Left Company"))
currentEmployees$promoFactor <- factor(fiveYrPrmo,levels=c(0,1),
labels=c("Did Not Get Promoted","Did Get Promoted"))
currentEmployees$wrkAcdntFactor <- factor(wrkAcdnt,levels=c(0,1),
labels=c("No Accident","Accident"))
attach(currentEmployees)
currentEmployees <- cbind (currentEmployees, dummy(job), dummy(salary))
#Log Transforms for: satLevel, lastEval
currentEmployees$satLevelLog <- log(satLevel)
currentEmployees$lastEvalLog <- log(lastEval)
#SQRT Transforms for: satLevel, lastEval
currentEmployees$satLevelSqrt <- sqrt(satLevel)
currentEmployees$lastEvalSqrt <- sqrt(lastEval)
#Scale Transforms for: satLevel, lastEval
currentEmployees$satLevelScale <- scale(satLevel)
currentEmployees$lastEvalScale <- scale(lastEval)
currentEmployees$greatEvalLowSat <- ifelse(lastEval>0.8 & satLevel <0.2, 1, 0)
currentEmployees$Predictions <- predict(model.bkwd.lr1, currentEmployees) # test predictions
head(currentEmployees)
#Calculate Overage Hours
currentEmployees$predictedOverageHours = currentEmployees$Predictions - 160
#In case the employees worked less than 160 hours per month set negative values to 0
currentEmployees$netPredictedOverageHours <- currentEmployees$predictedOverageHours
currentEmployees$netPredictedOverageHours[currentEmployees$netPredictedOverageHours<0] <- 0
head(currentEmployees)
#Calculate total hours worked by each category
lowWageTotal = sum(currentEmployees[currentEmployees$salary=='low',]$predictedOverageHours)
mediumWageTotal = sum(currentEmployees[currentEmployees$salary=='medium',]$predictedOverageHours)
highWageTotal = sum(currentEmployees[currentEmployees$salary=='high',]$predictedOverageHours)
#assume low salary makes $10/hr, medium is $25/hr and high is $50/hr and assume they get time and a half during overtime
#total paid is sum of all predictions of low * 10, medium * 20, high * 40
TotalPaidPerMonth = lowWageTotal*10*1.5 + mediumWageTotal*25*1.5 + highWageTotal*50*1.5
TotalPaidPerMonth
TotalPaidPerYear = TotalPaidPerMonth * 12
TotalPaidPerYear