Simple Linear Regression

In building simple linear regression, it is important to have predictor variable that has high linear relationship with the target variables. Based on the correlation matrix, I will be using pol_exp as my predictor variable. Prior to model building, I will inspect if there is any outlier from our data.

# further EDA - outlier inspection
plot(crimex$pol_exp, crimex$crime_rate)

boxplot(crimex$pol_exp)

Based on the boxplot, there is one outlier in our pol_exp data. I will remove this outlier first before building our model because based on the scatterplot, it looks like it has a high influence on our model.

# remove outliers
crime.clean <- crimex %>% 
  filter(pol_exp<320)
         
boxplot(crime.clean$pol_exp)

plot(crime.clean$pol_exp, crime.clean$crime_rate)

# model building
slm.o <- lm(crime_rate ~ pol_exp, crimex) 
slm <- lm(crime_rate ~ pol_exp, crime.clean) # without outlier


# plotting the linear regression model
plot(crimex$pol_exp, crimex$crime_rate)
abline(slm.o$coefficients[1],slm.o$coefficients[2], col = "red")
abline(slm$coefficients[1],slm$coefficients[2], col = "green") # without outlier

summary(slm.o)

## 
## Call:
## lm(formula = crime_rate ~ pol_exp, data = crimex)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -581.07 -151.94   20.33  147.24  580.59 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 152.0674   128.5331   1.183    0.243    
## pol_exp       4.5573     0.7354   6.197 1.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 287.2 on 45 degrees of freedom
## Multiple R-squared:  0.4605, Adjusted R-squared:  0.4485 
## F-statistic:  38.4 on 1 and 45 DF,  p-value: 1.591e-07

summary(slm)

## 
## Call:
## lm(formula = crime_rate ~ pol_exp, data = crime.clean)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -650.8 -153.4   28.0  163.1  541.3 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  45.4783   130.6502   0.348    0.729    
## pol_exp       5.2941     0.7679   6.894 1.64e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.8 on 44 degrees of freedom
## Multiple R-squared:  0.5193, Adjusted R-squared:  0.5083 
## F-statistic: 47.53 on 1 and 44 DF,  p-value: 1.636e-08

From both of the regression model, pol_exp has positive corelation with crime.rate and is significantly affecting crime_rate (3 stars/p-value < 0.05).

If we look at the abline of the plot, we can see that the outliers of our data have a high influence over our model. Therefore, it is safer to remove the outliers from the data. The removal of outliers increase the R-squared of the model ~0.6 points. Based on this analysis, I will use slm as our simple linear model. Additionally, It also has higher R-squared as one of the tools for evaluating a linear regression model.

R-squared value tells us how well our model describes our data. It measures the extent to which the variance in our target variable (crime_rate) can be explained by the predictor variables (pol_exp, etc.). As we increase the number of predictor variables, our model’s R-squared value will also increase as it incorporates new legitimate information as well as the noise introduced by these extra variables. Because of this, it is better to evaluate the adjusted R-squared of our model. Adjusted R-squared does not increase the way regular multiple R-squared does because it is adjusted for the number of predictor variables in our model. It increases only when the new variable actually leads to a better prediction.

slm model can be written as:

crime_rate = 45.4783 + 5.2941(pol_exp)

which also means:

every 1 unit increase of police expenditure in the last 2 years, there will be an increase in crime rate at about 5.2941 times of the police expenditure.

Even so, our simple linear model only have an adjusted R-squared of 0.5083. An improvement can be made by adding more data/adding more variables to our model that may explain our target. On the constraint of only having this dataset, I will try to add more variables and use multiple linear regression for this case.

Multiple Linear Regression

On building multiple linear regression, I will use automatic feature selection using step-wise regression with backward elimination. I will use crime_rate as our target variable and all the remaining variables as our predictor, except the nonwhites_per1000 variable which has zero correlation with crime_rate.

mlm <- lm(crime_rate ~ percent_m + is_south + mean_education + labour_participation + 
    m_per1000f + state_pop + unemploy_m24 + 
    unemploy_m39 + gdp + inequality + prob_prison + time_prison + 
    pol_exp, crime.clean)
step(mlm, direction = "backward")

