Statistical Models Exercise Answers

Institute of Technology of Cambodia — I4STM

Author

Student Name 

options(repos = c(CRAN = "https://cloud.r-project.org"))

SECTION 1: DATA VISUALIZATION

Task 1: Galton Heredity

install.packages("UsingR")
Warning: unable to access index for repository https://cloud.r-project.org/src/contrib:
  cannot open URL 'https://cloud.r-project.org/src/contrib/PACKAGES'
Warning: package 'UsingR' is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
Warning: unable to access index for repository https://cloud.r-project.org/bin/macosx/big-sur-x86_64/contrib/4.4:
  cannot open URL 'https://cloud.r-project.org/bin/macosx/big-sur-x86_64/contrib/4.4/PACKAGES'
library(UsingR)
Loading required package: MASS
Loading required package: HistData
Loading required package: Hmisc

Attaching package: 'Hmisc'
The following objects are masked from 'package:base':

    format.pval, units
data(Galton)
head(Galton)
  parent child
1   70.5  61.7
2   68.5  61.7
3   65.5  61.7
4   64.5  61.7
5   64.0  61.7
6   67.5  62.2
# Load required libraries
library(ggplot2)
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:Hmisc':

    src, summarize
The following object is masked from 'package:MASS':

    select
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

(a). Contingency table between children and parents’ height

table_Galton <- table(round(Galton$parent, 1), round(Galton$child, 1))
table_Galton
      
       61.7 62.2 63.2 64.2 65.2 66.2 67.2 68.2 69.2 70.2 71.2 72.2 73.2 73.7
  64      1    0    2    4    1    2    2    1    1    0    0    0    0    0
  64.5    1    1    4    4    1    5    5    0    2    0    0    0    0    0
  65.5    1    0    9    5    7   11   11    7    7    5    2    1    0    0
  66.5    0    3    3    5    2   17   17   14   13    4    0    0    0    0
  67.5    0    3    5   14   15   36   38   28   38   19   11    4    0    0
  68.5    1    0    7   11   16   25   31   34   48   21   18    4    3    0
  69.5    0    0    1   16    4   17   27   20   33   25   20   11    4    5
  70.5    1    0    1    0    1    1    3   12   18   14    7    4    3    3
  71.5    0    0    0    0    1    3    4    3    5   10    4    9    2    2
  72.5    0    0    0    0    0    0    0    1    2    1    2    7    2    4
  73      0    0    0    0    0    0    0    0    0    0    0    1    3    0

(b) Scatter plot with regression line

plot(Galton$parent, Galton$child,
main = "Scatter plot of Children Height vs Parents Average Height",
xlab = "Average Height of Parents (inch)",
ylab = "Height of Children (inch)",
pch = 19, col = "skyblue")
abline(lm(child ~ parent, data = Galton), col = "red", lwd = 2)

Task 2: Munich Rent Index (1999)

install.packages("gamlss.data")
Warning: unable to access index for repository https://cloud.r-project.org/src/contrib:
  cannot open URL 'https://cloud.r-project.org/src/contrib/PACKAGES'
Warning: package 'gamlss.data' is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
Warning: unable to access index for repository https://cloud.r-project.org/bin/macosx/big-sur-x86_64/contrib/4.4:
  cannot open URL 'https://cloud.r-project.org/bin/macosx/big-sur-x86_64/contrib/4.4/PACKAGES'
library(gamlss.data)

Attaching package: 'gamlss.data'
The following object is masked from 'package:UsingR':

    grip
The following object is masked from 'package:datasets':

    sleep
library(dplyr)
data(rent99)
rent99 <- data.frame(rent99)
head(rent99)
      rent   rentsqm area yearc location bath kitchen cheating district
1 109.9487  4.228797   26  1918        2    0       0        0      916
2 243.2820  8.688646   28  1918        2    0       0        1      813
3 261.6410  8.721369   30  1918        1    0       0        1      611
4 106.4103  3.547009   30  1918        2    0       0        0     2025
5 133.3846  4.446154   30  1918        2    0       0        1      561
6 339.0256 11.300851   30  1918        2    0       0        1      541

(a) Histograms and Kernel Density Plots

par(mfrow=c(2,2))
hist(rent99$rent, prob=TRUE, main="Net Rent", xlab="Net rent (€)", col="lightblue")
lines(density(rent99$rent), col="red", lwd=2)

hist(rent99$rentsqm, prob=TRUE, main="Net Rent per sqm", xlab="€", col="lightgreen")
lines(density(rent99$rentsqm), col="red", lwd=2)

hist(rent99$area, prob=TRUE, main="Area (sqm)", xlab="sqm", col="lightyellow")
lines(density(rent99$area), col="red", lwd=2)

hist(rent99$yearc, prob=TRUE, main="Year of Construction", xlab="Year", col="lightgray")
lines(density(rent99$yearc), col="red", lwd=2)

(b) Scatter Plots

