Multiple Linear Regression Models

1 Introduction

In this assignment, we are tasked with creating different multiple regression models so we can incorporate more variables into the model. I want to take a look and see which variables are most predictive, and which model is the best.

2 Description of Data Set

The data set I chose was Railroad Accidents in the U.S. in 2010. The data was collected via safetydata.fra.dot.gov. In last week’s project, I used this data set to determine if the number of Cars Derailed positively correlated with Equipment Damage (in dollars). This week, I want to add Track Damage (in dollars) to my model because it had the second strongest correlation coefficient r from last week’s correlation matrix. I’m going to see how accurately Cars Derailed and Track Damage can predict Equipment Damage.

2.1 Variables

• Railroad (categorical) - name of railroad
• Month (categorical)
• Day (categorical) - days of the week labeled 1 through 7
• State (categorical)
• County (categorical)
• Track Type (categorical)
• Track Management (categorical)
• Accident Type (categorical)
• Accident Cause (categorical)
• Equipment Damage (numerical) - amount of money spent on damages
• Track Damage (numerical) - amount of money spent on track repair
• Killed (numerical) - number of people killed from accident
• Injured (numerical) - number of people injured from accident
• RREquip (categorical) - type of train
• Speed (numerical) - mph the train was going at during impact
• Locomotives Derailed (numerical) - the number of locomotives derailed
• Cars Derailed (numerical) - the number of Cars Derailed

3 Methodology and Analysis

I am going to create two multiple linear regression models: one with Cars Derailed and Track Damage being added together, and one with the two variables being multiplied (interacted) together. There is a possibility that a single derailed car causes more damage if the track is already badly damaged. If that’s the case, then the interaction model would be more appropriate to use. We will be dropping all of the other variables mentioned in #2.

4 Model 1: Main Effects Model

We are going to start with the MLR model without the interaction term.

4.1 Inferential Statistics

# Multiple Linear Regression
model1 <- lm(EqpDamg ~ CarsDer + TrkDamg, data = myData)

# Summary of the model
summary(model1)

Call:
lm(formula = EqpDamg ~ CarsDer + TrkDamg, data = myData)

Residuals:
     Min       1Q   Median       3Q      Max 
-1019731   -48772     9400    32944  3143734 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.988e+04  4.001e+03  -4.968 7.19e-07 ***
CarsDer      2.538e+04  7.260e+02  34.961  < 2e-16 ***
TrkDamg      5.387e-01  3.819e-02  14.108  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 180500 on 2618 degrees of freedom
Multiple R-squared:  0.4733,    Adjusted R-squared:  0.4729 
F-statistic:  1176 on 2 and 2618 DF,  p-value: < 2.2e-16

Looking at the inferential statistics, Cars Derailed and Track Damage are both highly significant. The intercept is not statistically significant. For each additional car derailed (25,380, p < 2e-16), Equipment Damage increases by about $25,380, holding Track Damage constant. For each additional dollar spent on Track Damage (0.54, p < 2e-16), Equipment Damage increases by about $0.54, holding Cars Derailed constant. This model explains 47.29% of the variation in Equipment Damage. The model is highly significant (F = 1176, p = 0)

4.2 Residual Analysis

# Fit model without interaction
model1 <- lm(EqpDamg ~ CarsDer + TrkDamg, data = myData)
# Plot residuals vs fitted to check assumptions
par(mfrow = c(2,2))
plot(model1)

par(mfrow = c(1,1))

The Residuals vs Fitted plot shows there is no constant variance in the residuals. The Q-Q Residuals shows multiple outliers and two strong tails, indicating non-normality. The Scale Location plot shows a lack of constant variance and the Residuals vs Leverage plot shows an outlier at 390. The assumptions for normality are not met.

5 Model 2: Interaction Term

Now, I am going to incorporate the interaction term, CarsDer:TrkDamg, exclude the main effects, and see if the R-squared value increased and if the interaction term is statistically significant. I’ll also be looking to see if the model as a whole is significant. If it’s not, then I’ll be removing the term and go forward with the MLR model.