## Start:  AIC=499.41
## crime_rate ~ percent_m + is_south + mean_education + labour_participation + 
##     m_per1000f + state_pop + unemploy_m24 + unemploy_m39 + gdp + 
##     inequality + prob_prison + time_prison + pol_exp
## 
##                        Df Sum of Sq     RSS    AIC
## - state_pop             1        22 1298610 497.42
## - is_south              1       113 1298701 497.42
## - time_prison           1      3277 1301864 497.53
## - labour_participation  1      3529 1302116 497.54
## - gdp                   1      8271 1306859 497.71
## - unemploy_m24          1     36151 1334738 498.68
## - m_per1000f            1     42497 1341084 498.90
## <none>                              1298587 499.41
## - unemploy_m39          1     88024 1386611 500.43
## - prob_prison           1    117246 1415833 501.39
## - percent_m             1    168873 1467461 503.04
## - mean_education        1    299178 1597766 506.95
## - inequality            1    432350 1730938 510.63
## - pol_exp               1   1067853 2366441 525.02
## 
## Step:  AIC=497.42
## crime_rate ~ percent_m + is_south + mean_education + labour_participation + 
##     m_per1000f + unemploy_m24 + unemploy_m39 + gdp + inequality + 
##     prob_prison + time_prison + pol_exp
## 
##                        Df Sum of Sq     RSS    AIC
## - is_south              1       114 1298724 495.42
## - time_prison           1      3523 1302132 495.54
## - labour_participation  1      3546 1302155 495.54
## - gdp                   1      8633 1307242 495.72
## - unemploy_m24          1     36360 1334970 496.69
## - m_per1000f            1     48877 1347487 497.12
## <none>                              1298610 497.42
## - unemploy_m39          1     88625 1387234 498.45
## - prob_prison           1    117555 1416165 499.40
## - percent_m             1    168904 1467514 501.04
## - mean_education        1    299444 1598054 504.96
## - inequality            1    464523 1763132 509.48
## - pol_exp               1   1163183 2461793 524.84
## 
## Step:  AIC=495.42
## crime_rate ~ percent_m + mean_education + labour_participation + 
##     m_per1000f + unemploy_m24 + unemploy_m39 + gdp + inequality + 
##     prob_prison + time_prison + pol_exp
## 
##                        Df Sum of Sq     RSS    AIC
## - time_prison           1      3507 1302231 493.54
## - labour_participation  1      4958 1303683 493.59
## - gdp                   1      9106 1307830 493.74
## - unemploy_m24          1     41820 1340544 494.88
## - m_per1000f            1     49582 1348307 495.14
## <none>                              1298724 495.42
## - unemploy_m39          1     93362 1392087 496.61
## - prob_prison           1    128301 1427026 497.75
## - percent_m             1    180480 1479204 499.41
## - mean_education        1    299371 1598095 502.96
## - inequality            1    597093 1895817 510.82
## - pol_exp               1   1168618 2467342 522.94
## 
## Step:  AIC=493.54
## crime_rate ~ percent_m + mean_education + labour_participation + 
##     m_per1000f + unemploy_m24 + unemploy_m39 + gdp + inequality + 
##     prob_prison + pol_exp
## 
##                        Df Sum of Sq     RSS    AIC
## - labour_participation  1      3805 1306036 491.68
## - gdp                   1     10889 1313120 491.93
## - unemploy_m24          1     44275 1346506 493.08
## - m_per1000f            1     46320 1348551 493.15
## <none>                              1302231 493.54
## - unemploy_m39          1    102971 1405202 495.04
## - percent_m             1    213507 1515738 498.53
## - prob_prison           1    220749 1522981 498.75
## - mean_education        1    295885 1598116 500.96
## - inequality            1    610210 1912441 509.22
## - pol_exp               1   1172569 2474800 521.08
## 
## Step:  AIC=491.68
## crime_rate ~ percent_m + mean_education + m_per1000f + unemploy_m24 + 
##     unemploy_m39 + gdp + inequality + prob_prison + pol_exp
## 
##                  Df Sum of Sq     RSS    AIC
## - gdp             1     10398 1316434 490.04
## - unemploy_m24    1     40481 1346517 491.08
## - m_per1000f      1     44751 1350787 491.23
## <none>                        1306036 491.68
## - unemploy_m39    1    104836 1410871 493.23
## - prob_prison     1    218502 1524538 496.79
## - percent_m       1    231590 1537626 497.19
## - mean_education  1    300525 1606560 499.20
## - inequality      1    609144 1915180 507.29
## - pol_exp         1   1208349 2514385 519.81
## 
## Step:  AIC=490.04
## crime_rate ~ percent_m + mean_education + m_per1000f + unemploy_m24 + 
##     unemploy_m39 + inequality + prob_prison + pol_exp
## 
##                  Df Sum of Sq     RSS    AIC
## - unemploy_m24    1     47775 1364208 489.68
## - m_per1000f      1     50927 1367360 489.79
## <none>                        1316434 490.04
## - unemploy_m39    1    118175 1434609 492.00
## - percent_m       1    221285 1537719 495.19
## - prob_prison     1    257245 1573678 496.25
## - mean_education  1    327854 1644287 498.27
## - inequality      1    828730 2145163 510.50
## - pol_exp         1   1736089 3052523 526.73
## 
## Step:  AIC=489.68
## crime_rate ~ percent_m + mean_education + m_per1000f + unemploy_m39 + 
##     inequality + prob_prison + pol_exp
## 
##                  Df Sum of Sq     RSS    AIC
## - m_per1000f      1     18044 1382252 488.29
## <none>                        1364208 489.68
## - unemploy_m39    1     86771 1450980 490.52
## - percent_m       1    225573 1589781 494.72
## - prob_prison     1    262958 1627167 495.79
## - mean_education  1    297172 1661380 496.75
## - inequality      1    952361 2316569 512.04
## - pol_exp         1   3019469 4383677 541.38
## 
## Step:  AIC=488.29
## crime_rate ~ percent_m + mean_education + unemploy_m39 + inequality + 
##     prob_prison + pol_exp
## 
##                  Df Sum of Sq     RSS    AIC
## <none>                        1382252 488.29
## - unemploy_m39    1    115767 1498020 489.99
## - prob_prison     1    257895 1640148 494.16
## - percent_m       1    265460 1647713 494.37
## - mean_education  1    530001 1912253 501.22
## - inequality      1   1020942 2403194 511.73
## - pol_exp         1   3001462 4383714 539.38