par(mfrow=c(2,2))
plot(rent99$area, rent99$rent, main="Net Rent vs Area", xlab="Area (sqm)", ylab="Net Rent (€)", pch=19, col="blue")
plot(rent99$area, rent99$rentsqm, main="Rent per sqm vs Area", xlab="Area (sqm)", ylab="Rent/sqm (€)", pch=19, col="darkgreen")
plot(rent99$yearc, rent99$rent, main="Net Rent vs Year", xlab="Year", ylab="Net Rent (€)", pch=19, col="purple")
plot(rent99$yearc, rent99$rentsqm, main="Rent per sqm vs Year", xlab="Year", ylab="Rent/sqm (€)", pch=19, col="orange")

(c) Average and SD Plot

rent_summary <- rent99 %>%
group_by(area) %>%
summarise(avg_rent = mean(rent), sd_rent = sd(rent),
avg_rentsqm = mean(rentsqm), sd_rentsqm = sd(rentsqm))

plot(rent_summary$area, rent_summary$avg_rent, type="l",
main="Average Net Rent vs Area", xlab="Area", ylab="Average Rent")
lines(rent_summary$area, rent_summary$avg_rent + rent_summary$sd_rent, lty=2)
lines(rent_summary$area, rent_summary$avg_rent - rent_summary$sd_rent, lty=2)

(d) Boxplot and Density by Location

boxplot(rentsqm ~ location, data=rent99,
main="Net Rent per sqm by Location",
xlab="Location", ylab="Rent per sqm (€)", col=c("skyblue","lightgreen","lightpink"))

plot(density(rent99$rentsqm), main="Kernel Density: Net Rent per sqm", col="red", lwd=2)

Task 3: Fuel Consumption

install.packages("alr4")
Warning: unable to access index for repository https://cloud.r-project.org/src/contrib:
  cannot open URL 'https://cloud.r-project.org/src/contrib/PACKAGES'
Warning: package 'alr4' is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
Warning: unable to access index for repository https://cloud.r-project.org/bin/macosx/big-sur-x86_64/contrib/4.4:
  cannot open URL 'https://cloud.r-project.org/bin/macosx/big-sur-x86_64/contrib/4.4/PACKAGES'
library(alr4)
Loading required package: car
Loading required package: carData

Attaching package: 'car'
The following object is masked from 'package:dplyr':

    recode
Loading required package: effects
lattice theme set by effectsTheme()
See ?effectsTheme for details.

Attaching package: 'alr4'
The following objects are masked from 'package:UsingR':

    florida, rat, twins
data(fuel2001)
fuel2001 <- data.frame(fuel2001)

(a) Create new variables

fuel2001$Fuel <- 1000 * fuel2001$FuelC / fuel2001$Pop
fuel2001$Dlic <- 1000 * fuel2001$Drivers / fuel2001$Pop
fuel2001$logMiles <- log(fuel2001$Miles)

(b) Summary and Correlation Plot

summary(fuel2001)
    Drivers             FuelC              Income          Miles       
 Min.   :  328094   Min.   :  148769   Min.   :20993   Min.   :  1534  
 1st Qu.: 1087128   1st Qu.:  737361   1st Qu.:25323   1st Qu.: 36586  
 Median : 2718209   Median : 2048664   Median :27871   Median : 78914  
 Mean   : 3750504   Mean   : 2542786   Mean   :28404   Mean   : 77419  
 3rd Qu.: 4424256   3rd Qu.: 3039932   3rd Qu.:31208   3rd Qu.:112828  
 Max.   :21623793   Max.   :14691753   Max.   :40640   Max.   :300767  
      MPC             Pop                Tax             Fuel      
 Min.   : 6556   Min.   :  381882   Min.   : 7.50   Min.   :317.5  
 1st Qu.: 9391   1st Qu.: 1162624   1st Qu.:18.00   1st Qu.:575.0  
 Median :10458   Median : 3115130   Median :20.00   Median :626.0  
 Mean   :10448   Mean   : 4257046   Mean   :20.15   Mean   :613.1  
 3rd Qu.:11311   3rd Qu.: 4845200   3rd Qu.:23.25   3rd Qu.:666.6  
 Max.   :17495   Max.   :25599275   Max.   :29.00   Max.   :842.8  
      Dlic           logMiles     
 Min.   : 700.2   Min.   : 7.336  
 1st Qu.: 864.1   1st Qu.:10.507  
 Median : 909.1   Median :11.276  
 Mean   : 903.7   Mean   :10.914  
 3rd Qu.: 943.0   3rd Qu.:11.634  
 Max.   :1075.3   Max.   :12.614  
pairs(fuel2001[, c("Fuel", "Dlic", "Income", "logMiles", "Tax")], main="Correlation Plot")

SECTION 2: nassCDS DATA ANALYSIS

install.packages("DAAG")
Warning: unable to access index for repository https://cloud.r-project.org/src/contrib:
  cannot open URL 'https://cloud.r-project.org/src/contrib/PACKAGES'
