Can Track Damage and Cars Derailed Predict Equipment Damage in Train Accidents?

1 Introduction

In this project, I will be using both bootstrap sampling methods, bootstrap cases and bootstrap residuals, to find the confidence intervals for my main effects and their interaction term, which are number of Cars Derailed and Track Damage. The problem I am trying to solve is the Equipment Damage in train accidents. I want to see if Track Damage and number of Cars Derailed can predict the amount of Equipment Damage in these accidents.

2 Materials

The data set I chose was Railroad Accidents in the U.S. in 2010. The data was collected via safetydata.fra.dot.gov.

2.1 Description of Data Set

These are all of the variables in the data set. There is sufficient information for addressing the research question because our sample size is very large and we have all of the variables we need. Our “created” variable will be CarsDer:TrkDamg, which is an interaction term between Cars Derailed and Track Damage.

• 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

2.2 Potential Data Challenges

This data set is from 2010, so these statistics may be outdated. We could face the potential problem of our findings not translating to that of train accidents in 2025.

3 Methodology and Analysis

To begin, I’ll make a pairwise scatter plot of all the numerical variables in order to nail down wihc variables to test. I’ll perform a forward variable selection to make sure all of the necessary variables are used. With this forward selection, I’ll be finding the inferential statistics of my regression coefficients. This allows us to know if all of the variables are statistically significant, and provide us with our multiple linear regression equation. Then, I’ll be performing a residual analysis to determine if the assumptions for normality are met. If the assumptions are not met, I’ll make a Box-Cox transformation to fit the model and I will use bootstrap sampling as a remedy for this. We will then look at different goodness-of-fit measures for my model and make a selection on the final model. With the final model, I’ll perform bootstrap records and bootstrap residuals. We will then compare the results of the regular model and the bootstrap results.

4 Results and Conclusions

In this note, I will be performing all of the appropriate tests that I mentioned above. We will begin with a pairwise scatter plot.

4.1 Pairwise Scatter Plot of Numerical Variables

This scatter plot allows us to nail down which variables to focus in on because we see which ones correlate with Equipment Damage.

library(psych)

# Keep only numeric columns
numData <- myData[, sapply(myData, is.numeric)]

# Pairwise scatterplot
pairs.panels(
  numData,
  pch = 21,
  main = "Pair-wise Scatter Plot of numerical variables",
  cex.main = 1.1
)

The only two variables with a decent correlation with Equipment Damage are Cars Derailed (r = .66) and Track Damage (r = .48). These will be the variables that we perform our forward variable selection with.

4.2 Forward Variable Selection

First, I’ll start by finding the inferential statistics of Track Damage because it was the weakest of the two correlations. The goal is to find if it’s statistically significant and if it’s R-squared is good enough to explain our variation in Equipment Damage.

4.2.1 Track Damage

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

# Summary of the model
summary(model1)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-2243114   -36613   -25161    -8330  3164839 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.516e+04  4.454e+03   7.895 4.24e-15 ***
TrkDamg     1.144e+00  4.121e-02  27.763  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 218500 on 2619 degrees of freedom
Multiple R-squared:  0.2274,    Adjusted R-squared:  0.2271 
F-statistic: 770.8 on 1 and 2619 DF,  p-value: < 2.2e-16

Track Damage is statistically significant. (p = 0). Each dollar of Track Damage increases Equipment Damage by $1.14. We will keep this variable in our model, but we will go ahead and test Cars Derailed as well, because Track Damage only explains 22.74% of the variation in Equipment Damage (R² = 0.2274).

4.2.2 Cars Derailed and Track Damage

Now, we’ll take a look at the inferential statistics of the Cars Derailed variable, with our previous variable.

# Multiple Linear Regression
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 
-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

This model looks better than that of 4.1.1. Both variables are statistically significant, and this new model explains 47.29% of the variation in Equipment Damage (R² = .4729). Each additional car derailed (p = 0) increases Equipment Damage by about $25k on average, holding Track Damage constant. Each dollar of Track Damage (p = 0) increases Equipment Damage by about $0.54 on average, holding Cars Derailed constant.

