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.

This data set has 17 variables:

1. Railroad (categorical) - name of railroad
2. Month (categorical)
3. Day (categorical) - days of the week labeled 1 through 7
4. State (categorical)
5. County (categorical)
6. Track Type (categorical)
7. Track Management (categorical)
8. Accident Type (categorical)
9. Accident Cause (categorical)
10. Equipment Damage (numerical) - amount of money spent on damages
11. Track Damage (numerical) - amount of money spent on track repair
12. Killed (numerical) - number of people killed from accident
13. Injured (numerical) - number of people injured from accident
14. Railroad Equipment (categorical) - type of equipment
15. Speed (numerical) - mph the train was going at during impact
16. Locomotives Derailed (numerical) - the number of locomotives derailed
17. Cars Derailed (numerical) - the number of Cars Derailed

3 Testing the Models

We are 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 because they are irrelevant in this study.

3.1 Main Effects Model

The first model we are going to look at is the one with just the main effects.

3.1.2 Inferential Statistics

First, we are going to find the inferential statistics for the main effects model. This will show us if the variables and model are statistically significant, and show how much of the variation in Equipment Damage is explained by the two variables. We will be comparing these findings with the interaction model’s results to determine which model is sufficient.

# 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).

3.1.3 Residual Analysis

Next, we will be looking at the residual analysis for this model to determine if the residual assumptions are met. If the assumptions are not met and we decide to move forward with this model, then we will use a Box–Cox transformation and bootstrap confidence intervals for robust inference.

# 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 show 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 residual assumptions are not met.

3.2 Interaction Term

Now, we will take a look at the interaction model. If this model has a higher R-squared value and is statistically significant, we will move forward with this model instead of the previous one.

3.2.1 Inferential Statistics

Let’s take a look at the inferential statistics so we can find out if this model is better than the main effect model.

# 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) and the interaction (p < 2e-16) are statistically significant. The Track Damage main effect by itself is not (0.073, p = 0.095). Because the interaction is significant, the effect of Cars Derailed on Equipment Damage increases as Track Damage increases. When both predictors are zero, the expected Equipment Damage is about $3,400 (intercept not significant). Each additional derailed car increases expected Equipment Damage by about $17,350 when Track Damage = 0; each additional dollar of Track Damage is associated with about $0.07 of additional Equipment Damage at Cars Derailed = 0. The model explains 53.39% of the variation (F = 997.3, p < 0.001). Because the model includes an interaction, main-effect slopes are conditional (e.g., the Cars Derailed slope is evaluated at Track Damage = 0). Centering Track Damage at its mean would make interpretation more intuitive.

3.2.2 Residual Analysis

Now that we know that we’re going to move forward with this model, let’s see if it meets the residual assumptions. If not, we will be using a Box-Cox transformation to come to our conclusion.

# 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.

4 Remedies and Robust Inference

Since the residual assumptions were not met, we are going to move forward with a Box-Cox transformation, followed by 95% bootstrap confidence intervals.

4.1 Box-Cox Transformation

We are going to start with a Box-Cox transformation.

library(MASS)

k <- -min(myData$EqpDamg, na.rm = TRUE) + 1
myData$EqpDamg_pos <- myData$EqpDamg + k
model3_pos <- lm(EqpDamg_pos ~ CarsDer * TrkDamg, data = myData)

# Profile over a tight grid so λ≈0.20 is obvious
bc <- boxcox(model3_pos, lambda = seq(0.10, 0.35, 0.001), plotit = FALSE)
lam <- bc$x; ll <- bc$y
cutoff <- max(ll) - qchisq(0.95, df = 1) / 2
ci <- range(lam[ll >= cutoff])
lam_hat <- lam[which.max(ll)]

plot(lam, ll, type="l", xlab=expression(lambda), ylab="Log-Likelihood",
     xlim=c(0.10, 0.35))
abline(h = cutoff, lty = 3)
abline(v = ci, lty = 3)
abline(v = 0.20, lty = 2, lwd = 2)
mtext(sprintf("λ̂ = %.2f; 95%% CI = [%.2f, %.2f]", lam_hat, ci[1], ci[2]),
      side = 3, line = -1)

We checked whether a simple power transform would make our interaction model more reliable, since a few huge accidents create strong skew and uneven noise. The Box–Cox check pointed clearly to a gentle transform to about a fifth-root (λ ≈ 0.20) with a very tight 95% confidence interval [.2, .21] that rules out both “no transformation” (λ = 1) and a pure log (λ = 0). In plain terms, this transform slightly compresses the biggest damage values so they don’t dominate, which evens out the errors and makes the model’s uncertainty more trustworthy. We fit the regression on this transformed scale for better statistical behavior, then will convert predicted values back into dollars so results stay interpretable. The substantive story doesn’t change: Cars Derailed and the interaction remain reliably positive, while the Track Damage main effect alone is not clearly different from zero; the transformation primarily stabilized variance and improved diagnostics.

