options(scipen=999)
library(wooldridge)
library(socviz)
load("Murders96_viz.RData")LAB 6
Example Code for FIFA with depvar=Wage
DATA
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,]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) # (class note --> remove flextable) 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 NAInterpretation
Arrests, the dependent variable, has a mean of 5.21 and a standard deviation of 18.76. The variables statename, statecode, and region all have missing values. ### 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")
Interpretation
There are many boxplots in this chart. Many if them of outliers and a small difference between Q1 and Q3.
library(DataExplorer)
plot_boxplot(train, by="region") 
Interpretation
Including a region section to the boxplots gives us a more indepth look at these variables. Again, these boxplots have mutliple outliers, and a smaller range between Q1 and Q3.
Histograms
plot_histogram(train)
Interpretation of Histograms
The histograms above generally follow a bell-curve shape, or are skewed to the right. Percmale has a bell-curve shape and countyfips is skewed to the right.
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")
Interpretation
From these charts, I have decided to choose density, percmale, popul, percblack, percmale,and murdrate as the independent variables in the linear regression.
Correlation Matrix
plot_correlation(train, type="continuous")
| #### Interpretation This heat map has a large gray chunk in it because of the variable year. In this dataset only the year 1996 is being analyzed. |
| ## MODEL |
| ### Linear Regression Model |
| ### Estimate Coefficients and show coefficients table |
| ::: {.cell} |
{.r .cell-code} model <- lm(arrests ~ murdrate + popul + percblack + percmale + density, data=train) model |
| ::: {.cell-output .cell-output-stdout} ``` |
| Call: lm(formula = arrests ~ murdrate + popul + percblack + percmale + density, data = train) |
| Coefficients: (Intercept) murdrate popul percblack percmale density -4.64932537 4.23264983 0.00004872 0.10248119 0.03834277 0.00283916 ``` ::: |
{.r .cell-code} model %>% as_flextable |
| ::: {.cell-output-display} |
| ::: {.cell} |
{.r .cell-code} # 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) |
::: {.cell-output-display}  ::: ::: #### INterpretation This chart, which uses OLS, affirms what was interpreted above. The lenght of the indicates the magnitude of the coefficent. Murdrate and percmale had the highest positive and negative correlation, respectively. |
| ### Check for predictor independence |
| Using Variance Inflation Factors (VIF) |
| ::: {.cell} |
{.r .cell-code} library(car) vif(model) |
| ::: {.cell-output .cell-output-stdout} |
| ### Residual Analysis |
| #### Residual Range |
| ::: {.cell} |
{.r .cell-code} quantile(round(residuals(model, type = "deviance"),2)) |
| ::: {.cell-output .cell-output-stdout} |
| #### 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 |
| ::: {.cell} |
{.r .cell-code} plot(model) |
::: {.cell-output-display}  ::: |
::: {.cell-output-display}  ::: |
::: {.cell-output-display}  ::: |
::: {.cell-output-display}  ::: ::: |
| Plot Fitted Value by Actual Value |
| #### Interpretation The Residuals versus Fitted plot, is quite clustered together. There is similar variability across the different levels of the fitted values. The second graph presented above, displays that the variables are normally distributed.The points are clustered again in this third graph as well, indicating that this model’s errors might not have a consistent spread. Lastly, we see that there is minimal influential observations in this final chart. ::: {.cell} |
| ```{.r .cell-code} 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() ``` |
::: {.cell-output-display}  ::: ::: |
| Plot Residuals by Fitted Values |
| ::: {.cell} |
{.r .cell-code} p <- ggplot(train2, mapping = aes(y=.resid, x=.fitted)) p + geom_point() |
::: {.cell-output-display}  ::: ::: |
| ### Performance Evaluation |
#### Use Model to Score test dataset (Display First 10 values - depvar and fitted values only) |
| ::: {.cell} |
| ```{.r .cell-code} 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(murdrate, arrests, popul, percblack, percmale, arrestrate, density, fit) %>% as_flextable(show_coltype = FALSE) ``` |
| ::: {.cell-output-display} |
#### Plot Actual vs Fitted (test) |
| ::: {.cell} |
{.r .cell-code} plot(test_w_predict$arrests, test_w_predict$fit, pch=16, cex=1) |
::: {.cell-output-display}  ::: ::: |
| #### Performance Metrics |
| ::: {.cell} |
{.r .cell-code} 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) |
| ::: {.cell-output-display} |
| #### Interpretation The MAE value is 5.6806, suggesting the model’s predictions deviate from the actual values by 5.68 units. RMSE deviates from the actual values by approximately 39.74 units. The MAPE does not a value because |
| ::: {.cell} |
{.r .cell-code} # more in depth plots of performance p <- ggplot(data = subset(test_w_predict, region %in% c("AF", "SA")), aes(x = arrests, y = murdrate, ymin = percmale, ymax = percblack, color = arrests, fill = arrests, group = arrests)) p + geom_point(alpha = 0.5) + geom_line() + geom_ribbon(alpha = 0.2, color = FALSE) + labs(title="Actual vs Fitted with Upper and Lower CI", subtitle="Africa and South America", caption="{datasetsICR}") + xlab(NULL) + ylab("fitted") + theme(legend.position = "bottom") |
::: {.cell-output-display}  ::: |
| ```{.r .cell-code} library(ggplot2) |
| coefficients <- coef(model) |
| upper_bound <- coefficients[“murdrate”] + coefficients[“percblack”] lower_bound <- coefficients[“murdrate”] - coefficients[“percmale”] |
| # Create the plot ggplot(data = train, aes(x = arrests, y = murdrate)) + geom_point(alpha = 0.5) + geom_line(color = “blue”) + geom_ribbon(aes(ymin = lower_bound, ymax = upper_bound), alpha = 0.2, color = “grey”) + labs(title = “Actual vs Fitted with arrests and murdrate”, subtitle = “There is significant variability in the points”, x = NULL, y = “Fitted”) + theme(legend.position = “bottom”) ``` |
::: {.cell-output-display}  ::: ::: |
Interpretation
There are two peaks in the fitted line, one around values slightly above 100 and another around values above 200 on the x-axis. These peaks might correspond to distinct patterns or underlying factors related to arrests and murder rates.
ggplot2 explore (Intro Section)
These graphs here show that as the murder rate increases, the number of arrests tends to increase. The last graph shows an interesting peak at about 7, and then a decrease after that. This means when the murder rate is around 7, more arrests tend to occur.
SUMMARY ASSESSMENT AND EVALUATION OF THE MODEL
It is important to be able to visualize the linear regression model because just interpreting the numbers and looking at charts together creates a coherent message. This helps when I would present information to a client or boss, because I would like my information to be as reliable and accurate as possible. Visualization enhances my understanding of the regression model, and it identifies areas of improvement. According to the results, I would have said the linear regression was decent at predicting the outcome of arrests, but after analyzing the visuals, I would say that I needed to edit some of my variables in the code to make it more easily understood.