4.2.3 Main Effects with Interaction Term

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 an interaction term would be needed in this model. I am going to perform one last inferential statistics test, this time with an interaction term.

# 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

Our interaction model is definitely the best out of the three. For each additional car derailed (p = 0), Equipment Damage increases by about $17,350 on average, holding Track Damage constant. This is a strong effect.For each dollar of Track Damage (p = .095), Equipment Damage increases by about $0.073 on average, holding Cars Derailed constant. The effect is weak and only borderline significant. The interaction (p = 0) is highly significant. This means the effect of Track Damage on Equipment Damage depends on the number of Cars Derailed (and vice versa). This will be the model we proceed with.

4.3 Residual Analysis

In this section, we’ll be performing a residual analysis with our chosen multiple linear regression (MLR) model.

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

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

The Residuals vs Fitted plot does not have a flat red line. Residuals are larger and more spread out at higher fitted values, suggesting heteroscedasticity (not constant variance). There are extreme outliers at 981 and 1858

The residuals in the Q-Q plot deviate strongly from the straight line, especially in the upper tail. This suggests non-normality of the residuals. There are outliers at 793 and 1858.

In the Scale-Location plot, residual spread increases with fitted values, confirming no constant variance once again. The model fits smaller damages better than larger damages.

In the Residuals vs Leverage plot, 77, 1414, and 2586 have high leverage and high residuals. These cases have influence on the model fit.

In summary, the assumptions are not met. The residuals are non-normal, there is not constant variance, and there are influential outliers. It looks like the issue of heteroscedasticity outweighs the issue of the few outliers, so we will perform a Box-Cox transformation to fit the model.

4.4 Box-Cox Transformation

In this note, I’ll create a Box-Cox transformation with Model 3 (the main effects and interaction term). This will stabilize the variance and improve our residuals.

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, 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

Here is a quick summary of what this plot shows:

• x-axis (λ): the power transformation applied to your response variable Y

• y-axis: log-likelihood of the regression model under each λ

• solid line: shows how well the model fits at each λ.

• dotted vertical line: the estimated “best” λ (maximizes log-likelihood).

• horizontal dashed line (95%): confidence interval range — any λ within that range is considered statistically reasonable.

Our result was λ = 0.222. This suggests that the response variable is highly skewed, and a log transformation is necessary. In this next code chunk, we will be applying the transformation and refitting the model.

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

# Apply Box-Cox transformation with lambda = 0.22
lambda_opt <- 0.22
myData$EqpDamg_trans <- myData$EqpDamg_pos^lambda_opt

# Fit interaction model with transformed response
model3_trans <- lm(EqpDamg_trans ~ CarsDer * TrkDamg, data = myData)

# Model summary
summary(model3_trans)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-25.5956  -1.8140   0.4437   1.9688  18.7263 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      7.047e+00  8.108e-02  86.909  < 2e-16 ***
CarsDer          4.042e-01  1.656e-02  24.414  < 2e-16 ***
TrkDamg          3.508e-06  8.974e-07   3.909 9.52e-05 ***
CarsDer:TrkDamg -2.464e-08  4.628e-08  -0.532    0.594    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.467 on 2617 degrees of freedom
Multiple R-squared:  0.3147,    Adjusted R-squared:  0.3139 
F-statistic: 400.6 on 3 and 2617 DF,  p-value: < 2.2e-16

In our transformed model, our R² is much smaller than the raw model, with 31.29% of the variation in transformed Equipment Damage being explained by the variables. But, our residual standard error, 3.467, is significantly smaller compared to the previous model. All of the coefficients are still significant. Each additional car derailed (p = 0) increases transformed Equipment Damage by about 0.40 units (holding Track Damage constant). Track Damage (p = 0) has a small but statistically significant positive effect on transformed damage. In the transformed model, our interaction term (p = 0.6) is not statistically significant. Despite this term not being significant, we will be keeping it in the model for interpretability purposes. Even though the R² value is much lower, we will stick with this model in order to keep variance stabilized.

4.5 Bootstrap Methods

We will now be using two different bootstrap methods: cases and residuals.