4.2 Transforming Our Response

Now that we know λ ≈ 0.20, we will be transforming our response, refitting the interaction model, and then back-transform predictions so results remain interpretable on the dollar scale.

# Transform response with lambda = 0.20 ---
k <- -min(myData$EqpDamg, na.rm=TRUE) + 1
lam <- 0.20
y_pos <- myData$EqpDamg + k
y_bc  <- (y_pos^lam - 1) / lam                 # Box–Cox transform

# Refit the interaction model on transformed scale ---
model_bc <- lm(y_bc ~ CarsDer * TrkDamg, data = myData)

# Inference ---
summary(model_bc)                              # on transformed scale

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

Residuals:
    Min      1Q  Median      3Q     Max 
-96.569  -6.530   1.981   7.651  66.265 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      2.422e+01  3.092e-01  78.323  < 2e-16 ***
CarsDer          1.505e+00  6.314e-02  23.845  < 2e-16 ***
TrkDamg          1.238e-05  3.422e-06   3.616 0.000305 ***
CarsDer:TrkDamg -1.620e-07  1.765e-07  -0.918 0.358749    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.22 on 2617 degrees of freedom
Multiple R-squared:  0.2971,    Adjusted R-squared:  0.2963 
F-statistic: 368.6 on 3 and 2617 DF,  p-value: < 2.2e-16

# Back-transform fitted & predictions to dollars ---
inv_boxcox <- function(z, lam) if (abs(lam) < 1e-12) exp(z) else (lam*z + 1)^(1/lam)

# Naive back-transform of fitted values (good for medians/typicals)
fit_bc   <- fitted(model_bc)
fit_dols <- inv_boxcox(fit_bc, lam) - k

Using the previously selected λ = 0.20, I first added a small constant to Equipment Damage so all values were positive, then applied the power transformation and fit the interaction model (Cars Derailed, Track Damage, and their product) on that transformed outcome. This step reduces the influence of very large losses and produces a more even residual spread, which makes standard errors and diagnostic checks behave better. After estimating the model, I converted the fitted values back to dollars with the inverse transform and removed the shift so results can be read in their original units. Because I used a direct back-transform, these fitted dollars represent typical outcomes rather than unbiased means.

4.3 Bootstrap Confidence Intervals

We then constructed bootstrap confidence intervals to give robust, data-driven intervals for our interaction model, and to see if it matches our Box-Cox transformed model.

# 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

Using 1,000-sample bootstrap confidence intervals on the untransformed model, the effect of Cars Derailed is clearly positive: we estimate an increase of $12.8k–$21.7k per additional derailed car, holding Track Damage at its reference value. The Track Damage main effect is not statistically distinguishable from zero (CI = [−0.124, 0.322]), while the interaction is positively estimated (CI = [0.03, 0.06]), meaning the impact of Cars Derailed grows as Track Damage increases and vice versa. These conclusions match the Box–Cox (λ≈0.20) refit in terms of signs and which effects are detectable; the transformation mainly stabilized variance and improved diagnostics.

5 Conclusion

Equipment Damage in these railroad accidents is driven primarily by the number of Cars Derailed, and that result holds across methods. In the original, raw model, 1,000-sample bootstrap confidence intervals show a clearly positive Cars Derailed effect (about $12.8k–$21.7k per additional car), while the Cars Derailed × Track Damage interaction is also reliably positive (confidence interval entirely above zero), meaning damages rise faster when derailments coincide with higher Track Damage. The Track Damage main effect alone is not consistently different from zero once the interaction is included (its CI spans zero). To improve reliability, we applied a Box–Cox transformation to the response with λ ≈ 0.20 (95% CI [0.20, 0.21]), which stabilized variance. Fitted values were then back-transformed to dollars for interpretation.

These estimates should be read with a few caveats. The data is from 2010, so it may not reflect current technology and safety practices. Coefficients in the transformed model are on a transformed scale (we therefore interpret via back-transformed predictions), and R² on the transformed outcome is not directly comparable to R² on the dollar scale. Because the model includes an interaction, each main-effect slope is conditional; centering Track Damage at its mean would make that interpretation more intuitive.

Looking ahead, updating the analysis with recent accidents and adding operational and engineering covariates like train speed, consist weight, weather, track class/geometry, would likely improve explanatory power. Practically, the results point to two levers: reduce derailments overall, and prioritize maintenance and inspection where car derailment risk and Track Damage risk coincide, since their combination is associated with disproportionately higher losses.

6 References

safetydata.fra.dot.gov