Warning: package 'DAAG' is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages
Warning: unable to access index for repository https://cloud.r-project.org/bin/macosx/big-sur-x86_64/contrib/4.4:
  cannot open URL 'https://cloud.r-project.org/bin/macosx/big-sur-x86_64/contrib/4.4/PACKAGES'
library(DAAG)

Attaching package: 'DAAG'
The following object is masked from 'package:alr4':

    ais
The following object is masked from 'package:car':

    vif
The following object is masked from 'package:MASS':

    hills
data(nassCDS)

Question 1: Seatbelt Use and Mortality

(1). How many indiv used seatbelts

sum(nassCDS$seatbelt == "belted", na.rm=TRUE)
[1] 18573

(2) table of seatbelt usage vs accident outcome

table_seat <- table(nassCDS$seatbelt, nassCDS$dead)
table_seat
        
         alive  dead
  none    6964   680
  belted 18073   500

(3) Proportion test

test_proportion <- function(tab){
  p1 <- tab[1,2]/sum(tab[1,])
  p2 <- tab[2,2]/sum(tab[2,])
  p <- (tab[1,2]+tab[2,2])/sum(tab)
  se <- sqrt(p*(1-p)*(1/sum(tab[1,])+1/sum(tab[2,])))
  z <- (p1 - p2)/se
  pval <- 2*pnorm(-abs(z))
  return(data.frame(Z_statistic=z, P_value=pval))
}
test_proportion(table_seat)
  Z_statistic       P_value
1    22.02001 1.852181e-107

(4) Visualization

barplot(table_seat, beside=TRUE, col=c("skyblue","salmon"),
main="Seatbelt Usage vs Accident Outcome", ylab="Count")

Question 2: Age of Occupants and Accident Outcome

tapply(nassCDS$ageOFocc, nassCDS$dead, mean, na.rm=TRUE)
   alive     dead 
36.85701 44.61525 
tapply(nassCDS$ageOFocc, nassCDS$dead, sd, na.rm=TRUE)
   alive     dead 
17.65621 21.32238 
boxplot(ageOFocc ~ dead, data=nassCDS, main="Occupant Age by Outcome",
xlab="Outcome", ylab="Age", col=c("orange","lightblue"))
stripchart(ageOFocc ~ dead, data=nassCDS, vertical=TRUE, add=TRUE, pch=16, col="blue")

alive <- nassCDS$ageOFocc[nassCDS$dead=="alive"]
dead  <- nassCDS$ageOFocc[nassCDS$dead=="dead"]
t.test(alive, dead, var.equal=FALSE)

    Welch Two Sample t-test

data:  alive and dead
t = -12.302, df = 1256.4, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -8.995522 -6.520963
sample estimates:
mean of x mean of y 
 36.85701  44.61525

(1) mean and standard deviation

# Load the nassCDS dataset
data(nassCDS, package = "DAAG")

# Check the structure of the dataset
str(nassCDS)
'data.frame':   26217 obs. of  15 variables:
 $ dvcat      : Ord.factor w/ 5 levels "1-9km/h"<"10-24"<..: 3 2 2 3 3 4 5 5 2 2 ...
 $ weight     : num  25.1 25.1 32.4 495.4 25.1 ...
 $ dead       : Factor w/ 2 levels "alive","dead": 1 1 1 1 1 1 1 2 1 1 ...
 $ airbag     : Factor w/ 2 levels "none","airbag": 1 2 1 2 1 1 1 1 1 1 ...
 $ seatbelt   : Factor w/ 2 levels "none","belted": 2 2 1 2 2 2 2 1 2 2 ...
 $ frontal    : num  1 1 1 1 1 1 1 1 0 1 ...
 $ sex        : Factor w/ 2 levels "f","m": 1 1 1 1 1 1 2 2 2 1 ...
 $ ageOFocc   : num  26 72 69 53 32 22 22 32 40 18 ...
 $ yearacc    : num  1997 1997 1997 1997 1997 ...
 $ yearVeh    : num  1990 1995 1988 1995 1988 ...
 $ abcat      : chr  "unavail" "deploy" "unavail" "deploy" ...
 $ occRole    : chr  "driver" "driver" "driver" "driver" ...
 $ deploy     : num  0 1 0 1 0 0 0 0 0 0 ...
 $ injSeverity: num  3 1 4 1 3 3 3 4 1 0 ...
 $ caseid     : chr  "2:3:1" "2:3:2" "2:5:1" "2:10:1" ...
# Look at the variables we need
head(nassCDS[, c("dead", "ageOFocc")])
   dead ageOFocc
1 alive       26
2 alive       72
3 alive       69
4 alive       53
5 alive       32
6 alive       22
# mean and standard deviation of age by accident outcome
age_stats <- nassCDS %>%
  group_by(dead) %>%
  summarise(
    mean_age = mean(ageOFocc, na.rm = TRUE),
    sd_age = sd(ageOFocc, na.rm = TRUE),
    n = n()
  )

print(age_stats)
# A tibble: 2 × 4
  dead  mean_age sd_age     n
  <fct>    <dbl>  <dbl> <int>