## 
## Call:
## lm(formula = crime_rate ~ percent_m + mean_education + unemploy_m39 + 
##     inequality + prob_prison + pol_exp, data = crime.clean)
## 
## Coefficients:
##    (Intercept)       percent_m  mean_education    unemploy_m39  
##      -4605.850           8.724          16.721           7.064  
##     inequality     prob_prison         pol_exp  
##          7.051       -3863.787           6.796

This step-wise regression method will return the best formula that gives lowest AIC (Akaike Information Criterion), whereas the lower the AIC, the lower the information loss that our model has.

# assigning multiple linear regression model to an object
mlm.step <- lm(formula = crime_rate ~ percent_m + mean_education + unemploy_m39 + 
    inequality + prob_prison + pol_exp, data = crime.clean)

summary(mlm.step)

## 
## Call:
## lm(formula = crime_rate ~ percent_m + mean_education + unemploy_m39 + 
##     inequality + prob_prison + pol_exp, data = crime.clean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -575.90 -117.56   -1.68  125.80  448.44 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -4605.8496   860.8376  -5.350 4.12e-06 ***
## percent_m          8.7238     3.1876   2.737 0.009295 ** 
## mean_education    16.7212     4.3241   3.867 0.000407 ***
## unemploy_m39       7.0635     3.9083   1.807 0.078433 .  
## inequality         7.0513     1.3138   5.367 3.91e-06 ***
## prob_prison    -3863.7869  1432.3620  -2.697 0.010267 *  
## pol_exp            6.7957     0.7385   9.202 2.55e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 188.3 on 39 degrees of freedom
## Multiple R-squared:  0.7985, Adjusted R-squared:  0.7676 
## F-statistic: 25.77 on 6 and 39 DF,  p-value: 3.888e-12

If we compare between our simple linear regression and multiple linear regression model, we can see that our multiple linear regression gave better adjusted R-squared of 0.7676 from previously 0.5083. Therefore, I will use the mlm.step model for our prediction.

From this model, we can see that probability have negative correlation with crime_rate, while other variable positively correlate with crime_rate. The formula of our final model is:

crime_rate = 8.7238(percent_m) + 16.7212(mean_education) + 7.0635(unemploy_m39) + 7.0513(inequality) - 3863.7869(prob_prison) + 6.7957(pol_exp) - 4605.8496

Linearity

This assumption is already checked prior to model building using correlation matrix in EDA step, and through scatterplot between our predictor and our target variable.

Normality

hist(mlm.step$residuals, breaks = 5)

shapiro.test(mlm.step$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mlm.step$residuals
## W = 0.96654, p-value = 0.2051

P-value > 0.05; H0 accepted. This result means that the residuals of our model distributed normally, therefore our model will have error around its mean.

Heteroscedasticity

plot(crime.clean$crime_rate, mlm.step$residuals)
abline(h = 0, col = "red")

bptest(mlm.step)

## 
##  studentized Breusch-Pagan test
## 
## data:  mlm.step
## BP = 8.3082, df = 6, p-value = 0.2164

P-value > 0.05; H0 accepted. This means that our model residuals do not have a pattern that cannot be captured by our model.

Variance Inflation Factor (Multicollinearity)

vif(mlm.step)

##      percent_m mean_education   unemploy_m39     inequality    prob_prison 
##       1.970723       3.035998       1.412245       3.486175       1.343397 
##        pol_exp 
##       1.955926

There are no values equal to or more than 10 so Multicollinearity between variables is not found (predictor variables are mutually independent).

Based on the analysis, our mlm.step model passed all of the assumption evaluation.