lab6

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 Murder Data

DATA

options(scipen=999)

library(socviz) # use for the round_df() function

#load("FIFA_viz.RData")
#load("Catholic_viz.RData")
#load("NBA_viz.RData")
load("Murders96_viz.RData")
data <- Murders96_viz

Partition 60% Train / 40% Test

library(caret)
Warning: package 'caret' was built under R version 4.1.2
Loading required package: ggplot2
Loading required package: lattice
set.seed(1)
index <- createDataPartition(Murders96_viz$arrests, p=0.6, list=FALSE)
train <- Murders96_viz[index,]
test  <- Murders96_viz[-index,]

Descriptive Statistics

library(psych)

Attaching package: 'psych'
The following objects are masked from 'package:ggplot2':

    %+%, alpha
library(tidyverse)
Warning: package 'tidyverse' was built under R version 4.1.2
Warning: package 'tibble' was built under R version 4.1.2
Warning: package 'forcats' was built under R version 4.1.2
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ psych::%+%()    masks ggplot2::%+%()
✖ psych::alpha()  masks ggplot2::alpha()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ purrr::lift()   masks caret::lift()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(flextable)

Attaching package: 'flextable'

The following object is masked from 'package:purrr':

    compose
psych::describe(train, fast=TRUE)  # (class note --> remove flextable)
Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf
Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf
Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf
            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      NA

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")

library(DataExplorer)
plot_boxplot(train, by="region") 

plot_histogram(train)

plot_scatterplot(test, by="arrests")

Correlation Matrix

