STA321 Assignment 2

1 Introduction

For Assignment 2, we are tasked with finding a data set and selecting one explanatory variable to correlate with the response variable. We will be using the least squares method and 95% bootstrap confidence intervals if applicable.

2 Description of my 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. This data set has a handful of response variables, so I’m going to see which explanatory variables correlate best with number of people killed, number of people injured, number of locomotives derailed, and number of cars derailed. My guess is that speed, equipment damage, and number equipment damage will correlate the strongest with these response variables.Before even digging into the data, it looks like the set does have enough information to answer these questions.

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 Pairwise Scatterplot

library(psych)

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

# Optional: check the structure
str(numData)

tibble [2,621 × 10] (S3: tbl_df/tbl/data.frame)
 $ Acident : num [1:2621] 1 2 3 4 5 6 7 8 9 10 ...
 $ Month   : num [1:2621] 1 1 1 1 1 1 1 1 1 1 ...
 $ Day     : num [1:2621] 1 2 2 2 2 2 2 3 3 3 ...
 $ EqpDamg : num [1:2621] 8485 1500000 103833 15000 9964 ...
 $ TrkDamg : num [1:2621] 333700 0 17615 50000 750 ...
 $ Killed  : num [1:2621] 0 0 0 0 0 0 0 0 0 0 ...
 $ Injured : num [1:2621] 0 0 0 0 0 0 0 0 0 0 ...
 $ Speed   : num [1:2621] 24 0 10 4 45 2 4 9 3 10 ...
 $ LocosDer: num [1:2621] 0 0 0 5 0 0 0 0 0 0 ...
 $ CarsDer : num [1:2621] 13 0 6 0 0 1 3 6 0 0 ...

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

3.1 Finding a Strong Correlation

I wanted to see which variables had the strongest correlation to make the next step easier. As you can see, I had to pivot from my original plan of finding a variable that correlates with number of people killed and number of people injured. This is because that both of these numbers were both very little throughout all of the accidents, with the results typically being 0 or 1. I had to change my regression model to find two other variables that had a correlation. My new explanatory variable is cars derailed and my new response variable is equipment damage. The correlation coefficient for these two variables is 0.66.

# Get only numeric columns
numeric_data <- myData[sapply(myData, is.numeric)]

# Correlation matrix
cor_matrix <- cor(numeric_data, use = "complete.obs")

cor_matrix

              Acident        Month          Day      EqpDamg      TrkDamg
Acident   1.000000000  0.996205159  0.026732926 -0.007057778 -0.003522700
Month     0.996205159  1.000000000 -0.056019706 -0.005761464 -0.004673623
Day       0.026732926 -0.056019706  1.000000000 -0.011989520  0.015826863
EqpDamg  -0.007057778 -0.005761464 -0.011989520  1.000000000  0.476854734
TrkDamg  -0.003522700 -0.004673623  0.015826863  0.476854734  1.000000000
Killed   -0.023449754 -0.023510018  0.004927620  0.023794277  0.004064826
Injured  -0.008792426 -0.006753119 -0.026061616  0.117567797 -0.010711459
Speed    -0.017698831 -0.016700409 -0.002005197  0.239990735  0.162821086
LocosDer -0.022181280 -0.022486719  0.011468588  0.085812854  0.046443288
CarsDer  -0.006447788 -0.005438397 -0.011992079  0.658219030  0.453507160
               Killed      Injured        Speed     LocosDer      CarsDer
Acident  -0.023449754 -0.008792426 -0.017698831 -0.022181280 -0.006447788
Month    -0.023510018 -0.006753119 -0.016700409 -0.022486719 -0.005438397
Day       0.004927620 -0.026061616 -0.002005197  0.011468588 -0.011992079
EqpDamg   0.023794277  0.117567797  0.239990735  0.085812854  0.658219030
TrkDamg   0.004064826 -0.010711459  0.162821086  0.046443288  0.453507160
Killed    1.000000000  0.055314747  0.247128656  0.016628834 -0.048867444
Injured   0.055314747  1.000000000  0.214750720  0.072628280 -0.024335555
Speed     0.247128656  0.214750720  1.000000000  0.007642235  0.126906896
LocosDer  0.016628834  0.072628280  0.007642235  1.000000000 -0.018006320
CarsDer  -0.048867444 -0.024335555  0.126906896 -0.018006320  1.000000000

4 Results and Conclusions

4.1 Least Squares Regression Model

There is a positive linear association between cars derailed and equipment damage.