1 alive     36.9   17.7 25037
2 dead      44.6   21.3  1180

(2) Visualize the distribution of the occupants’ age by accident outcome with data points

# Create boxplot with data points
ggplot(nassCDS, aes(x = dead, y = ageOFocc, fill = dead)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.3, size = 0.8) +
  labs(title = "Distribution of Occupant Age by Accident Outcome",
       x = "Accident Outcome",
       y = "Age of Occupant") +
  theme_minimal() +
  scale_fill_manual(values = c("alive" = "lightblue", "dead" = "lightcoral"))

The boxplot shows the distribution of ages for both groups (alive and dead) with individual data points overlaid. This visualization helps compare the central tendency, spread, and potential outliers in age distribution between survivors and fatalities.

(3) 95% confidence interval for the mean difference of the age of occupant using t distribution.

# Calculate 95% confidence interval for mean age difference using t-distribution

# Extract ages for each group
age_alive <- nassCDS$ageOFocc[nassCDS$dead == "alive"]
age_dead <- nassCDS$ageOFocc[nassCDS$dead == "dead"]

# Perform t-test to get confidence interval
t_test_result <- t.test(age_dead, age_alive, conf.level = 0.95)

# Extract the confidence interval
ci_lower <- t_test_result$conf.int[1]
ci_upper <- t_test_result$conf.int[2]
mean_diff <- t_test_result$estimate[1] - t_test_result$estimate[2]

cat("Mean age difference (dead - alive):", round(mean_diff, 2), "\n")
Mean age difference (dead - alive): 7.76 
cat("95% Confidence Interval: [", round(ci_lower, 2), ",", round(ci_upper, 2), "]\n")
95% Confidence Interval: [ 6.52 , 9 ]
cat("t-statistic:", round(t_test_result$statistic, 2), "\n")
t-statistic: 12.3 
cat("p-value:", format.pval(t_test_result$p.value), "\n")
p-value: < 2.22e-16

Description:

The 95% confidence interval for the mean difference in age between occupants who died and those who survived is [lower bound, upper bound]. This means we are 95% confident that the true mean age difference in the population falls within this interval.

Question 3: Sex, Age, and Safety Devices

(1)

boxplot(ageOFocc ~ sex, data=nassCDS, main="Occupant Age by Sex", col=c("pink","lightblue"))

(2)

sum(nassCDS$ageOFocc > 50 & nassCDS$dead=="alive", na.rm=TRUE)
[1] 5174

(3)

nassCDS$AgeOFocc_class <- ifelse(nassCDS$ageOFocc > 50, 1, 0)

(4-7)

nassCDS_o50 <- nassCDS %>%
filter(ageOFocc > 50) %>%
select(dead, airbag, weight, injSeverity) %>%
na.omit()

dim(nassCDS_o50)
[1] 5573    4
barplot(table(nassCDS_o50$dead, nassCDS_o50$airbag),
beside=TRUE, col=c("skyblue","tomato"),
main="Airbag vs Outcome (Age > 50)", ylab="Count")

barplot(table(nassCDS_o50$injSeverity, nassCDS_o50$airbag),
beside=TRUE, col=c("lightgreen","pink"),
main="Airbag vs Injury Severity (Age > 50)", ylab="Count")

Question 4: Summary Function for Fatal Accidents

nass_analysis <- function(data){
  d <- subset(data, dead=="dead")
  pct_dead <- 100 * nrow(d) / nrow(data)
  pct_male <- 100 * sum(d$sex=="male", na.rm=TRUE)/nrow(d)
  pct_female <- 100 * sum(d$sex=="female", na.rm=TRUE)/nrow(d)
  mode_sev <- names(sort(table(d$injSeverity), decreasing=TRUE))[1]
  age_min <- min(d$ageOFocc, na.rm=TRUE)
  age_max <- max(d$ageOFocc, na.rm=TRUE)
  table_out <- data.frame(pct_dead, pct_male, pct_female, mode_sev, age_min, age_max)
  hist(d$injSeverity, main="Injury Severity among Deaths", xlab="Injury Severity", col="lightblue")
  return(table_out)
}
nass_analysis(nassCDS)

  pct_dead pct_male pct_female mode_sev age_min age_max
1 4.500896        0          0        4      16      97

(1) which of the accident outcome is “dead”?

fatal_cases <- nassCDS %>% filter(dead == "dead")

# Check the result
cat("Total observations in original dataset:", nrow(nassCDS), "\n")
Total observations in original dataset: 26217 
cat("Fatal cases found:", nrow(fatal_cases), "\n")
Fatal cases found: 1180 
cat("Percentage of fatal cases:", round(nrow(fatal_cases)/nrow(nassCDS)*100, 2), "%\n")
Percentage of fatal cases: 4.5 %
# Preview the fatal cases
head(fatal_cases)
    dvcat weight dead airbag seatbelt frontal sex ageOFocc yearacc yearVeh
