options(scipen=999)
library(socviz) # use for the round_df() function
library(datasetsICR)
#load("Catholic_viz.RData")
#load("NBA_viz.RData")
load("Murders96_viz.RData")DAT-4313 - Data Viz in Model Development
Murder Data
FURTHER INSTRUCTIONS
Use the dataset assigned to you in the LAB Instructions.
Use the depvar indicated for your dataset.
Select 5-8 independent variables from your dataset as predictors, and estimate a linear regression model.
Note that not all variables in your dataset make sense as predictors. For example, categorical variables with many levels or all components of a derived variable are unreasonable.
Work through the EDA to identified good candidates. Your model should be linear in the parameters. Consider alternative functional forms for your depvar or continuous variables.
Include interpretations and comments after each chunk.
Conclude with a final assessment of the model and a reflection of the value of data visualization to aid model development.
Code for Murder96 with depvar=arrests
DATA
This code loads the data in.
Partition 60% Train / 40% Test
library(caret)
set.seed(1)
index <- createDataPartition(Murders96_viz$arrests, p=0.6, list=FALSE)
train <- Murders96_viz[index,]
test <- Murders96_viz[-index,]This code allows us to partition the data into a training set and a test set.
Training set number of observations: 1320
Test dataset number of observations: 877
EDA
Descriptive Statistics
library(psych)
library(tidyverse)
library(flextable)
psych::describe(train, fast=TRUE) vars n mean sd min max range se
arrests 1 1320 5.21 18.76 0.00 320.00 320.00 0.52
countyid 2 1320 32550.57 15488.80 1001.00 56043.00 55042.00 426.32
density 3 1320 265.34 1817.83 0.08 54058.77 54058.69 50.03
popul 4 1320 94031.70 223537.24 141.00 3126966.00 3126825.00 6152.66
perc1019 5 1320 15.11 1.88 7.60 24.54 16.94 0.05
perc2029 6 1320 12.05 3.19 5.62 34.75 29.14 0.09
percblack 7 1320 8.34 13.93 0.00 86.28 86.28 0.38
percmale 8 1320 49.33 1.63 44.83 61.64 16.82 0.04
rpcincmaint 9 1320 207.78 116.83 19.14 1281.10 1261.96 3.22
rpcpersinc 10 1320 12569.38 3031.90 5322.20 41094.22 35772.02 83.45
rpcunemins 11 1320 52.90 34.85 0.00 287.10 287.10 0.96
year 12 1320 1996.00 0.00 1996.00 1996.00 0.00 0.00
murders 13 1320 5.88 26.27 0.00 420.00 420.00 0.72
murdrate 14 1320 0.42 0.77 0.00 11.01 11.01 0.02
arrestrate 15 1320 0.48 0.94 0.00 13.99 13.99 0.03
statefips 16 1320 32.45 15.46 1.00 56.00 55.00 0.43
countyfips 17 1320 99.81 109.51 1.00 840.00 839.00 3.01
execs 18 1320 0.01 0.14 0.00 3.00 3.00 0.00
lpopul 19 1320 10.40 1.35 4.95 14.96 10.01 0.04
execrate 20 1320 0.00 0.04 0.00 1.18 1.18 0.00
statename 21 1320 NaN NA Inf -Inf -Inf NA
statecode 22 1320 NaN NA Inf -Inf -Inf NA
region 23 1320 NaN NA Inf -Inf -Inf NAAnalyzing the descriptive statistics for each variable allows us to see the different aspects of each variable individually, but it does not really teach anything about the categorical variables.
Boxplots – All Numeric
data_long <- train %>%
select(where(is.numeric)) %>%
pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value")
ggplot(data_long, aes(x = Variable, y = Value)) +
geom_boxplot() +
coord_flip() +
facet_wrap(~ Variable, scales = "free", ncol=3) +
theme_minimal() +
theme(axis.text.x = element_blank()) +
labs(title = "Horizontal Boxplot for Each Numeric Variable")
The boxplots are very insightful, and can teach us some interesting information. First, arrests are typically very low, but there are a lot of outliers. Murders is pretty much the same as arrests. It has a very low mean, but there are a lot of outliers.
library(DataExplorer)
plot_boxplot(train, by="region") 
By keeping the variables on the x axis, but putting the regions on the y axis, each variable can be be seen split up by the different regions. By looking extremely closely, arrests and murder rate both have higher means in the Northeast. Overall, the results are fairly even for each region.
Histograms
plot_histogram(train)
The histograms show that arrests, murders, and executions are extremely uncommon relative to the amount of people.
Scatterplots (depvar ~ all x)
Look at the shape of the relationship between the dependent variable and all of the continuous potential independent variables.
plot_scatterplot(test, by="arrests")
By viewing the scatterplots of each variable by the arrests variable, there can be some observations. It appears that the places that had a high amount of arrests also had a fairly high population. This means that the amount of arrests relative to the populations is pretty even. It can also be seen that the amount of arrests goes down in the places that have had executions (excluding an outlier).
Correlation Matrix
train <- train[, !(names(train) == "year")]
plot_correlation(train, type="continuous")
The correlation matrix gives some super helpful information. It shows strong correlation between arrest and density, arrests and popul, arrests and murders, and arrests and lpopul. This will be helpful in creating the linear regression model as these variables can be used as predictors along with some other fairly high correlations.
MODEL
Linear Regression Model
Estimate Coefficients and show coefficients table
model <- lm((arrests) ~ lpopul + execs + arrestrate + murdrate + murders + rpcpersinc + popul + density, data=train)
model %>% as_flextable()Estimate | Standard Error | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|
(Intercept) | -3.303 | 2.338 | -1.413 | 0.1579 |
|
lpopul | 0.827 | 0.236 | 3.506 | 0.0005 | *** |
execs | 4.342 | 1.719 | 2.526 | 0.0117 | * |
arrestrate | 4.614 | 0.317 | 14.578 | 0.0000 | *** |
murdrate | -3.615 | 0.422 | -8.561 | 0.0000 | *** |
murders | 0.445 | 0.014 | 31.261 | 0.0000 | *** |
rpcpersinc | -0.000 | 0.000 | -4.858 | 0.0000 | *** |
popul | 0.000 | 0.000 | 10.439 | 0.0000 | *** |
density | 0.002 | 0.000 | 12.098 | 0.0000 | *** |
Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05 | |||||
Residual standard error: 8.232 on 1311 degrees of freedom | |||||
Multiple R-squared: 0.8086, Adjusted R-squared: 0.8074 | |||||
F-statistic: 692.2 on 1311 and 8 DF, p-value: 0.0000 | |||||
I used a model with 8 predictors. This allowed my model to explain 80.86% of the variability in arrests. The predictors that I selected are all significant at a confidence level of 95%, and all but one are significant at a confidence level of 99.9%.
Coefficient Magnitude Plot
# note this simple code using coefplot package produces chart like Figure 6.6 in one short line.
library(coefplot)
coefplot(model, soft="magnitude", intercept=FALSE)
This visualization shows the magnitude of the coefficient predictors. It appears that arrestrate and execs have the biggest impact.
Check for predictor independence
Using Variance Inflation Factors (VIF)
library(car)
vif(model) lpopul execs arrestrate murdrate murders rpcpersinc popul
1.976168 1.080618 1.709901 2.080418 2.720242 1.553693 3.127875
density
1.404460 The VIF shows if we need to worry about multicollinearity in the model. While there are no values over 5, the popul variable is fairly high, which makes sense as lpopul is just the log of popul. The same goes for murders and murdrate.
Residual Analysis
Residual Range
quantile(round(residuals(model, type = "deviance"),2)) 0% 25% 50% 75% 100%
-173.940 -1.340 -0.165 1.040 132.210 The residual range of our model is very good in the interquartile range. Outside of that, it might begin to be a problem with such high errors.
Residual Plots
We are looking for: - Random distribution of residuals vs fitted values - Normally distributed residuals : Normal Q-Q plot with values along line - Homoskedasticity with a Scale-Location line that is horizontal and no residual pattern - Minimal influential obs - that is, those outside the borders of Cook’s distance
plot(model)
The first plot does not appear to have a completely random distribution pattern for the residuals. The second plot does appear to be normally distributed as it follows the line fairly closely. The third plot has a fairly horizontal line and I cannot see a residual pattern. The fourth plot shows very few values outside of Cook’s distance which means there are minimal influential observations.
Plot Fitted Value by Actual Value
library(broom)
train2 <- augment(model, data=train) # this appends predicted to original dataset to build plots on your own
p <- ggplot(train2, mapping = aes(y=.fitted, x=arrests))
p + geom_point()
When fitted by arrests, it appears to be pretty linear.
Plot Residuals by Fitted Values
p <- ggplot(train2, mapping = aes(y=.resid, x=.fitted))
p + geom_point()
When fitted by residuals, it appears that it starts out linear, but then breaks off.
Performance Evaluation
Use Model to Score test dataset (Display First 10 values - depvar and fitted values only)
pred <- predict(model, newdata = test, interval="predict") #score test dataset wtih model
test_w_predict <- cbind(test, pred) # append score to original test dataset
test_w_predict %>%
head(10) %>%
select(countyid, arrests, fit, lwr,upr) %>%
as_flextable(show_coltype = FALSE) countyid | arrests | fit | lwr | upr |
|---|---|---|---|---|
1,003 | 6 | 5.8 | -10.3 | 22.0 |
1,005 | 1 | 1.6 | -14.5 | 17.8 |
1,009 | 0 | 0.7 | -15.5 | 16.8 |
1,015 | 20 | 13.5 | -2.7 | 29.6 |
1,019 | 2 | 5.3 | -10.8 | 21.5 |
1,021 | 0 | 1.2 | -14.9 | 17.4 |
1,029 | 2 | 6.9 | -9.3 | 23.0 |
1,031 | 8 | 7.1 | -9.1 | 23.3 |
1,041 | 1 | 3.4 | -12.7 | 19.6 |
1,043 | 5 | 5.4 | -10.8 | 21.6 |
n: 10 | ||||
We can now implement the test dataset and score it. It shows the arrests along with the fit.
Plot Actual vs Fitted (test)
plot(test_w_predict$arrests, test_w_predict$fit, pch=16, cex=1)
By comparing the actual vs the fitted, we are able to see a linear relationship, which means that it is predicting fairly accurately.
Performance Metrics
library(Metrics)
metric_label <- c("MAE","RMSE", "MAPE")
metrics <- c(round(mae(test_w_predict$arrests, test_w_predict$fit),4),
round(rmse(test_w_predict$arrests, test_w_predict$fit),4),
round(mape(test_w_predict$arrests, test_w_predict$fit),4))
pmtable <- data.frame(Metric=metric_label, Value = metrics)
flextable(pmtable)Metric | Value |
|---|---|
MAE | 3.8583 |
RMSE | 29.3548 |
MAPE | Inf |
-MAE (Mean Absolute Error): This metric indicates the average absolute difference between the predicted values and the actual values. In this case, an MAE of 3.8583 suggests that, on average, the model’s predictions deviate from the actual values by approximately 3.8583 units.
-RMSE (Root Mean Squared Error): RMSE is another measure of the differences between predicted and actual values. It’s calculated by taking the square root of the average of the squared differences between the predicted and actual values. An RMSE of 29.3548 indicates the square root of the average squared differences, which is approximately 29.3548 units. RMSE tends to penalize larger errors more heavily than MAE.
-MAPE (Mean Absolute Percentage Error): This metric calculates the average of the absolute percentage differences between predicted and actual values. An MAPE of “Inf” (infinity) typically means that there are instances where the actual values are zero or very close to zero, causing division by zero errors when calculating percentages. This could suggest that the model is inaccurately predicting values close to zero or that there’s an issue with the data. Alternatively, it could mean that the model is performing well for most values but fails to handle zero values effectively. This makes sense for our model as there are a lot of arrest variables that are zero.
Used ChatGPT to help with the interpretations.
Model Fit by Region
# more in depth plots of performance
p <- ggplot(data = subset(test_w_predict, region %in% c("West", "Northeast")),
aes(x = arrests,
y = fit, ymin = lwr, ymax = upr,
color = region,
fill = region,
group = region))
p + geom_point(alpha = 0.5) +
geom_line() + geom_ribbon(alpha = 0.2, color = FALSE) +
labs(title="Actual vs Fitted by Region",
subtitle="West and Northeast",
caption="Murders96{datasetsICR}") +
xlab(NULL) + ylab("fitted") +
theme(legend.position = "bottom")
Comparing the actual and fitted values for the regions Northeast and West, we can see a very similar pattern.
SUMMARY ASSESSMENT AND EVALUATION OF THE MODEL
The progression of visualization in my model helped me understand what was happening step by step. In the end I was able see why the results showed what they did. With my dataset (Murders96), it was difficult to understand or see relationship with the beginning visualizations, but that is not always the case for every dataset. The correlation matrix was where I was able to start seeing and understanding what was happening and what was important. In the end, I was pretty happy with how my predictions turned out given my confusion with the beginning interpretations.