5.1 Inferential Statistics

# Multiple Linear Regression with Interaction
model2 <- lm(EqpDamg ~ CarsDer:TrkDamg, data = myData)

# Summary of the model
summary(model2)

Call:
lm(formula = EqpDamg ~ CarsDer:TrkDamg, data = myData)

Residuals:
     Min       1Q   Median       3Q      Max 
-1799092   -44483   -33634   -15771  3153755 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     4.624e+04  3.635e+03   12.72   <2e-16 ***
CarsDer:TrkDamg 7.002e-02  1.508e-03   46.42   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 184200 on 2619 degrees of freedom
Multiple R-squared:  0.4514,    Adjusted R-squared:  0.4512 
F-statistic:  2155 on 1 and 2619 DF,  p-value: < 2.2e-16

For each additional unit increase in CarsDer:TrkDamg, Equipment Damage increases by $0.07. The interaction term is highly significant (F = 2155, p = 0). This model explains 45.12% of the variation in Equipment Damage.

5.2 Residual Analysis

# Fit the interaction model
model2 <- lm(EqpDamg ~ CarsDer:TrkDamg, data = myData)

# Residual diagnostics
par(mfrow = c(2,2))  # show 4 plots together
plot(model2)

par(mfrow = c(1,1))  # reset layout

The Residuals vs Fitted plot shows there is no constant variance in the residuals. The Q-Q Residuals shows multiple outliers and one strong tail, showing non-normality. The Scale Location plot shows a lack of constant variance and the Residuals vs Leverage plot shows an outlier at 1414. All in all, the assumptions for normality are not met.

6 Model 3: Interaction Term with Main Effects

Finally, I will be using both main effects and their interaction term to create Model 3.

6.1 Inferential Statistics

# Model 3: CarsDer + TrkDamg + CarsDer:TrkDamg
model3 <- lm(EqpDamg ~ CarsDer + TrkDamg + CarsDer:TrkDamg, data = myData)

# Summary of the model
summary(model3)

Call:
lm(formula = EqpDamg ~ CarsDer + TrkDamg + CarsDer:TrkDamg, data = myData)

Residuals:
     Min       1Q   Median       3Q      Max 
-1084820   -37763    -3774    11799  3144550 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     3.401e+03  3.974e+03   0.856    0.392    
CarsDer         1.735e+04  8.114e+02  21.382   <2e-16 ***
TrkDamg         7.346e-02  4.398e-02   1.670    0.095 .  
CarsDer:TrkDamg 4.165e-02  2.268e-03  18.363   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 169900 on 2617 degrees of freedom
Multiple R-squared:  0.5334,    Adjusted R-squared:  0.5329 
F-statistic: 997.3 on 3 and 2617 DF,  p-value: < 2.2e-16

Cars Derailed (17,350, p < 2e-16), Track Damage (0.073, p = 0.095), and their interaction term (p = 0) are all statistically significant. Since the interaction term is significant, this means that the variables depend on each other. The effect of Cars Derailed on Equipment Damage grows as Track Damage grows. When both Cars Derailed and Track Damage are zero, the expected Equipment Damage is $3400. The intercept is not statistically significant. For each additional car derailed, Equipment Damage increases by about $17,350 on average, holding Track Damage constant at zero. Each additional dollar of Track Damage is associated with about $0.07 of additional Equipment Damage. This model explains 53.39% of the variation in Equipment Damage. The model is highly significant (F = 997.3, p = 0).

6.2 Residual Analysis

# Residual diagnostic plots
par(mfrow = c(2,2))
plot(model3)

par(mfrow = c(1,1))  # reset

The assumptions of normality are not met. There is no constant variance in the Residuals vs Fitted plot, there is a strong right tail in the Q-Q Residuals plot, the Scale Location plot does not show constant variance, and there is one outlier, 1414, in the Residuals vs Leverage plot.

7 ANOVA

anova(model1, model2, model3)

Analysis of Variance Table