4.5.1 Bootstrap Cases

Let’s run bootstrap cases with the transformed model to find the confidence intervals for our regression coefficients.

## Fit final model (with interaction; remove * if you only want main effects)
mlr_model <- lm(EqpDamg_trans ~ CarsDer * TrkDamg, data = myData)

## Number of bootstrap replicates
B <- 1000       

## Extract dimensions
num.p <- length(coef(mlr_model))  # number of parameters in the model
smpl.n <- nrow(myData)            # sample size

## Zero matrix to store bootstrap coefficients 
coef.mtrx <- matrix(0, nrow = B, ncol = num.p)

## Bootstrap loop
for (i in 1:B){
  bootc.id <- sample(1:smpl.n, smpl.n, replace = TRUE)  # resample cases
  mlr_model_btc <- lm(EqpDamg_trans ~ CarsDer * TrkDamg, 
                      data = myData[bootc.id, ])     
  coef.mtrx[i, ] <- coef(mlr_model_btc)  # store bootstrap coefficients    
}

## Label columns with coefficient names
colnames(coef.mtrx) <- names(coef(mlr_model))

## Save original OLS coefficients
coef_orig <- coef(mlr_model)

## Define histogram function
boot.hist <- function(cmtrx, bt.coef.mtrx, var.id, var.nm){
  x_seq <- seq(min(bt.coef.mtrx[, var.id]),
               max(bt.coef.mtrx[, var.id]), length = 300)

  y_norm <- dnorm(x_seq,
                  mean(bt.coef.mtrx[, var.id]),
                  sd(bt.coef.mtrx[, var.id]))

  highestbar <- max(hist(bt.coef.mtrx[, var.id], plot = FALSE)$density)
  ylimit <- max(c(y_norm, highestbar))

  hist(bt.coef.mtrx[, var.id],
       probability = TRUE,
       main = var.nm,
       xlab = "",
       col = "azure1",
       ylim = c(0, ylimit),
       border = "lightseagreen")

  lines(x_seq, y_norm, col = "red3", lwd = 2)
  lines(density(bt.coef.mtrx[, var.id], adjust = 2), col = "blue", lwd = 2)

  abline(v = cmtrx[var.id], col = "darkgreen", lwd = 2, lty = 2)

  legend("topright",
         legend = c("Normal approx", "Kernel density", "OLS coefficient"),
         col = c("red3", "blue", "darkgreen"),
         lty = c(1,1,2), lwd = 2, bty = "n")
}

## Plot histograms for each coefficient
par(mfrow = c(2,2))  # 2x2 grid for 4 coefficients
boot.hist(coef_orig, coef.mtrx, 1, "Intercept")
boot.hist(coef_orig, coef.mtrx, 2, "CarsDer")
boot.hist(coef_orig, coef.mtrx, 3, "TrkDamg")
boot.hist(coef_orig, coef.mtrx, 4, "CarsDer:TrkDamg")

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

## Bootstrap confidence intervals
num.p <- ncol(coef.mtrx)
btc.ci <- character(num.p)
btc.wd <- numeric(num.p)

for (i in 1:num.p){
  lci.025 <- quantile(coef.mtrx[, i], 0.025, type = 2)
  uci.975 <- quantile(coef.mtrx[, i], 0.975, type = 2)
  btc.wd[i] <- uci.975 - lci.025
  btc.ci[i] <- paste0("[", round(lci.025, 4), ", ", round(uci.975, 4), "]")
}

btc_results <- data.frame(
  Coefficient = colnames(coef.mtrx),
  Estimate = round(apply(coef.mtrx, 2, mean), 4),
  `95% CI` = btc.ci,
  Width = round(btc.wd, 4)
)

kable(btc_results, caption = "Bootstrap Coefficient Estimates and 95% CIs")

Bootstrap Coefficient Estimates and 95% CIs

	Coefficient	Estimate	X95..CI	Width
(Intercept)	(Intercept)	7.0283	[6.8055, 7.27]	0.4645
CarsDer	CarsDer	0.4097	[0.3234, 0.4879]	0.1645
TrkDamg	TrkDamg	0.0000	[0, 0]	0.0000
CarsDer:TrkDamg	CarsDer:TrkDamg	0.0000	[0, 0]	0.0000

