Variables that Predict Train Equipment Damage

1 Introduction

The goal of this study is to determine whether the number of Cars Derailed and the amount of Track Damage can be used to reliably predict Equipment Damage in U.S. railroad accidents. Specifically, I want to answer the research question:

Can the number of Cars Derailed and the amount of Track Damage be used to reliably predict Equipment Damage in train accidents?

To investigate this, I use accident data collected in 2010 from the Federal Railroad Administration (safetydata.fra.dot.gov). The dataset includes information such as accident cause, speed, injuries, fatalities, track type, and several measures of damage. From these variables, I focus on Cars Derailed and Track Damage because exploratory analysis shows they have the strongest correlation with Equipment Damage. I also create an interaction term (Cars Derailed × Track Damage) to test whether the effect of derailments depends on track conditions.

The approach will include exploratory analysis of the data set, construction of several candidate regression models (main effects, interaction model, and transformed model), residual diagnostics, and nonparametric techniques such as bootstrap confidence intervals. By comparing these models, I will identify which predictors meaningfully explain variation in Equipment Damage and evaluate the stability of the results.

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

The main finding of this analysis is that Equipment Damage in railroad accidents is overwhelmingly explained by the number of Cars Derailed. Both classical regression and bootstrap methods consistently showed that Cars Derailed is a strong and highly significant predictor, with narrow confidence intervals that exclude zero. In contrast, Track Damage alone and its interaction with Cars Derailed were not reliable predictors once Cars Derailed was included, as both residual and bootstrap intervals collapsed to zero. This indicates that derailments themselves, rather than track repair costs, are the primary drivers of equipment damage.

A few drawbacks limit this study. First, the data set comes from 2010, so the results may not fully generalize to today’s accidents given technological and safety improvements since then. Second, the residual analysis showed clear heteroscedasticity and outliers, which we attempted to address with a Box-Cox transformation, but that came at the cost of a lower R². Third, the model was limited to variables available in the dataset and did not include other potentially important predictors such as train speed, weather conditions, or train load.

There are several future improvements worth considering. Collecting more recent accident data would make the model more relevant to current conditions. Adding additional predictors such as weather, maintenance practices, or cargo weight may improve explanatory power. Exploring nonlinear models or robust regression could also help address heteroscedasticity and outliers more effectively.

In terms of applications, this model provides clear guidance for rail safety analysts and policymakers. Since the number of Cars Derailed is the most consistent predictor of costly equipment damage, mitigation strategies should focus on reducing derailments through improved inspection protocols, stronger track maintenance, and upgraded braking systems. While Track Damage itself was not a strong predictor in this analysis, it remains a real financial burden, and future models could integrate engineering assessments to capture its true costs more fully. Overall, this project demonstrates how statistical modeling can guide data-driven decision-making to reduce accident losses.

10 References

safetydata.fra.dot.gov