8     55+ 27.078 dead   none     none       1   m       32    1997    1987
14  25-39 89.627 dead airbag   belted       0   f       54    1997    1994
83    55+ 27.078 dead   none   belted       1   m       67    1997    1992
84    55+ 27.078 dead   none   belted       1   f       64    1997    1992
265   55+ 13.374 dead   none     none       1   m       23    1997    1986
311 25-39 12.383 dead   none   belted       0   f       77    1997    1992
       abcat occRole deploy injSeverity caseid AgeOFocc_class
8    unavail  driver      0           4 2:13:2              0
14  nodeploy  driver      0           4 2:17:1              1
83   unavail  driver      0           4 2:79:1              1
84   unavail    pass      0           4 2:79:1              1
265  unavail  driver      0           4 4:58:1              0
311  unavail    pass      0           4 4:96:2              1

This task filters the dataset to include only cases where the occupant died

(2) The percentage of deaths

# 4.2 Calculate percentage of deaths out of the overall number of observations
total_observations <- nrow(nassCDS)
fatal_count <- nrow(fatal_cases)
death_percentage <- (fatal_count / total_observations) * 100

cat("Overall death statistics:\n")
Overall death statistics:
cat("Total observations:", total_observations, "\n")
Total observations: 26217 
cat("Number of deaths:", fatal_count, "\n")
Number of deaths: 1180 
cat("Percentage of deaths:", round(death_percentage, 2), "%\n")
Percentage of deaths: 4.5 %

(3) Calculate gender percentages among fatalities

# 4.3 Calculate the percentages of females and males among the occupants who died
gender_breakdown <- fatal_cases %>%
  count(sex) %>%
  mutate(percentage = (n / fatal_count) * 100)

print("Gender distribution among fatalities:")
[1] "Gender distribution among fatalities:"
print(gender_breakdown)
  sex   n percentage
1   f 464   39.32203
2   m 716   60.67797
# Extract specific percentages
female_pct <- gender_breakdown$percentage[gender_breakdown$sex == "f"]
male_pct <- gender_breakdown$percentage[gender_breakdown$sex == "m"]

cat("\nDetailed breakdown:\n")

Detailed breakdown:
cat("Female fatalities:", ifelse(length(female_pct) > 0, round(female_pct, 2), 0), "%\n")
Female fatalities: 39.32 %
cat("Male fatalities:", ifelse(length(male_pct) > 0, round(male_pct, 2), 0), "%\n")
Male fatalities: 60.68 %

(4) Find most frequent injury severity

# 4 Show the most frequent severity of their injuries
severity_analysis <- fatal_cases %>%
  count(injSeverity) %>%
  arrange(desc(n))  # Sort by frequency in descending order

print("Injury severity frequency among fatalities:")
[1] "Injury severity frequency among fatalities:"
print(severity_analysis)
  injSeverity    n
1           4 1084
2           3   87
3           2    7
4           1    1
5           5    1
most_frequent_severity <- severity_analysis$injSeverity[1]
most_frequent_count <- severity_analysis$n[1]

cat("\nMost frequent injury severity:", most_frequent_severity, "\n")

Most frequent injury severity: 4 
cat("Number of cases with this severity:", most_frequent_count, "\n")
Number of cases with this severity: 1084 
cat("This represents", round(most_frequent_count/fatal_count*100, 2), "% of all fatalities\n")
This represents 91.86 % of all fatalities

This identifies the most common injury severity level among those who died. We count occurrences of each injury severity level and sort them to find the most frequent one. This reveals what injury severity is most associated with fatal outcomes.

(5) Calculate age range of fatalities

# 4.5 Calculate the minimum and maximum age of the occupant
min_age <- min(fatal_cases$ageOFocc, na.rm = TRUE)
max_age <- max(fatal_cases$ageOFocc, na.rm = TRUE)

cat("Age range:", max_age - min_age, "years\n")
Age range: 81 years

(6) Create injury severity vs age visualization

# histogram with severity of injuries on x-axis and occupant's age on y-axis
ggplot(fatal_cases, aes(x = factor(injSeverity), y = ageOFocc)) +
  geom_boxplot(aes(fill = factor(injSeverity)), alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5, size = 1, color = "darkred") +
  labs(title = "Relationship Between Injury Severity and Occupant Age in Fatal Accidents",
       subtitle = paste("Analysis of", fatal_count, "fatal cases"),
       x = "Injury Severity Level",
       y = "Occupant Age (years)",
       caption = "Each point represents one fatal case") +
  theme_minimal() +
  theme(legend.position = "none") +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "blue")

This creates a visualization showing how age distributes across different injury severity levels.

(7) Complete SINGLE-for-all function