The bootstrap analysis reinforces the earlier regression results. The intercept is stable, with a narrow 95% confidence interval well above zero, suggesting a reliable baseline level of transformed damage. The number of Cars Derailed is the dominant predictor, with an estimate around 0.41 and a tight confidence interval [0.33, 0.49], providing strong evidence that derailments consistently increase damage. In contrast, Track Damage hows an extremely small effect that rounds to zero in the summary, indicating its influence is negligible compared to Cars Derailed. Likewise, the interaction term’s bootstrap confidence interval collapses to 0, confirming it does not contribute meaningfully. Overall, the bootstrap results highlight that Equipment Damage is driven almost entirely by the number of Cars Derailed, with little to no additional explanatory power from Track Damage or the interaction.

4.5.2 Residual Bootstrap Regression

Now, we will perform the bootstrap regression method, Compared to the cases technique, this method is much more robust. It resamples rows of the data set with replacement.

hist(sort(mlr_model$residuals), n = 40,
     xlab = "Residuals",
     col = "lightblue",
     border = "navy",
     main = "Histogram of Residuals")

Here is the histogram of our residuals. As expected, it is not normally distributed.

## Final model (with interaction; remove * if you only want main effects)
mlr_model <- lm(EqpDamg_trans ~ CarsDer * TrkDamg, data = myData)

## Residuals
model.resid <- mlr_model$residuals

## Number of bootstrap replicates
B <- 1000
num.p <- length(coef(mlr_model))   # number of parameters
samp.n <- nrow(myData)             # sample size

## Zero matrix to store bootstrap coefficients
btr.mtrx <- matrix(0, nrow = B, ncol = num.p)
colnames(btr.mtrx) <- names(coef(mlr_model))

## Bootstrap loop
for (i in 1:B){
  # Bootstrap response values = fitted + resampled residuals
  bt.response <- mlr_model$fitted.values +
                 sample(model.resid, samp.n, replace = TRUE)

  # Add bootstrap response to data
  myData$bt.response <- bt.response

  # Refit model with bootstrap response
  btr.model <- lm(bt.response ~ CarsDer * TrkDamg, data = myData)

  # Store bootstrap coefficients
  btr.mtrx[i, ] <- coef(btr.model)
}

## Inspect results
head(btr.mtrx)

     (Intercept)   CarsDer      TrkDamg CarsDer:TrkDamg
[1,]    6.963507 0.4278439 4.468685e-06   -8.602684e-08
[2,]    6.892562 0.4313595 2.316918e-06   -1.209087e-08
[3,]    7.104123 0.4295606 2.657571e-06   -4.243206e-08
[4,]    7.215222 0.3905643 2.615681e-06    1.940970e-08
[5,]    7.035613 0.4097809 4.631840e-06   -8.736683e-08
[6,]    7.026921 0.4023421 3.482149e-06   -7.678797e-10

apply(btr.mtrx, 2, mean)   # bootstrap mean estimates

    (Intercept)         CarsDer         TrkDamg CarsDer:TrkDamg 
   7.046363e+00    4.046779e-01    3.518416e-06   -2.475642e-08 

apply(btr.mtrx, 2, sd)     # bootstrap standard errors

    (Intercept)         CarsDer         TrkDamg CarsDer:TrkDamg 
   8.019326e-02    1.641608e-02    9.145254e-07    4.724174e-08 