Model 1: EqpDamg ~ CarsDer + TrkDamg
Model 2: EqpDamg ~ CarsDer:TrkDamg
Model 3: EqpDamg ~ CarsDer + TrkDamg + CarsDer:TrkDamg
  Res.Df        RSS Df   Sum of Sq      F    Pr(>F)    
1   2618 8.5269e+13                                    
2   2619 8.8817e+13 -1 -3.5477e+12 122.91 < 2.2e-16 ***
3   2617 7.5536e+13  2  1.3281e+13 230.06 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The ANOVA table shows that adding the interaction term to the main effects significantly improves the model. As a result, Model 3 is the best model because it reduces residual variation the most. Model 2 is the worst of the three models because it doesn’t take into account the influence that the independent main effects have on the model. Model 2 also has the highest residual sum of squares, meaning there is more unexplained variation in this model compared to the other two.

8 Nonparametric Tests

We are going to move forward with Model 3 because we established it as the best of the three. I will be conducting a Box-Cox transformation and Bootstrap confidence intervals because the assumptions for normality were not met.

8.1 Box-Cox Transformation

library(MASS)

# Shift the response to make all values positive
myData$EqpDamg_pos <- myData$EqpDamg + abs(min(myData$EqpDamg)) + 1

# Fit Model 3: CarsDer + TrkDamg + interaction
model3_pos <- lm(EqpDamg_pos ~ CarsDer + TrkDamg + CarsDer:TrkDamg, data = myData)

# Box-Cox transformation
bc <- boxcox(model3_pos,
             lambda = seq(-2, 2, 0.1),
             lambda.lab = "Lambda",
             ylab = "Log-Likelihood",
             main = "Box-Cox Transformation for Model 3")

# Identify lambda that maximizes log-likelihood
lambda_opt <- bc$x[which.max(bc$y)]
lambda_opt

[1] 0.2222222

I conducted a Box-Cox transformation to assess whether transforming Equipment Damage would improve the model assumptions. Since λ = 1 is within the confidence interval, no transformation is needed.

8.2 Bootstrap Confidence Intervals

library(boot)

# Define function to extract coefficients
boot_coef <- function(data, indices) {
  d <- data[indices, ]  # resample rows
  fit <- lm(EqpDamg ~ CarsDer + TrkDamg + CarsDer:TrkDamg, data = d)
  return(coef(fit))
}

# Run bootstrap
set.seed(123)  # for reproducibility
boot_results <- boot(data = myData, statistic = boot_coef, R = 1000)  # 1000 resamples

# Create a tidy table of 95% percentile CIs for all coefficients
coef_names <- names(coef(model3))
ci_list <- lapply(1:length(coef_names), function(i) {
  ci <- boot.ci(boot_results, type = "perc", index = i)$percent[4:5]
  data.frame(Coefficient = coef_names[i],
             Lower95 = ci[1],
             Upper95 = ci[2])
})

ci_table <- do.call(rbind, ci_list)
print(ci_table)

      Coefficient       Lower95      Upper95
1     (Intercept) -5542.0358228 1.329135e+04
2         CarsDer 12827.1587568 2.172305e+04
3         TrkDamg    -0.1240326 3.215833e-01
4 CarsDer:TrkDamg     0.0289688 6.334345e-02

The Cars Derailed interval is entirely positive. This indicates strong evidence of a positive effect. We are 95% confident that the slope of Cars Derailed is between 12,800 and 21,700. The interval for Track Damage includes 0, which means the effect alone may not have a consistent effect on Equipment Damage. We are 95% confident that the slope of Track Damage is between -0.124 and 0.322. The interval for the interaction term is entirely positive, meaning the interaction is meaningful. We are 95% confident that the slope of the interaction term is between 0.03 and 0.06.

9 Conclusion

In conclusion, the model with the interaction term and both main effects, Model 3, was the best model of the three. It explained the highest variation in Equipment Damage (53.29%). We are certain that the main effect of Cars Derailed is significant, but are unsure if Track Damage alone is significant because its 95% bootstrap interval contains 0. The interaction term is statistically significant, meaning that the effect of Cars Derailed depends on Track Damage and vice versa.

10 References

safetydata.fra.dot.gov