plot_correlation(train, type="continuous")
Warning in cor(x = structure(list(arrests = c(8L, 1L, 2L, 4L, 1L, 1L, 2L, : the
standard deviation is zero
Warning: Removed 38 rows containing missing values (`geom_text()`).

MODEL

Linear Regression Model

From the above EDA, we chose the following variables to start. Note that, given the shape of the relationship between arrests and density.
From experience and prior research, it is common to specify a dependent variable that is a currency variable (e.g., sales, revenue, wages) in log form.

Estimate Coefficients and show coefficients table

train <- na.omit(train[, c("arrests", "density", "popul", "perc1019", "perc2029", 
                            "percblack", "percmale", "rpcincmaint", "rpcpersinc", 
                            "rpcunemins", "murders", "murdrate", "arrestrate")])

model <- lm((arrests) ~ density + popul + perc1019 + perc2029 + percblack + percmale + 
              rpcincmaint + rpcpersinc + rpcunemins + murders + murdrate + arrestrate  
              , data = train)

# Display summary of the model
summary(model)

Call:
lm(formula = (arrests) ~ density + popul + perc1019 + perc2029 + 
    percblack + percmale + rpcincmaint + rpcpersinc + rpcunemins + 
    murders + murdrate + arrestrate, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-174.856   -1.036   -0.154    0.715  132.276 

Coefficients:
               Estimate  Std. Error t value             Pr(>|t|)    
(Intercept)  4.86227099  8.25435681   0.589              0.55593    
density      0.00170543  0.00015264  11.173 < 0.0000000000000002 ***
popul        0.00002176  0.00000162  13.431 < 0.0000000000000002 ***
perc1019     0.16763131  0.13710881   1.223              0.22170    
perc2029     0.01061716  0.08098726   0.131              0.89572    
percblack   -0.00255221  0.02142381  -0.119              0.90519    
percmale    -0.08570352  0.15279780  -0.561              0.57497    
rpcincmaint  0.00136038  0.00269721   0.504              0.61409    
rpcpersinc  -0.00030567  0.00011005  -2.778              0.00555 ** 
rpcunemins  -0.00641095  0.00696006  -0.921              0.35716    
murders      0.44833356  0.01417638  31.625 < 0.0000000000000002 ***
murdrate    -3.52570980  0.43171133  -8.167 0.000000000000000739 ***
arrestrate   4.55892299  0.32676661  13.952 < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.294 on 1307 degrees of freedom
Multiple R-squared:  0.8063,    Adjusted R-squared:  0.8045 
F-statistic: 453.4 on 12 and 1307 DF,  p-value: < 0.00000000000000022
# Show coefficients table
library(flextable)
model %>% as_flextable

Estimate

Standard Error

t value

Pr(>|t|)

(Intercept)

4.862

8.254

0.589

0.5559

density

0.002

0.000

11.173

0.0000

***

popul

0.000

0.000

13.431

0.0000

***

perc1019

0.168

0.137

1.223

0.2217

perc2029

0.011

0.081

0.131

0.8957

percblack

-0.003

0.021

-0.119

0.9052

percmale

-0.086

0.153

-0.561

0.5750

rpcincmaint

0.001

0.003

0.504

0.6141

rpcpersinc

-0.000

0.000

-2.778

0.0056

**

rpcunemins

-0.006

0.007

-0.921

0.3572

murders

0.448

0.014

31.625

0.0000

***

murdrate

-3.526

0.432

-8.167

0.0000

***

arrestrate

4.559

0.327

13.952

0.0000

***

Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

Residual standard error: 8.294 on 1307 degrees of freedom

Multiple R-squared: 0.8063, Adjusted R-squared: 0.8045

F-statistic: 453.4 on 1307 and 12 DF, p-value: 0.0000

Coefficient Magnitude Plot

# note this simple code using coefplot package produces chart like Figure 6.6 in one short line. 
library(coefplot)
Warning: package 'coefplot' was built under R version 4.1.2
coefplot(model, soft="magnitude", intercept=FALSE)

Check for predictor independence

Using Variance Inflation Factors (VIF)

library(car)
Warning: package 'car' was built under R version 4.1.2
Loading required package: carData
Warning: package 'carData' was built under R version 4.1.2

Attaching package: 'car'
The following object is masked from 'package:dplyr':

    recode
The following object is masked from 'package:purrr':

    some
The following object is masked from 'package:psych':

    logit
vif(model)
    density       popul    perc1019    perc2029   percblack    percmale 
   1.476366    2.515093    1.275677    1.277430    1.707980    1.186051 
rpcincmaint  rpcpersinc  rpcunemins     murders    murdrate  arrestrate 
   1.903997    2.134696    1.128059    2.659054    2.141794    1.795401

Residual Analysis

Residual Range

quantile(round(residuals(model, type = "deviance"),2))
       0%       25%       50%       75%      100% 
-174.8600   -1.0325   -0.1550    0.7125  132.2800

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)

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()

Plot Residuals by Fitted Values

p <- ggplot(train2, mapping = aes(y=.resid, x=.fitted))
p  + geom_point()

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( arrests, fit, lwr, upr, countyid, density, perc1019, perc2029, percblack) %>% 
  as_flextable(show_coltype = FALSE)

arrests

fit

lwr

upr

countyid

density

perc1019

perc2029

percblack

6

5.0

-11.3

21.2

1,003

77.1

13.9

11.6

13.5

1

1.5

-14.8

17.9

1,005

29.9

15.1

13.7

46.2

0

0.2

-16.1

16.5

1,009

67.2

14.2

13.0

1.4

20

12.7

-3.6

29.0

1,015

186.5

14.7

14.3

19.1

2

5.1

-11.2

21.4

1,019

38.3

14.1

12.1

7.0

0

0.9

-15.4

17.2

1,021

50.9

14.8

12.6

11.9

2

7.1

-9.2

23.4

1,029

24.0

14.3

13.3

4.9

8

6.8

-9.5

23.1

1,031

61.7

14.3

13.8

18.2

1

3.9

-12.4

20.2

1,041

22.2

14.8

11.8

27.2

5

4.5

-11.8

20.8

1,043

99.2

13.8

12.8

0.9

n: 10

Plot Actual vs Fitted (test)

plot(test_w_predict$arrests, test_w_predict$fit, pch=16, cex=1)

Performance Metrics

library(Metrics)

Attaching package: 'Metrics'
The following objects are masked from 'package:caret':

    precision, recall
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.5386

RMSE

28.9149

MAPE

Inf

Model Fit by Region

# more in depth plots of performance 
p <- ggplot(data = subset(test_w_predict, region %in% c("South", "Midwest")),
            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 with Upper and Lower CI",
              subtitle="South and Midwest Region",
              caption="Countymurders {wooldridge}") +
       xlab(NULL) + ylab("fitted") +
       theme(legend.position = "bottom")

ggplot2 explore (Intro Section)

This is just replication of the first section within Chapter 6 (as FYI)

#pip <- lm(log(arrests) ~ density, data=train)
#summary(pip)

# knitr does not like autoplot
#autoplot(lm(log(Wage.n) ~ Potential, data=train), data=train, colour = 'Nationality_Continent')


library(ggplot2)
p <- ggplot(data=train, mapping=aes(y=log(arrests), x=density))

p +  geom_point(alpha = 0.2) + 
       geom_smooth(method = "lm", aes(color = "OLS", fill = "OLS"))
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 623 rows containing non-finite values (`stat_smooth()`).

p + geom_point(alpha=0.1) +
    geom_smooth(color = "tomato", fill="tomato", method = MASS::rlm) +
    geom_smooth(color = "steelblue", fill="steelblue", method = "lm")
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 623 rows containing non-finite values (`stat_smooth()`).
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 623 rows containing non-finite values (`stat_smooth()`).

p + geom_point(alpha=0.1) +
    geom_smooth(color = "tomato", method = "lm", linewidth = 1.2, 
                formula = y ~ splines::bs(x, 3), se = FALSE)
Warning: Removed 623 rows containing non-finite values (`stat_smooth()`).

p + geom_point(alpha=0.1) +
    geom_quantile(color = "tomato", size = 1.2, method = "rqss",
                  lambda = 1, quantiles = c(0.20, 0.5, 0.85))
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Warning: Removed 623 rows containing non-finite values (`stat_quantile()`).
Smoothing formula not specified. Using: y ~ qss(x, lambda = 1)

SUMMARY ASSESSMENT AND EVALUATION OF THE MODEL

Visualizations are valuable for understanding and evaluating how well the linear regression model works. They help us see patterns and trends in the data more easily. For example, boxplots, histograms, and scatterplots show us how the data is spread out and if there are any unusual points. We can also use these plots to pick out which features might be most important for our model. The correlation matrix tells us how strongly different features are related to each other, which helps us see if some features might be too similar and cause problems for our model.

Looking at the residuals, or the differences between the predicted values and the actual values, helps us see if our model is doing a good job of predicting outcomes. If the residuals have a pattern, it means our model might not be capturing everything in the data. By looking at the fitted versus actual values and the residuals versus fitted values, we can see if there are any patterns or if our model is working well.

In conclusion, the MAE and RMSE values seem relatively low, indicating that our model’s predictions are reasonably close to the actual values on average. However, the MAPE value being infinity could be a cause for concern, as it suggests potential issues with the model’s performance, especially when with instances where the actual value is zero. Further investigation and refinement of the model may be needed to address this issue and improve its reliability.

END