# 4.7 Complete function that returns both numerical and graphical outputs
analyze_fatal_accidents <- function(nassCDS) {
  # 4.1 - Filter fatal cases
  fatal_cases <- nassCDS %>% filter(dead == "dead")
  
  # 4.2 - Calculate death percentage
  total_obs <- nrow(nassCDS)
  fatal_count <- nrow(fatal_cases)
  death_pct <- (fatal_count / total_obs) * 100
  
  # 4.3 - Gender percentages
  gender_pct <- fatal_cases %>%
    count(sex) %>%
    mutate(percentage = (n / fatal_count) * 100)
  
  # 4.4 - Most frequent injury severity
  most_severe <- fatal_cases %>%
    count(injSeverity) %>%
    arrange(desc(n)) %>%
    slice(1)
  
  # 4.5 - Age range
  min_age <- min(fatal_cases$ageOFocc, na.rm = TRUE)
  max_age <- max(fatal_cases$ageOFocc, na.rm = TRUE)
  
  # Create numerical output table
  numerical_output <- data.frame(
    Metric = c("Death Percentage", "Female Fatalities %", "Male Fatalities %", 
               "Most Frequent Severity", "Min Age", "Max Age"),
    Value = c(round(death_pct, 2),
              round(gender_pct$percentage[gender_pct$sex == "f"], 2),
              round(gender_pct$percentage[gender_pct$sex == "m"], 2),
              most_severe$injSeverity,
              round(min_age, 1),
              round(max_age, 1))
  )
  
  # 4.6 - Create plot
  graphical_output <- ggplot(fatal_cases, aes(x = factor(injSeverity), y = ageOFocc)) +
    geom_boxplot(aes(fill = factor(injSeverity)), alpha = 0.7) +
    labs(title = "Injury Severity vs Age in Fatal Accidents",
         x = "Injury Severity", y = "Occupant Age") +
    theme_minimal() +
    theme(legend.position = "none")
  
  return(list(numerical = numerical_output, plot = graphical_output))
}

# Execute the complete function
results <- analyze_fatal_accidents(nassCDS)
print(results$numerical)
                  Metric Value
1       Death Percentage  4.50
2    Female Fatalities % 39.32
3      Male Fatalities % 60.68
4 Most Frequent Severity  4.00
5                Min Age 16.00
6                Max Age 97.00
print(results$plot)

Question 5: Seatbelt Users Summary

seatbelt_users <- nassCDS[nassCDS$seatbelt=="belted", ]
nrow(seatbelt_users)
[1] 18573
table(seatbelt_users$dead)

alive  dead 
18073   500 
table(seatbelt_users$dead, seatbelt_users$occRole)
       
        driver  pass
  alive  14550  3523
  dead     363   137
seatbelt_users <- seatbelt_users[order(seatbelt_users$injSeverity, seatbelt_users$ageOFocc), ]
head(seatbelt_users[order(-seatbelt_users$weight), ], 25)
        dvcat   weight  dead airbag seatbelt frontal sex ageOFocc yearacc
6888    10-24 57871.59 alive airbag   belted       1   f       16    1998
6889    10-24 57871.59 alive airbag   belted       1   m       16    1998
6844    10-24 57871.59 alive airbag   belted       1   f       59    1998
2630    10-24 47463.09 alive airbag   belted       0   f       24    1997
25422   10-24 31694.04 alive airbag   belted       1   m       33    2002
25423   10-24 31694.04 alive airbag   belted       1   f       47    2002
10596   10-24 29533.13 alive airbag   belted       1   m       19    1999
24332   10-24 29301.14 alive   none   belted       1   m       23    2002
10572   10-24 28281.11 alive   none   belted       0   m       43    1999
20993   25-39 28215.77 alive   none   belted       1   m       35    2001
20994   25-39 28215.77 alive   none   belted       1   f       37    2001
11092   10-24 26789.34 alive airbag   belted       0   m       30    1999
10538   10-24 25988.01 alive   none   belted       0   f       34    1999
2507  1-9km/h 25688.37 alive airbag   belted       0   m       36    1997
19732   10-24 25029.20 alive airbag   belted       1   f       31    2001
19733   10-24 25029.20 alive   none   belted       1   m       37    2001
19704 1-9km/h 23428.61 alive   none   belted       0   f       31    2001
15683 1-9km/h 23131.14 alive airbag   belted       0   m       42    2000
15684 1-9km/h 23131.14 alive airbag   belted       0   f       43    2000
15682   10-24 23131.14 alive airbag   belted       1   m       18    2000
19612   10-24 22232.64 alive   none   belted       1   m       36    2001
10521 1-9km/h 20296.81 alive   none   belted       0   m       16    1999
10522 1-9km/h 20296.81 alive   none   belted       0   m       17    1999
24243   10-24 20081.63 alive airbag   belted       0   f       21    2002
24242   10-24 20081.63 alive airbag   belted       0   m       26    2002
      yearVeh    abcat occRole deploy injSeverity   caseid AgeOFocc_class