boot.hist <- function(bt.coef.mtrx, var.id, var.nm){
  ## bt.coef.mtrx = matrix of bootstrap coefficient estimates
  ## var.id = column index of coefficient (1 = Intercept, 2 = CarsDer, etc.)
  ## var.nm = variable name for plot title

  # sequence for normal curve
  x_seq <- seq(min(bt.coef.mtrx[, var.id]),
               max(bt.coef.mtrx[, var.id]),
               length = 300)

  # normal density based on bootstrap mean and sd
  y_norm <- dnorm(x_seq,
                  mean(bt.coef.mtrx[, var.id]),
                  sd(bt.coef.mtrx[, var.id]))

  # determine y-axis limit
  highestbar <- max(hist(bt.coef.mtrx[, var.id], plot = FALSE)$density)
  ylimit <- max(c(y_norm, highestbar))

  # histogram
  hist(bt.coef.mtrx[, var.id],
       probability = TRUE,
       main = var.nm,
       xlab = "",
       col = "azure1",
       ylim = c(0, ylimit),
       border = "lightseagreen")

  # add normal density curve
  lines(x_seq, y_norm, col = "red3", lwd = 2)

  # add kernel density curve
  lines(density(bt.coef.mtrx[, var.id], adjust = 2),
        col = "blue", lwd = 2)
}
par(mfrow = c(2,2))  # 2x2 grid for 4 coefficients
boot.hist(bt.coef.mtrx = btr.mtrx, var.id = 1, var.nm = "Intercept")
boot.hist(bt.coef.mtrx = btr.mtrx, var.id = 2, var.nm = "CarsDer")
boot.hist(bt.coef.mtrx = btr.mtrx, var.id = 3, var.nm = "TrkDamg")
boot.hist(bt.coef.mtrx = btr.mtrx, var.id = 4, var.nm = "CarsDer:TrkDamg")

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

The bootstrap histograms indicate that the model’s coefficient estimates are stable across resamples. The intercept is centered near 7, providing a consistent baseline for the transformed measure of Equipment Damage. The coefficient for CarsDer is tightly clustered around 0.40, reinforcing that the number of Cars Derailed is the primary and most reliable predictor of damage. In contrast, the coefficient for TrkDamg is very close to zero, and the interaction term (CarsDer:TrkDamg) is also centered at zero, suggesting that these variables do not contribute meaningfully once CarsDer is included. Overall, the bootstrap results confirm that Equipment Damage is driven mainly by the number of derailed cars, while Track Damage and its interaction with CarsDer have little to no explanatory value. Let’s finish off our boostrapping with regression confidence intervals.

## number of parameters (coefficients)
num.p <- ncol(btr.mtrx)

## store CIs and widths
btr.ci <- character(num.p)
btr.wd <- numeric(num.p)

for (i in 1:num.p){
  lci.025 <- quantile(btr.mtrx[, i], 0.025, type = 2)
  uci.975 <- quantile(btr.mtrx[, i], 0.975, type = 2)
  btr.wd[i] <- uci.975 - lci.025
  btr.ci[i] <- paste0("[", round(lci.025, 4), ", ", round(uci.975, 4), "]")
}

## original OLS coefficients
coef_orig <- coef(mlr_model)

## results table
btr_results <- data.frame(
  Coefficient = names(coef_orig),
  Estimate = round(coef_orig, 4),
  `95% CI` = btr.ci,
  Width = round(btr.wd, 4)
)

kable(btr_results, caption = "Regression Coefficients with 95% Residual Bootstrap CI")

Regression Coefficients with 95% Residual Bootstrap CI

	Coefficient	Estimate	X95..CI	Width
(Intercept)	(Intercept)	7.0467	[6.8948, 7.2126]	0.3178
CarsDer	CarsDer	0.4042	[0.3705, 0.4359]	0.0654
TrkDamg	TrkDamg	0.0000	[0, 0]	0.0000
CarsDer:TrkDamg	CarsDer:TrkDamg	0.0000	[0, 0]	0.0000

The residual bootstrap results confirm the main findings of your model. The intercept is estimated at about 7.05, with a fairly narrow 95% confidence interval [6.89, 7.21], showing it is stable and well estimated. The coefficient for CarsDer is also highly precise, averaging around 0.40 with a very tight confidence interval [0.37, 0.44]; this reinforces that the number of Cars Derailed is the strongest and most reliable predictor of Equipment Damage. In contrast, both TrkDamg and the CarsDer:TrkDamg interaction collapse to zero, with their bootstrap confidence intervals also equal to [0, 0], meaning they provide no meaningful explanatory power in this model. Overall, the residual bootstrap analysis underscores that Equipment Damage is overwhelmingly explained by the number of derailed cars, while Track Damage and its interaction contribute essentially nothing once CarsDer is included.