Cars_derailed <- myData$CarsDer
Equipment_Damage <- myData$EqpDamg
plot(Cars_derailed, Equipment_Damage, pch = 21, col ="red",
     main = "Relationship between Cars Derailed and Equipment Damage")

4.2 Plotting the Residuals

The Residuals vs Fitted plot shows there is not a constant variance in the residuals. The QQ-Plot violates the normality assumption as well. The Scale Location plot shows the same assumption violation that the Residuals vs Fitted plot showed. There are three outliers: accident observation numbers 77, 273, and 981.

parametric.model <- lm(Cars_derailed ~ Equipment_Damage)
par(mfrow = c(2,2))
plot(parametric.model)

par(mar = c(4, 4, 2, 1))
plot(Cars_derailed, Equipment_Damage)

4.3 Bootstrap Regression

Since the normality assumptions were not met, we will be using bootstrap confidence intervals to determine the slope. We will be using a 95% confidence interval.

boot.beta <- boot(data = myData, statistic = boot_fn, R = 1000)


boot.beta

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = myData, statistic = boot_fn, R = 1000)


Bootstrap Statistics :
     original  bias    std. error
t1* -16757.54       0           0
t2*  30025.80       0           0

vec.id <- 1:length(Equipment_Damage)
boot.id <- sample(vec.id, length(Equipment_Damage), replace = TRUE)
boot.Equipment_Damage <- Equipment_Damage[boot.id]
boot.Cars_derailed <- Cars_derailed[boot.id]

B <- 1000
boot.beta0 <- NULL
boot.beta1 <- NULL
vec.id <- 1:length(Equipment_Damage)
for(i in 1:B){
  boot.id <- sample(vec.id, length(Equipment_Damage), replace = TRUE)
boot.Equipment_Damage <- Equipment_Damage[boot.id]
boot.Cars_derailed <- Cars_derailed[boot.id]
boot.reg <- lm(Equipment_Damage[boot.id] ~ Cars_derailed[boot.id])
boot.beta0[i] <- coef(boot.reg)[1]
boot.beta1[i] <- coef(boot.reg)[2]}
boot.beta0.ci <- quantile(boot.beta0, c(0.025, 0.975), type = 2)
boot.beta1.ci <- quantile(boot.beta1, c(0.025, 0.975), type = 2)
boot.coef <- data.frame(rbind(boot.beta0.ci, boot.beta1.ci))
names(boot.coef) <- c("2.5%", "97.5%")
pander(boot.coef, caption="Bootstrap Confidence Intervals of Regression Coefficients")

Bootstrap Confidence Intervals of Regression Coefficients Using bootstrap re sampling with 1,000 replications, we obtained the following 95% confidence intervals for the regression coefficients:

	2.5%	97.5%
boot.beta0.ci	-27647	-5178
boot.beta1.ci	25254	34519

Intercept (boot.beta0.ci): The 95% bootstrap confidence interval ranges from –27,327 to –6,691. This suggests that when the number of cars derailed is zero, the expected equipment damage is significantly below zero. The negative intercept indicates that the model’s prediction at very low values of cars derailed may not be meaningful outside the observed data range.

Slope (boot.beta1.ci): The 95% bootstrap confidence interval ranges from 25,819 to 34,771. This indicates a strong positive relationship: for each additional car derailed, equipment damage is estimated to increase by between $25,800 and 34,800. Because the interval does not include zero, the effect of cars derailed on equipment damage is statistically significant.

4.4 P-Values

We ran inferential statistics to find the estimate and t-value of our slope and intercept.

library(pander)

parametric.model <- lm(Equipment_Damage ~ Cars_derailed, data = myData)
reg.table <- coef(summary(parametric.model))

response_var <- all.vars(formula(parametric.model))[1]
predictors <- all.vars(formula(parametric.model))[-1]

caption_text <- paste0(
  "Inferential statistics for predicting ", response_var,
  " (y) from ", paste(predictors, collapse = ", "), " (x)"
)

pander(reg.table, caption = caption_text)

Inferential statistics for predicting Equipment_Damage (y) from Cars_derailed (x)

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	-16758	4143	-4.045	5.386e-05
Cars_derailed	30026	671	44.74	0

Our slope and intercept are highly significant. For each additional car derailed, equipment damage increases by $30,026.

6 Conclusion

We should be using the bootstrap method because the residual assumptions are not met. We can see that there is not constant variance within the residuals. It also appears that none of the variables had a strong correlation, with the only two variables with even a moderate correlation being Equipment Damage and Cars Derailed.