6888     1998 nodeploy  driver      0           0 48:139:1              0
6889     1998 nodeploy    pass      0           0 48:139:1              0
6844     1990   deploy  driver      1           2  48:97:1              1
2630     1995 nodeploy  driver      0           0  48:92:2              0
25422    1994 nodeploy  driver      0           0  75:24:2              0
25423    1994 nodeploy    pass      0           0  75:24:2              0
10596    1995   deploy  driver      1           0 43:191:1              0
24332    1984  unavail  driver      0           1  48:96:1              0
10572    1992  unavail  driver      0           0 43:174:1              0
20993    1987  unavail  driver      0           0 76:113:2              0
20994    1987  unavail    pass      0           1 76:113:2              0
11092    1996   deploy  driver      1           3 48:131:1              0
10538    1993  unavail  driver      0           0 43:151:1              0
2507     1993 nodeploy  driver      0           0 45:169:1              0
19732    1994   deploy  driver      1           0 48:117:1              0
19733    1994  unavail    pass      0           0 48:117:1              0
19704    1993  unavail  driver      0           3  48:94:2              0
15683    1998 nodeploy  driver      0           0 48:164:2              0
15684    1998 nodeploy    pass      0           0 48:164:2              0
15682    2000 nodeploy  driver      0           2 48:164:1              0
19612    1991  unavail  driver      0           0   48:8:2              0
10521    1990  unavail  driver      0           0 43:135:2              0
10522    1990  unavail    pass      0           0 43:135:2              0
24243    1996 nodeploy    pass      0           0  48:34:2              0
24242    1996 nodeploy  driver      0           0  48:34:2              0

(1) New df containing only occupants who used seatbel

# Create a new data frame which contains only occupants who used seatbelt
seatbelt_users <- nassCDS %>% filter(seatbelt == "belted")

# Check the result
cat("Created new data frame: seatbelt_users\n")
Created new data frame: seatbelt_users
cat("Dimensions:", dim(seatbelt_users), "(rows × columns)\n")
Dimensions: 18573 16 (rows × columns)
cat("First few rows:\n")
First few rows:
head(seatbelt_users)
  dvcat  weight  dead airbag seatbelt frontal sex ageOFocc yearacc yearVeh
1 25-39  25.069 alive   none   belted       1   f       26    1997    1990
2 10-24  25.069 alive airbag   belted       1   f       72    1997    1995
4 25-39 495.444 alive airbag   belted       1   f       53    1997    1995
5 25-39  25.069 alive   none   belted       1   f       32    1997    1988
6 40-54  25.069 alive   none   belted       1   f       22    1997    1985
7   55+  27.078 alive   none   belted       1   m       22    1997    1984
    abcat occRole deploy injSeverity caseid AgeOFocc_class
1 unavail  driver      0           3  2:3:1              0
2  deploy  driver      1           1  2:3:2              1
4  deploy  driver      1           1 2:10:1              1
5 unavail  driver      0           3 2:11:1              0
6 unavail  driver      0           3 2:11:2              0
7 unavail  driver      0           3 2:13:1              0

(2) How many occupants used seatbelt ?

total_seatbelt_users <- nrow(seatbelt_users)

cat("Occupants who used seatbelt:", total_seatbelt_users, "\n")
Occupants who used seatbelt: 18573

(3) Survival outcomes among seatbelt users

# Among seatbelt users, how many died and how many survived?
survival_breakdown <- seatbelt_users %>%
  count(dead) %>%
  mutate(percentage = (n / total_seatbelt_users) * 100)

# Calculate fatality rate
fatalities <- survival_breakdown$n[survival_breakdown$dead == "dead"]
survivors <- survival_breakdown$n[survival_breakdown$dead == "alive"]

cat(
  "Seatbelt users who survived:", 
  survivors,
  "(", round(survivors/total_seatbelt_users*100, 2), 
  "%)\n"
)
Seatbelt users who survived: 18073 ( 97.31 %)
cat(
  "Seatbelt users who died:", 
  fatalities, 
  "(", round(fatalities/total_seatbelt_users*100, 2), 
  "%)\n"
)
Seatbelt users who died: 500 ( 2.69 %)

(4) Role analysis among seatbelt users

# Analyze drivers vs passengers among different outcomes
role_analysis <- seatbelt_users %>%
  group_by(dead, occRole) %>%
  summarise(count = n(), .groups = 'drop') %>%
  group_by(dead) %>%
  mutate(percentage = (count / sum(count)) * 100)

cat("Role distribution among seatbelt users by outcome:\n")
Role distribution among seatbelt users by outcome:
print(role_analysis)
# A tibble: 4 × 4
# Groups:   dead [2]
  dead  occRole count percentage
  <fct> <chr>   <int>      <dbl>
1 alive driver  14550       80.5
2 alive pass     3523       19.5
3 dead  driver    363       72.6
4 dead  pass      137       27.4
# Specific answers for the question
drivers_died <- seatbelt_users %>%
  filter(dead == "dead" & occRole == "driver") %>%
  nrow()

passengers_survived <- seatbelt_users %>%
  filter(dead == "alive" & occRole == "passenger") %>%
  nrow()

cat("\nSpecific counts requested:\n")

Specific counts requested:
cat("Drivers among seatbelt users who died:", drivers_died, "\n")
Drivers among seatbelt users who died: 363 
cat("Passengers among seatbelt users who survived:", passengers_survived, "\n")
Passengers among seatbelt users who survived: 0 
# Additional insight: compare role distribution
cat("\nAdditional insight - Role distribution:\n")