#5 Combining and Comparing

Now that we’ve created a handful of models, let’s put all the inferential statistics into a single table and compare.

## OLS coefficients with p-values
cmtrx <- summary(mlr_model)$coefficients  # Estimate, Std. Error, t value, Pr(>|t|)

## Combine with bootstrap CIs
final_results <- cbind(
  round(cmtrx, 4),                # OLS inferential stats
  `Case Bootstrap 95% CI` = btc.ci,
  `Residual Bootstrap 95% CI` = btr.ci
)

## Output clean table
kable(as.data.frame(final_results),
      caption = "Final Combined Inferential Statistics: p-values and Bootstrap CIs")

Final Combined Inferential Statistics: p-values and Bootstrap CIs

	Estimate	Std. Error	t value	Pr(>|t|)	Case Bootstrap 95% CI	Residual Bootstrap 95% CI
(Intercept)	7.0467	0.0811	86.9086	0	[6.8055, 7.27]	[6.8948, 7.2126]
CarsDer	0.4042	0.0166	24.4144	0	[0.3234, 0.4879]	[0.3705, 0.4359]
TrkDamg	0	0	3.9086	1e-04	[0, 0]	[0, 0]
CarsDer:TrkDamg	0	0	-0.5324	0.5945	[0, 0]	[0, 0]

This combined table shows how the classical regression results and both bootstrap approaches line up. The intercept is highly significant (t = 86.9, p < 0.001), with both the case bootstrap CI [6.81, 7.28] and residual bootstrap CI [6.89, 7.21] confirming a stable and positive baseline for transformed Equipment Damage. The coefficient for Cars Derailed is also very strong and precise (estimate ≈ 0.40, t = 24.4, p < 0.001), and both bootstrap methods give tight confidence intervals that exclude zero, reinforcing it as the key predictor. By contrast, Track Damage has a small but statistically significant t-test result (p = 0.0001), yet both bootstrap methods collapse its confidence interval to [0, 0], suggesting the classical test may be picking up noise rather than a meaningful effect. Finally, the Cars Derailed × Track Damage interaction is clearly not significant (p = 0.59) and both bootstrap CIs confirm it contributes nothing. Overall, all three inferential approaches converge on the conclusion that Cars Derailed is the only variable with a consistent and meaningful effect, while Track Damage and its interaction with derailments are negligible once the model is properly checked with bootstrapping.

6 Conclusion

The main finding of this project is that Equipment Damage in train accidents is overwhelmingly explained by the number of Cars Derailed. Both traditional regression and bootstrap methods confirmed that Cars Derailed is a highly significant predictor with tight confidence intervals, while Track Damage and its interaction with derailments consistently collapsed to zero in the bootstrap analysis. This means that, although Track Damage shows a weak effect in classical tests, it does not meaningfully improve prediction once more robust methods are used.

A drawback of this analysis is that the data are from 2010, meaning the findings may not fully represent the current landscape of train accidents. In addition, the models are limited to the variables available in the dataset and do not account for other potential factors such as maintenance practices, weather, or train load. Residual analysis also showed influential outliers and evidence of heteroscedasticity, which led us to a Box-Cox transformation. While this improved variance stabilization, it came at the cost of a lower R², showing that some complexity in the data remains unexplained. Future work could improve the model by incorporating more recent data, testing additional predictors, and exploring nonlinear or mixed-effects approaches to better capture variability.

In terms of application, this model can be used to help rail safety analysts and policymakers prioritize accident prevention and resource allocation. Since the number of Cars Derailed is the clearest driver of Equipment Damage, mitigation strategies should focus on reducing derailments, such as through improved inspection protocols, better track maintenance, and enhanced braking systems. Although Track Damage itself was not a reliable predictor in this study, it remains an important cost component in practice, so combining predictive models with engineering assessments could give decision-makers a fuller picture. Overall, the analysis shows how statistical modeling can inform data-driven strategies to minimize losses from train accidents.

7 References

safetydata.fra.dot.gov