Additional insight - Role distribution:
role_summary <- seatbelt_users %>% count(occRole)
print(role_summary)
  occRole     n
1  driver 14913
2    pass  3660

(5) Sort by injury severity and age

# Sort the data frame according to injury severity and occupant age
seatbelt_sorted <- seatbelt_users %>%
  arrange(desc(injSeverity), desc(ageOFocc))

cat("Data frame sorted by injury severity (descending) and age (descending)\n")
Data frame sorted by injury severity (descending) and age (descending)
cat("First 10 rows of sorted data:\n")
First 10 rows of sorted data:
print(seatbelt_sorted[1:10, c("injSeverity", "ageOFocc", "dead", "occRole", "weight")])
      injSeverity ageOFocc  dead occRole  weight
8618            6       68 alive  driver  25.089
1354            6       57 alive  driver  31.342
23855           5       90 alive  driver  42.379
20225           5       84 alive    pass 190.091
21345           5       84 alive    pass  28.693
17673           5       83 alive  driver 772.047
17954           5       82 alive  driver 225.811
21925           5       81 alive  driver 134.640
21856           5       80 alive  driver 314.908
21872           5       79 alive  driver  44.780
# Alternative: ascending order for age
seatbelt_sorted_alt <- seatbelt_users %>%
  arrange(desc(injSeverity), ageOFocc)

cat("\nAlternative: Sorted by injury severity (desc) and age (asc)\n")

Alternative: Sorted by injury severity (desc) and age (asc)
cat("First 10 rows:\n")
First 10 rows:
print(seatbelt_sorted_alt[1:10, c("injSeverity", "ageOFocc", "dead", "occRole", "weight")])
      injSeverity ageOFocc  dead occRole   weight
1354            6       57 alive  driver   31.342
8618            6       68 alive  driver   25.089
5130            5       16  dead  driver   49.483
10028           5       16 alive    pass  258.285
2185            5       17 alive  driver  404.961
17397           5       17 alive  driver   65.680
655             5       18 alive    pass   33.288
8982            5       18 alive  driver 1408.990
21482           5       18 alive  driver  367.552
25994           5       18 alive    pass   24.828
# Show range of injury severity
cat("\nInjury severity range in seatbelt users:\n")

Injury severity range in seatbelt users:
cat("Min severity:", min(seatbelt_users$injSeverity, na.rm = TRUE), "\n")
Min severity: 0 
cat("Max severity:", max(seatbelt_users$injSeverity, na.rm = TRUE), "\n")
Max severity: 6 
cat("Mean severity:", round(mean(seatbelt_users$injSeverity, na.rm = TRUE), 2), "\n")
Mean severity: 1.5

This sorts the seatbelt users dataframe first by injury severity (most severe first) and then by age (oldest first within each severity level)

(6) Top 25 occupants by weight

# Print the 25 occupants with the highest weight
top_25_weight <- seatbelt_users %>%
  arrange(desc(weight)) %>%
  head(25)

cat("25 seatbelt users with highest weight:\n")
25 seatbelt users with highest weight:
print(top_25_weight[, c("weight", "ageOFocc", "dead", "injSeverity", "occRole", "sex")])
        weight ageOFocc  dead injSeverity occRole sex
6844  57871.59       59 alive           2  driver   f
6888  57871.59       16 alive           0  driver   f
6889  57871.59       16 alive           0    pass   m
2630  47463.09       24 alive           0  driver   f
25422 31694.04       33 alive           0  driver   m
25423 31694.04       47 alive           0    pass   f
10596 29533.13       19 alive           0  driver   m
24332 29301.14       23 alive           1  driver   m
10572 28281.11       43 alive           0  driver   m
20993 28215.77       35 alive           0  driver   m
20994 28215.77       37 alive           1    pass   f
11092 26789.34       30 alive           3  driver   m
10538 25988.01       34 alive           0  driver   f
2507  25688.37       36 alive           0  driver   m
19732 25029.20       31 alive           0  driver   f
19733 25029.20       37 alive           0    pass   m
19704 23428.61       31 alive           3  driver   f
15682 23131.14       18 alive           2  driver   m
15683 23131.14       42 alive           0  driver   m
15684 23131.14       43 alive           0    pass   f
19612 22232.64       36 alive           0  driver   m
10521 20296.81       16 alive           0  driver   m
10522 20296.81       17 alive           0    pass   m
24242 20081.63       26 alive           0  driver   m
24243 20081.63       21 alive           0    pass   f

Question 6: Age, Seatbelt, and Airbag Relationship

par(mfrow=c(1,2))
with(nassCDS, {
boxplot(ageOFocc ~ injSeverity + seatbelt,
main="Age by Injury Severity (Seatbelt)", col=c("skyblue","lightgreen"))
boxplot(ageOFocc ~ injSeverity + airbag,
main="Age by Injury Severity (Airbag)", col=c("pink","lightyellow"))
})