מבוא

כלים דו גלגליים תופסים תאוצה רבה בימינו ( בדגש על כלים חשמליים), בין היתר עקב צפיפות האוכלוסין העולה בערים גדולות והעלייה בנגישות הכלים בזכות שירותים שיתופיים. אמצעים אלו מהווים פתרון תחבורתי המקל על חיי משתמשיהם ושל יתר משתמשי הדרך. אך לצד היתרונות התחבורתיים ניצבים חסרונות וסיכונים. בנושא המחקר, החלטנו לבדוק גורמים אפשריים לפציעות וגורמים נוספים הנקשרים לתאונות בקרב כלים דו גלגליים. נשאף בזכות מחקר זה להצביע על גורמים אשר באמצעותם נוכל לחזות או לשלוט בסיכונים ובחסרונות הנובעים משימוש בכלים דו גלגליים.
לטובת ביצוע מחקר זה, ניעזר במאגרי מידע אודות תאונות שנאספו בעיר ניו-יורק בין השנים 2012-2024, ובעיר לידס שבבריטניה בשנים 2015 ו- 2016. בכל שאלה יפורט אודות השימוש המותאם במאגרים אלו.

האם המשתנים גיל הרוכב, הימצאות רישיון נהיגה וגורם תורם לתאונה משפיעים על הסיכוי לפציעה בתאונה

נבו בטש - 319065033


במטרה למענה על שאלה זו נרכיב מודל לרגרסיה לוגיסטית בינארית על מנת לבחון האם מתקיים קשר בין המשתנים המסבירים לבין המשתנה ‘פציעה’ אשר הוגדר כ: ‘פציעה’ (1), ללא פציעה (0).

לטובת השאלה ניעזר בשלושת מאגרי התאונות (Crashes, Persons, Vehicles):

  • Crashes – בעזרתו נחלץ את המס”ד של כל תאונה בה היה מעורב כלי דו גלגלי.

  • Persons - ממנו נחלץ את משתנה הגיל, מצב רישיון הנהיגה, ונצליב לרוכבים הנכונים ביתר המאגרים.

  • Vehicles - חילוץ משתנה הגורם התורם וחישוב המשתנה המוסבר על סמך הפצועים בתאונה

שלושת המאגרים מסופקים על ידי עיריית ניו-יורק דרך פורטל NYC OpenData, הנתונים נאספו על ידי משטרת ניו יורק בין השנים 2012-2024.


תיאור המודל הפרמטרי


\[\begin{aligned} &\log(\frac{p}{1-p})=\hat{\beta_0}+\hat{\beta_1}x_1+\hat{\beta_2}x_2+\hat{\beta_3}x_3)\\ &\\ &\hat{y}=\frac{1}{1+exp(-(\hat{\beta_0}+\hat{\beta_1}x_1+\hat{\beta_2}x_2+\hat{\beta_3}x_3))} \end{aligned}\]
הפרמטרים להלן:

  • \(p\) - ההסתברות להתרחשות אירוע (1) - פציעה

  • \(\hat{y}\) - אמידת הקשר להתסברות

  • \(\hat{\beta_0}\) - החותך

  • \(\{(\hat{\beta_1},\hat{\beta_2},...,\hat{\beta_k}) | \hat{\beta_k} \in\mathbb{R},1\le k\le 3\}\) - מקדמי הרגרסיה

  • \(\{(x_1,x_2,...x_k) | x_k \in\mathbb{R},1\le k\le 3\}\) - משתנים מסבירים הכוללים - כמותי, בינארי, נומינלי


שיוך הפרמטרים למשתנים מסבירים:

  • \(\hat{\beta_1}\) - Driver License Status - סטאטוס רישיון נהיגה (קיים/לא קיים)

  • \(\hat{\beta_2}\) - Contributing Factor - גורם תורם לתאונה(אי שמירה על מרחק, סינוור, השפעת אלכוהול וכו’)

  • \(\hat{\beta_3}\) - Person Age - גיל הרוכב


השערות המודל

\[\begin{aligned} &H_0=\hat{\beta_1}=\hat{\beta_2}=\hat{\beta_3}=0\\ &\\ &H_1=else \end{aligned}\]
אנו נשער כי נדחה את השערת האפס וכי מתקיים קשר בין המשתנים המסבירים שנקבעו למודל לבין הסיכוי לפציעה בתאונה.


הנחות המודל העיקריות

  1. משתנה מוסבר בינארי או דיכוטומי.

  2. אי תלות בין המשתנים המסבירים - (מולטיקולניאריות נמוכה עד ללא קיימת).

  3. גודל מדגם בגודל מספק.

לפני שנוכל להתאים את המודל הלוגיסטי עלינו קודם כל לטפל בריקים ובחריגים. בנתונים שלרשותנו עולים ריקים במשתנים Person Age ו- Driver License Status.


חבילות בשימוש המחקר

library(tidyverse) # Includes dplyr, ggplot2
library(vroom) # For quick loadtimes of large datasets
library(plotly) # For interactive plots
library(mice) # Iterational imputation of PERSON_AGE
library(glmtoolbox) # GVIF Multicolinearity Test
library(DescTools) # Odds Ratio and R^2
library(forcats) # releveling CONTRIBUTING_FACTOR_1 after level dropping
library(performance) # HL Test
library(pROC) # ROC Curve and model perfomance
library(ggeffects) # for prediction plots


טיפול בריקים

סטאטוס רישיון נהיגה - Driver Licesnse Status

נרצה לברר האם נוכל לשייך את הערכים הריקים לדפוס מסוים על פי ניתוח נתונים רשמיים. לשם כך נעזרנו במספר מאמרים, דוחות מוניציפליים ופדרליים.

ראשית, נגדיר כי אנשים עם היתר נהיגה (כלומר ללא רישיון אך לוקחים שיעורי נהיגה ) שווי ערך לאנשים ללא רישיון נהיגה. לפי Census 2022 ארה”ב, בניו-יורק נכון לשנת 2022 אוכלוסייה המונה 8,335,897 בני אדם. מתוך כלל האוכלוסייה המתועדת -  20.6% מתחת לגיל 18. מאחר ובארה”ב הגיל המותר להוצאת רישיון נהיגה הוא 16 נניח כי 19.5% מכלל תושבי העיר ניו-יורק מתחת לגיל 16. לפיכך על סמך חישוב כ- 6,710,397 בני אדם בעלי זכות להוצאת רישיון. יש לקחת בחשבון כי, מאז החלה מגמת עזיבה באוכלוסיית העיר עקב מגפת הקורונה.

בנוסף, על פי מאמר מבקר העיר ניו-יורק, בשנת 2017 כ- 57% מתושבי העיר הזכאים החזיקו ברישיון נהיגה. על פי נתוני מדינת ניו-יורק העדכניים,  עולה כי כ – 3,032,672 מתושבי העיר ניו-יורק מחזיקים ברישיון (מעודכן לתחילת 2024). נקבל כי כ- 45% מתוך האוכלוסייה הזכאית מחזיקה ברישיון נהיגה. נציין כי פרק הזמן בו תקף רישיון במדינת ניו-יורק עומד על- 8 שנים ומכך נניח כי הנתונים המחושבים תקפים בין השנים 2017-2024. ניקח את שיעור מחזיקי הרישיון שקיבלנו לחומרא (45% אל מול 57% ) ונניח כי שיעורם עומד על 50%.

df_license<-vroom("G:/Data Science/DSTemp/Driver_License__Permit__and_Non-Driver_Identification_Cards_Issued_20240308.csv")
city_type<-c(City=unique(df_license$City[grep("NEW YORK|QUEENS|BROOKLYN|BRONX|STATEN ISLAND", df_license$City, ignore.case=TRUE)]))
df_newyork <- df_license %>% filter(df_license$City %in% city_type & df_license$`Year of Birth`>=1949 & `License Class`!="ID ONLY")

df_valid<- df_newyork %>% filter(df_newyork$Status=="VALID" & df_newyork$`Year of Birth`>=1949) # 1949 - as 75 is the avg age to stop driving in the US
eligible<-8335897*(1-0.195) # Based on Census 2022 - New-Yorkers above the age of 16 
total_valid<-nrow(df_valid)
total_invalid<-nrow(df_newyork)-nrow(df_valid)
precent<-as.numeric(count(df_valid)/(eligible))*100
cat("Residents with valid DL: ", nrow(df_valid),"\nPercentage of License holders: ", precent,"%")
## Residents with valid DL:  3032672 
## Percentage of License holders:  45.19363 %

מסידור נתוני מאגרינו עולה כי שיעור מחזיקי הרישיון בקרב רוכבי כלים דו גלגליים המעורבים בתאונות הוא כ- 65% בעוד שיעורם של אלו שאינם מחזיקים ברישיון עומד על כ-17%. אחוז הריקים הינו כ-17%. לכן, על סמך נתוני מבקר העיר ניו-יורק וחישובנו, הוחלט על שיוך הריקים לקבוצת חסרי הרישיון, וזאת בניסיון לעמוד בפרופורציה המחושבת.

את החלטה זו נוכל לסמוך גם בנתוני דוח ‘Cycling in the City’ המפורסם על ידי עיריית ניו יורק, בו ניתן להבחין בעלייה של 94% ברוכבי אופניים. על פי דוח זה, נכון לשנת 2022 כ-2 מיליון תושבי ניו יורק (מעל גיל 18) משתמשים באופניים וכי יותר מ- 900 אלף רוכבים באופן יום-יומי (שיעור הרוכבים עולה לעומת שיעור מחזיקי רישיון יורד).

להלן תרשים עוגה המתאר את הפרופורציה שבין מחזיקי תושבי ניו יורק המחזיקים רישיון תקף לבין אלו שרישיונם אינו תקף/נשלל:

מחזיקי רישיון / אינם מחזיקים ברישיון מקרב האוכלוסיה הזכאית (מעל גיל 16)


גיל הרוכב - Person Age

למשתנה הגיל נדרש לטפל ראשית בערכים חריגים (שליליים, נמוכים וגבוהים באופן חריג). הוחלט על להגביל את גיל הרוכב לטווח הנע בין 5 ל-75, כאשר ערכים החורגים מן הכלל יהפכו לריקים. מאחר ומדובר במשתנה כמותי נוכל לבצע ‘זקיפה’ (Imputation) ולהשלים את הערכים הריקים באמצעות תחזית המחושבת על סמך תצפיות דומות.

נבצע הליך זה באמצעות חבילת MICE (Rubin 1987) – נשתמש בשיטת  ‘pmm’ אשר תשלים את ערכי המשתנה לאורך העמודה בסדרת איטרציות של מודלי חיזוי, בכל איטרציה יושלם ערכו של Person Age על סמך יתר המשתנים במאגר. האיטרציות ימשכו עד להתכנסות (השלמת כל הערכים הריקים).

#impute by multivariate with MICE - machine learning - predict PERSON_AGE based on the other vars - PERSON_AGE ~ INJURY + CONTRIBUTING_FACTOR_1 + NO_INJURIES
age_constraint <- function(ages) # Constraint of age to prevent impossible values
{
  valid_range <- c(5, 75)
  replace(ages, ages < valid_range[1] | ages > valid_range[2], NA) # replaces impossible values with NA's - for imputation
}
res$PERSON_AGE <- age_constraint(res$PERSON_AGE) #applies constraint on column
imputation_method <- mice(res, method = "pmm", m = 10)  #Setting MICE
res <- complete(imputation_method)# Perform imputation


סטטיסטיקה תיאורית

לפני שניגש לבדוק את הנחות המודל, נתאר באמצעות תרשים Lollipopאת פילוח התאונות לפי סטאטוס החזקת רישיון נהיגה לפי כל אחת מקטגוריות הגורמים התורמים לתאונה: לבדיקה

נבחין לדוגמא כי כאשר בוצע מעבר לא תקין בין נתיבים ,בכ- 500 תאונות הרוכב היה בעל רישיון לעומת 250 ללא, ההפרש בין ערכים אלו עומד על 176 תאונות.

בדיקת הנחות המודל

  1. משתנה מוסבר בינארי או דיכוטומי

    המשתנה המוסבר - Injury מחושב תחילה באמצעות סכימת כלל הפצועים המשויכים למס’ד תאונה, במידה ומספר הפציעות לתאונה גדול מ- 0 נגדיר את ערך המשתנה כ- (1) אחרת, יוגדר כ-(0). מדובר במבנה בינארי ולכן הנחה זו מתקיימת.

  2. אי תלות בין המשתנים המסבירים

    ראשית, נבדוק את מדד GVIF (Generalized Variance Inflation Factors) – אשר מתייחס לקבוצות של משתנים מסבירים (במקרה זה Levels במשתנה Contributing Factor) כאחד, ומודד את רמת המולטיקוליניאריות שבמודל. מדד זה מחושב על ידי מדידת השונות המוסברת של משתנה מסביר אחד על ידי משתנה מסביר אחר. על מנת לעמוד בבדיקה זו נבדוק כי- \(GVIF\le5\) , אחרת ניאלץ לקבוע כי קיימת קורלציה חזקה בין המשתנים המסבירים.  

    ############### Setting up variables for model 
    q2_df$DRIVER_LICENSE_STATUS<-ifelse(q2_df$DRIVER_LICENSE_STATUS=="Licensed",1,0) # Converting strings to 0/1
    q2_df$PERSON_AGE<-as.numeric(q2_df$PERSON_AGE) #verifying numeric type
    q2_df$CONTRIBUTING_FACTOR_1 <- as.factor(q2_df$CONTRIBUTING_FACTOR_1)
    #####-----------

    #####---------- Multicolinearity test
    model <- glm(INJURY ~ ., family=binomial,data = q2_df)
    gvif(model)
##                         GVIF df GVIF^(1/(2*df))
## DRIVER_LICENSE_STATUS 1.0474  1          1.0234
## CONTRIBUTING_FACTOR_1 1.0330 33          1.0005
## PERSON_AGE            1.0241  1          1.0120

\(GVIF\) לכל משתנה לא עולים על 5 ולכן נוכל לקבוע כי ההנחה לאי תלות בין המשתנים המסבירים מתקיימת.

שנית, נבחין בתרשים הקורלציה הבא: לבדיקה


  1. גודל מדגם מספק

על מנת לבדוק את גודל המדגם הנדרש לקבלת תוצאות אמינות ניעזר בתוכנת GPower. על מנת לקבל את התוצאה המבוקשת נדרש לחשב מראש מספר פרמטרים:

  • \(\alpha\) - רמת מובהקות - תוגדר לרמה של 5%.

  • \(Power(1-\beta)\) - עוצמת המבחן - תוגדר בתהאם לרמה של 95%.

  • \(Odds Ratio\) - יחס הסיכויים - יחס ההסתברות לקבלת (1) - פציעה אל מול (0) - אי פציעה.

    \[ OddsRatio=\frac{\exp(\beta_0+\beta_1+\beta_2+\beta_3)}{\exp(\beta_0)} \]

  • \(R^2_\mathrm{Other X}\) - שונות המשתנה המסביר העיקרי המבוטאת על ידי יתר משתני המודל. מאחר ובמודל קיים משתנה נומינלי (Contributing Factor), נחשב שונות מדומה.

  • \(X Distribution\) - התפלגות המשתנים המסבירים - מאחר ו Person Age הוא המשתנה הכמותי היחיד, וגם מתפלג בקירוב לנורמלי, נבחר בהתפלגות נורמלית.


חישוב יחס הסיכויים

model_glm<- glm(INJURY~DRIVER_LICENSE_STATUS+CONTRIBUTING_FACTOR_1+PERSON_AGE,data=q2_df,family=binomial(link = "logit")) ## model fit
odds_ratios <- exp(coef(model_glm)) #get coefficients 
ratio_zero<-data.frame(odds_ratios) #sets as data frame

## proportion of variance by other predictors - DescTools
full_model<-model_glm

nagelkerke_other_x_r2 <- PseudoR2(full_model, which = "McFadden")#R^2 Calculation
## zero_event_prob function to calculate expected H0  
zero_event_prob<-function(odds_ratios)
  {
    return<-exp(odds_ratios)/(1+exp(odds_ratios))
  }
combined_OR <- mean(ratio_zero$odds_ratios<10) ## exclude impossible ratios caused by very low frequencies
Prob_Zero_Event <- zero_event_prob(combined_OR) ##  Ho Probabilty 
odds_probZero_r2<-matrix(c(combined_OR,Prob_Zero_Event,nagelkerke_other_x_r2),ncol=3,dimnames=list(c(),c("odds_Ratio","Zero Event","R^2"))) # to present values nicely
print(odds_probZero_r2)
##      odds_Ratio Zero Event        R^2
## [1,]  0.8333333  0.6970593 0.01773448

\[\begin{aligned} &Odds Ratio = 0.8333333\\ &Pr(Y=1|X=1)\rightarrow p_o=\frac{\exp(Odds Ratio)}{\exp(1+Odds Ratio)}=0.6970593\\ &R^2_\mathrm{Other X}=0.01773448 \end{aligned}\]


נציב בתוכנה את הנתונים:


נבחין כי גודל המדגם הנדרש לביצוע רגרסיה לוגיסטית בר”מ של 5% ובעוצמת מבחן של 95% עומד על 1,902 תצפיות. בפועל אנו פועלים עם מדגם בגודל 18,404 תצפיות – מספר תצפיות גדול באופן ניכר מן הנדרש. לכן אנו מקיימים גם את הנחת המודל הזו.


ניתוח המודל

לאחר שהמודל שנבנה עמד בהנחות, נפעיל את הרגרסיה הלוגיסטית, נבחן את מובהקות המשתנים ונקבע האם מתקיים קשר בין המשתנים המסבירים לבין המוסבר הבינארי - פציעה (1) / אי פציעה (0).

model_glm<- glm(INJURY~DRIVER_LICENSE_STATUS+CONTRIBUTING_FACTOR_1+PERSON_AGE,data=q2_df,family=binomial(link = "logit")) ## model fit
##Model Summary
summary(model_glm)
## 
## Call:
## glm(formula = INJURY ~ DRIVER_LICENSE_STATUS + CONTRIBUTING_FACTOR_1 + 
##     PERSON_AGE, family = binomial(link = "logit"), data = q2_df)
## 
## Coefficients:
##                                                             Estimate Std. Error
## (Intercept)                                                 1.138568   0.247324
## DRIVER_LICENSE_STATUS                                      -0.218620   0.039176
## CONTRIBUTING_FACTOR_1Alcohol Involvement                    0.295546   0.322594
## CONTRIBUTING_FACTOR_1Animals Action                        -0.133081   0.528681
## CONTRIBUTING_FACTOR_1Device Use                            -0.258109   1.249924
## CONTRIBUTING_FACTOR_1Driver Distraction                     0.534506   0.247911
## CONTRIBUTING_FACTOR_1Driver Fatigue                        -0.461157   0.690526
## CONTRIBUTING_FACTOR_1Driver Illness                        13.759493 254.563374
## CONTRIBUTING_FACTOR_1Driver Inexperience                    0.741554   0.272085
## CONTRIBUTING_FACTOR_1Driver Lost Consciousness             13.795007 278.643269
## CONTRIBUTING_FACTOR_1Driver Physical Disability           -15.409470 882.743409
## CONTRIBUTING_FACTOR_1Driverless/Runaway Vehicle            13.794031 623.930981
## CONTRIBUTING_FACTOR_1External Distraction                  13.584620 882.743409
## CONTRIBUTING_FACTOR_1Failure to Keep Right                  0.807992   0.595925
## CONTRIBUTING_FACTOR_1Failure to Yield                       0.854216   0.269633
## CONTRIBUTING_FACTOR_1Following Distance                    -0.041621   0.257892
## CONTRIBUTING_FACTOR_1Glare                                 -0.054952   0.899977
## CONTRIBUTING_FACTOR_1Limited Visibility                     1.007745   0.417992
## CONTRIBUTING_FACTOR_1Obstruction/Debris on Road             1.109273   0.375977
## CONTRIBUTING_FACTOR_1Other                                  1.134247   0.272779
## CONTRIBUTING_FACTOR_1Other Vehicular Issue                  0.567933   0.270682
## CONTRIBUTING_FACTOR_1Passing/Lane Usage Improper            0.333707   0.257517
## CONTRIBUTING_FACTOR_1Pavement Defective                     1.339946   0.348212
## CONTRIBUTING_FACTOR_1Pedestrian/Bicyclist Error/Confusion   0.753761   0.259912
## CONTRIBUTING_FACTOR_1Reaction to Uninvolved Vehicle         1.040790   0.286808
## CONTRIBUTING_FACTOR_1Slippery Road Conditions               0.611505   0.312876
## CONTRIBUTING_FACTOR_1Traffic Control Device Defect         13.759700 440.059182
## CONTRIBUTING_FACTOR_1Traffic Control Disregarded            0.715594   0.268327
## CONTRIBUTING_FACTOR_1Turning Improperly                     0.044491   0.327355
## CONTRIBUTING_FACTOR_1Unsafe Backing                        -1.894444   0.852891
## CONTRIBUTING_FACTOR_1Unsafe Lane Changing                   0.059761   0.308120
## CONTRIBUTING_FACTOR_1Unspecified                            0.203077   0.243802
## CONTRIBUTING_FACTOR_1Vehicle Defect (Accelerator)           0.367993   0.367579
## CONTRIBUTING_FACTOR_1Vehicle Defect (Steering)              0.749809   0.510287
## CONTRIBUTING_FACTOR_1Vehicle Defect (Tire)                 13.739716 207.886372
## PERSON_AGE                                                 -0.004029   0.001439
##                                                           z value Pr(>|z|)    
## (Intercept)                                                 4.604 4.15e-06 ***
## DRIVER_LICENSE_STATUS                                      -5.580 2.40e-08 ***
## CONTRIBUTING_FACTOR_1Alcohol Involvement                    0.916 0.359587    
## CONTRIBUTING_FACTOR_1Animals Action                        -0.252 0.801255    
## CONTRIBUTING_FACTOR_1Device Use                            -0.206 0.836401    
## CONTRIBUTING_FACTOR_1Driver Distraction                     2.156 0.031081 *  
## CONTRIBUTING_FACTOR_1Driver Fatigue                        -0.668 0.504239    
## CONTRIBUTING_FACTOR_1Driver Illness                         0.054 0.956894    
## CONTRIBUTING_FACTOR_1Driver Inexperience                    2.725 0.006421 ** 
## CONTRIBUTING_FACTOR_1Driver Lost Consciousness              0.050 0.960515    
## CONTRIBUTING_FACTOR_1Driver Physical Disability            -0.017 0.986073    
## CONTRIBUTING_FACTOR_1Driverless/Runaway Vehicle             0.022 0.982362    
## CONTRIBUTING_FACTOR_1External Distraction                   0.015 0.987722    
## CONTRIBUTING_FACTOR_1Failure to Keep Right                  1.356 0.175143    
## CONTRIBUTING_FACTOR_1Failure to Yield                       3.168 0.001535 ** 
## CONTRIBUTING_FACTOR_1Following Distance                    -0.161 0.871788    
## CONTRIBUTING_FACTOR_1Glare                                 -0.061 0.951312    
## CONTRIBUTING_FACTOR_1Limited Visibility                     2.411 0.015913 *  
## CONTRIBUTING_FACTOR_1Obstruction/Debris on Road             2.950 0.003174 ** 
## CONTRIBUTING_FACTOR_1Other                                  4.158 3.21e-05 ***
## CONTRIBUTING_FACTOR_1Other Vehicular Issue                  2.098 0.035892 *  
## CONTRIBUTING_FACTOR_1Passing/Lane Usage Improper            1.296 0.195022    
## CONTRIBUTING_FACTOR_1Pavement Defective                     3.848 0.000119 ***
## CONTRIBUTING_FACTOR_1Pedestrian/Bicyclist Error/Confusion   2.900 0.003731 ** 
## CONTRIBUTING_FACTOR_1Reaction to Uninvolved Vehicle         3.629 0.000285 ***
## CONTRIBUTING_FACTOR_1Slippery Road Conditions               1.954 0.050647 .  
## CONTRIBUTING_FACTOR_1Traffic Control Device Defect          0.031 0.975056    
## CONTRIBUTING_FACTOR_1Traffic Control Disregarded            2.667 0.007656 ** 
## CONTRIBUTING_FACTOR_1Turning Improperly                     0.136 0.891891    
## CONTRIBUTING_FACTOR_1Unsafe Backing                        -2.221 0.026337 *  
## CONTRIBUTING_FACTOR_1Unsafe Lane Changing                   0.194 0.846213    
## CONTRIBUTING_FACTOR_1Unspecified                            0.833 0.404869    
## CONTRIBUTING_FACTOR_1Vehicle Defect (Accelerator)           1.001 0.316767    
## CONTRIBUTING_FACTOR_1Vehicle Defect (Steering)              1.469 0.141728    
## CONTRIBUTING_FACTOR_1Vehicle Defect (Tire)                  0.066 0.947304    
## PERSON_AGE                                                 -2.800 0.005105 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 19728  on 18403  degrees of freedom
## Residual deviance: 19378  on 18368  degrees of freedom
## AIC: 19450
## 
## Number of Fisher Scoring iterations: 13

נבחין מתוך ניתוח המודל במספר נקודות עיקריות:

  • המשתנה המסביר של סטאטוס רישיון הנהיגה (מופיע בתור DRIVER_LICENSE_STATUS) מובהק במידה - חזקה, כאשר \(pValue=2.40e^{-08}\approx0<\alpha=5\%\)

  • המשתנה המסביר של גיל הרוכב (מופיע בתור PERSON_AGE) מובהק במידה מתונה, כאשר -\(pValue=0.005105=0.51\%<\alpha=5\%\)

  • חלק מרמות המשתנה של ה’גורם התורם’  (מופיעים כחלק מ CONTRIBUTING_FACTOR_1_X) מובהקים במידה חזקה, בינונית וחלשה. ניקח כדוגמא את –‘פגם במדרכה’ (Pavement Defective) אשר יצא מובהק במידה חזקה, כאשר - \(pValue=0.000119=0.02\%<\alpha=5\%\)

  • מנגד, חלק מרמות המשתנה של ה’גורם התורם’ אינם מובהקים, לדוגמא – ‘החלפת נתיב באופן לא בטוח’, אשר לא התקבל כמובהק מאחר ו- \(pValue=0.846213=84.62\%\nless\alpha=5\%\)

לפיכך, לפני שנבנה מודל חדש בו יכללו רק המשתנים המובהקים, נעדכן את משתנה הגורם התורם (Contributing Factor) כך שיכלול רק את הרמות המובהקות.

#filter out non significant levels observations
q2_df <- q2_df[q2_df$CONTRIBUTING_FACTOR_1 %in% significant_levels, ]
q2_df$CONTRIBUTING_FACTOR_1 <- droplevels(q2_df$CONTRIBUTING_FACTOR_1) # dropping unused levels
# Relevel the factor levels  - forcats
q2_df$CONTRIBUTING_FACTOR_1 <- fct_relevel(
  q2_df$CONTRIBUTING_FACTOR_1,
  significant_levels
)


לאחר התאמת המודל מחדש, נוכל לבחון את חיזוי המודל לסיכוי בפציעה בהתאם לגורם התורם לתאונה:



בתרשים העמודות המתאר את הקשר שבין הרמות השונות של ה’גורם התורם’ לבין הסיכויים החזויים לאירוע (1) פציעה ואירוע (0) אי פציעה, נוכל לראות כי במרבית סוגי הגורמים הסיכוי החזוי לפציעה עומד בטווח 80%-90%. כלומר, הסיכוי לפציעה גבוה במיוחד בכל הרמות מלבד ‘Unsafe Backing’


בדיקת טיב המודל

לבדיקת טיב ההתאמה של המודל שבנינו ניעזר במבחן Hosmer-Lemeshow (HL) המיועד לבחינת טיב ההתאמה רק כאשר משתנה התגובה (במקרה זה ‘פציעה’) בינארי. ניקח בחשבון כי מבחן זה רגיש למספר הקבוצות (bin) שנבחר (מספר הנע בין 10 ל-20):

השערות המבחן:

\[\begin{aligned} &H_0=\textrm{Specified model is correct}\\ &\\ &H_1=\textrm{Model is incomplete} \end{aligned}\]

##Fitting model after levels dropping 
model_glm<- glm(INJURY~DRIVER_LICENSE_STATUS+CONTRIBUTING_FACTOR_1+PERSON_AGE,data=q2_df,family=binomial(link = "logit")) ## model fit
#model fit -with performance package 
performance_hosmer(model_glm, n_bins = 10)
## # Hosmer-Lemeshow Goodness-of-Fit Test
## 
##   Chi-squared: 6.425
##            df: 8    
##       p-value: 0.600
## Summary: model seems to fit well.

תוצאות המבחן - \(pValue=0.6=60\%\nless\alpha=5\%\), מצביעות על כך שלא נדחה את \(H_0\) וכי יש התאמה של המודל. עקב רגישותו של מבחן זה נבדוק גם את עקומת ROC אשר השטח תחתיה יספק לנו מדד מצטבר של ביצועי המודל:

השטח תחת עקומת ה- ROC נקרא גם AUC (Area Under Curve), שטח זה משקף את יכולת המודל להבין בין אירוע ה- (1) פציעה לבין אירוע ה- (0) אי פציעה. ערכי AUC נעים במרבית המקרים בטווח שבין 0.5 לבין 1.0, כש 0.5 מצביע על כך שיכולת ההבחנה של המודל איננה טובה יותר מאשר מניחוש אקראי בין שני סוגי האירועים. ערכי AUC מעל 0.8 מצביעים על יעילות המודל.

לפיכך, נסיק כי קיים קשר בין גיל הרוכב, סטאטוס רישיון הנהיגה והגורם התורם לתאונה לבין הסיכוי לפציעה בתאונה. אך מאידך המודל שהוצע לא יכול לחזות באופן טוב את הסיכוי לפציעה בתאונה ( רק בקצת יותר מ- 50% מהמקרים המודל יחזה נכון). את נסיבות התוצאה נוכל לנסות להסביר ככל הנראה על ידי דרכי הטיפול בריקים, ואמצעי המיגון המועטים של רוכבי כלים דו גלגליים שעלולים להביא לפציעה במקרים רבים מתאונות הדרכים.

האם יש קשר בין גיל ומין הנהג לבין כמות כלי הרכב הניזוקים בתאונה?

עדי בן דור - 208061481


בכדי לענות על שאלה זו נבנה מודל של רגרסיה פואסונית ונבדוק האם יש קשר בין המשתנים המסבירים: גיל הרוכב – משתנה כמותי בדיד, מין הרוכב – משתנה נומינלי למשתנה המוסבר: כמות כלי רכב ניזוקים בתאונה – משתנה ספירה בדיד.

רגרסיה פואסונית משמשת לנתוני ספירה כאשר המשתנה המוסבר מייצג ספירות של אירועים המתרחשים במרווח זמן או מרחב קבוע. בהתחשב בכך שהמשתנה המוסבר שלנו הוא ספירה במרחב (מספר כלי הרכב שניזוקו בתאונה), רגרסיה פואסונית יכולה להיות בחירה מתאימה לניתוח נתונים אלו, עם זאת חשוב לוודא שהנחות המודל אכן מתקיימות.

לשם כך נשתמש ב2 בסיסי נתונים של העיר ‘ניו-יורק’:

  1. Vehicles – מתעסק בנתוני כלי רכב, ממנו נשלוף את המספר המזהה של התאונה ונחשב את כמות כלי הרכב הניזוקים על סמך נזק לכל רכב.
  2. Person – מתעסק בנתוני אנשים-נפגעים, ממנו נשלוף את משתנה הגיל, מין הנהג ונצליב לרוכבים מטבלת כלי הרכב.

הנתונים נאספו דרך משטרת ניו-יורק בין השנים 2012-2024.

תיאור פרמטרי של המודל


\(\log{λ_i}=\hat{β}_0+\hat{β}_1x_1+\hat{β}_2x_2\)


שיוך הפרמטרים למשתנים המסבירים

\(\hat{β}_1\) –age גיל הנהג

\(\hat{β}_2\) –sex מין הנהג


השערות המודל

\[\begin{aligned}&H_0:\hat{β}_1=\hat{β}_2=0\\ &\\ &H_1:else\end{aligned}\]


אנו נשער כי נדחה את השערת האפס וכי קיים קשר בין המשתנים המסבירים לבין המשתנה המסביר, כלומר שיש קשר בין גיל ומין הנהג לבין כמות כלי הרכב הניזוקים בתאונה.

הנחות המודל עבור רגרסיה פואסונית

  1. משתנה מוסבר פואסוני - המשתנה המוסבר הוא ספירה ליחידת זמן או מרחב המתואר על ידי התפלגות פואסון.
  2. אי תלות - התצפיות חייבות להיות בלתי תלויות זו בזו.
  3. תוחלת = שונות - הממוצע של משתנה מקרי פואסון חייב להיות שווה לשונות שלו.
  4. לינאריות – בדיקה שהקשר בין הספירות הנצפות והערכים החזויים הוא ליניארי בקירוב.

לפני שנתאים את המודל הפואסוני נטפל קודם בערכים ריקים ובחריגים. בנתונים שלרשותנו עולים ריקים במשתנים Person Age ו- Person Sex.

החבילות בהן נשתמש

library(vroom)  # For quick loadtimes of large datasets
library(tidyverse)
library(mice)# Iterational imputation of PERSON_AGE
library(plotly) #for pie
library(ggplot2)
library(plotrix) ##graph
library(ggcorrplot) #for cor plot
library(hrbrthemes) # for scatterplot


טיפול בריקים - Person Sex

#Turns NA data to "Unknown" and change "F" and "M"
df_p$PERSON_SEX<- ifelse(df_p$PERSON_SEX == ""|is.na(df_p$PERSON_SEX)|df_p$PERSON_SEX == "U", "Unknown", df_p$PERSON_SEX)
df_p$PERSON_SEX<- ifelse(df_p$PERSON_SEX == "F", "Female", df_p$PERSON_SEX)
df_p$PERSON_SEX<- ifelse(df_p$PERSON_SEX == "M", "Male", df_p$PERSON_SEX)
# View the summary statistics
print(summary_by_gender)
##    gender count
## 1  Female 10408
## 2    Male 79182
## 3 Unknown 21385

טיפול בריקים וחריגים - Person Age

תחילה נטפל בערכים החריגים למשתנה הגיל (שליליים, נמוכים וגבוהים באופן חריג). הוחלט על להגביל את גיל הרוכב לטווח הנע בין 5 ל-75, כאשר ערכים החורגים מן הכלל יהפכו לריקים. מאחר ומדובר במשתנה כמותי נוכל לבצע ‘זקיפה’ (Imputation) ולהשלים את הערכים הריקים באמצעות תחזית המחושבת על סמך תצפיות דומות.

נבצע הליך זה באמצעות חבילת MICE (Rubin 1987) – נשתמש בשיטת ‘pmm’ אשר תשלים את ערכי המשתנה לאורך העמודה בסדרת איטרציות של מודלי חיזוי, בכל איטרציה יושלם ערכו של Person Age על סמך יתר המשתנים במאגר. האיטרציות ימשכו עד להתכנסות (השלמת כל הערכים הריקים).

#Treatment of minus values at age
merged_inner$Age<- ifelse((merged_inner$Age)<=4 ,NA,merged_inner$Age)
#Treatment of big values at age
merged_inner$Age<- ifelse((merged_inner$Age)>75 ,NA,merged_inner$Age)
#Treatment of NA values at age
df_withoutk<-data.frame(Age=merged_inner$Age,Sex=merged_inner$Sex,Vehicles_damaged_per_accident=merged_inner$Vehicles_damaged_per_accident)
imputed_age <- mice(df_withoutk,method="pmm", m = 1)
## 
##  iter imp variable
##   1   1  Age
##   2   1  Age
##   3   1  Age
##   4   1  Age
##   5   1  Age
merged_inner <- complete(imputed_age)

summary(age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00   24.00   31.00   33.72   42.00   75.00


להלן תרשים עוגה המתאר את התפלגות הנהגים, הנהגות, ואלו שלא הזדהו:


בדיקת הנחות המודל

  1. משתנה מוסבר פואסוני - המשתנה המוסבר Vehicles damaged per accident, כמות כלי רכב ניזוקים, מחושב באמצעות סכימת רכבים ניזוקים והוא נע בין 0-4, משתנה כמותי בדיד הנספר ליחידת תאונה אחת.
  2. אי תלות בין המשתנים המסבירים – קורלציה נמדדת בדרך כלל באמצעות מקדם מתאם פירסון, המעריך את החוזק והכיוון של הקשר הליניארי בין שני משתנים רציפים. עם זאת, במקרה שלנו קיים משתנה בדיד אחד (גיל) ומשתנה קטגורי אחד (מין), לכן לא נוכל להשתמש במקדם פירסון. נמיר את המודל למטריצה שתגרום לפירוק מגדר הנהג ל3 משתני דמה (Female /NOT Female , Male /NOT Male , Unknown /NOT Unknown).


כעת נוכל לטעון כי קיים אי תלות בין המשתנה הכמותי גיל לבין המשתנה האיכותי מין.


  1. תוחלת = שונות -
# Calculate mean
mean_value <- mean(damaged)
print(mean_value)
## [1] 1.021212
# Calculate variance
variance_value <- var(damaged)
print(variance_value)
## [1] 1.273578


קיבלנו שהתוחלת שווה בקירוב לשונות . הם אמנם לא בדיוק שווים, אבל הם קרובים יחסית. בהתפלגות פואסון, הממוצע והשונות צריכים להיות שווים תיאורטית. הפער הקל עשוי לנבוע מגודל המדגם הסופי או השונות בנתונים.

  1. לינאריות

ניתוח המודל

לאחר שבדקנו את הנחות המודל וראינו כי הן מתקיימות, נפעיל רגרסיה פואסונית ונסכם את התוצאות לפיהן נקבע האם ישנו קשר בין המשתנים המסבירים שנבחרו למודל לבין המשתנה המוסבר אשר מתאר משתנה בדיד של כמות כלי רכב ניזוקים בתאונה.


summary(model)
## 
## Call:
## glm(formula = damaged ~ age + sex, family = poisson, data = merged_inner)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.0173043  0.0121644  -1.423    0.155    
## age          0.0025811  0.0002257  11.434   <2e-16 ***
## sexMale      0.1508284  0.0099882  15.101   <2e-16 ***
## sexUnknown  -1.9288507  0.0197537 -97.645   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 141344  on 110974  degrees of freedom
## Residual deviance: 112954  on 110971  degrees of freedom
## AIC: 274532
## 
## Number of Fisher Scoring iterations: 6


כפי שניתן לראות מן המודל : כל המשתנים המסבירים, מגדר וגיל הנהג התקבלו כמובהקים, לכן נוכל להסיק כי קיים קשר בין המשתנים המסבירים למשתנה המוסבר.

ישנה השפעה של המשתנים מסבירים, גיל ומין הנהג, על המשתנה המוסבר אשר מוגדר ככמות כלי הרכב הניזוקים בתאונה.


ויזואליזציה

תרשים כינור (צפיפות) הממחיש התפלגות של משתנה ‘כמות כלי רכב ניזוקים’ לפי קטגוריות המגדר.

רוחב הכינור מייצג את צפיפות הנתונים בערכים שונים, כאשר חלקים רחבים יותר מציינים תצפיות נוספות.

גרף קופסה המציג את המגדר לפי גיל הרוכב . תרשים קופסא מציג את התפלגות הנתונים על סמך חמישה נתונים סטטיסטיים מסכמים: מינימום, רבעון ראשון (Q1), חציון (רבעון שני או Q2), רבעון שלישי (Q3) ומקסימום.

ה”קופסה” מייצגת את הטווח הבין-רבעוני (IQR), שמתפרש מהרבעון הראשון (Q1) לרבעון השלישי (Q3). החציון מצוין בקו בתוך התיבה.

ה”חריצים” בתרשים מייצגים רווח סמך סביב החציון. החריצים של קופסאות המגדר אינן חופפות, זה מצביע על כך שהחציונים שונים באופן משמעותי.


תרשים פיזור של כמות כלי רכב ניזוקים בתאונה לפי גיל הרוכב.

הנקודות המופיעות בגרף הינן מגדר הרוכב לפי המקרא.

נראה את ניתוח הקשר בין מספר כלי רכב ניזוקים בתאונה לבין גיל, תוך התחשבות במגדר של הרוכבים. ניתן לראות כל מגמה, דפוס או מתאם בין המשתנים הללו. בנוסף,קידוד הצבע מאפשר לנו להשוות חזותית את התפלגות הנקודות בין קטגוריות מגדר שונות. ניתן לזהות למשל כי לאחר גיל 60 נשים פחות היו מעורבות בתאונות בהן נפגעו כלי רכב.

האם קיים קשר בין השימוש בקסדה לבין המצב ההכרתי של הרוכב לאחר התאונה?

יובל ברששת - 318745635

בשאלה זו נחקור את קיומו של הקשר שבין השימוש באמצעי בטיחות לבין המצב הקוגניטיבי של המעורבים בתאונה.


library(dplyr)
library(ggplot2)
library(vroom)
library(tidyverse)
library(summarytools)
library(ggthemes)


נטען את המאגרים Crashes, Person ו- Vehicles בעזרת חבילת Vroom.


crashes_df<- vroom("G:/Data Science/DSTemp/Motor_Vehicle_Collisions_-_Crashes_20240219.csv")
person_df<- vroom("G:/Data Science/DSTemp/Motor_Vehicle_Collisions_-_Person_20240219.csv")
Vehicle_df <- vroom("G:/Data Science/DSTemp/Motor_Vehicle_Collisions_-_Vehicles_20240219.csv")

הגדרת פונקציות

פונקציה לחילוץ סוגי כלים

  • Original - ה- Data Frame עליו תבוצע הפונקציה.

  • Name - מילת המפתח לחיפוש.


getType<-function(original,name)
  {
  rbind(original,data.frame(Name=unique(crashes_df$`VEHICLE TYPE CODE 1`[grep(name,crashes_df$`VEHICLE TYPE CODE 1`, ignore.case=TRUE)])))
  }

פונקציה לטיפול בריקים

נזין שתי עמודות של ‘גורם תורם לתאונה’ ונבדוק האם שתיהן ריקות. במידה והתנאי מתקיים נניח שלא היו כלים מעורבים ונחליף את הערכים הריקים ב – ‘None’.


selectEmpty<-function(col1,col2) 
    {
      selected_rows  <- crashes_df[rowSums(is.na(crashes_df[c(col1,col2)]) | crashes_df[c(col1,col2)] == "") == 2, ]
      crashes_df[rownames(selected_rows),col1]<-"None"
      crashes_df[rownames(selected_rows),col2]<-"None"
      return (crashes_df) 
    }

עיבוד מקדים של הנתונים

נטפל בערכים הריקים בעמודות ‘סוג כלי’ ו- ‘גורם תורם לתאונה’:


crashes_df<-selectEmpty("VEHICLE TYPE CODE 2","CONTRIBUTING FACTOR VEHICLE 2")
    crashes_df<-selectEmpty("VEHICLE TYPE CODE 3","CONTRIBUTING FACTOR VEHICLE 3")
    crashes_df<-selectEmpty("VEHICLE TYPE CODE 4","CONTRIBUTING FACTOR VEHICLE 4")
    crashes_df<-selectEmpty("VEHICLE TYPE CODE 5","CONTRIBUTING FACTOR VEHICLE 5")

נחלץ מתוך הנתונים רק את המילים בהן ‘סוג כלי’ מכיל מילות מפתח כדוגמת – ‘bike’, ‘scoo’ או- ‘cyc’ בעזרת פונקציית grepl. פעולה זו תוודא כי נתמקד רק בסוגי הכלים הרלוונטיים לניתוח.


twoWheeled<-data.frame(Name=unique(crashes_df$`VEHICLE TYPE CODE 1`[grep("bike", crashes_df$`VEHICLE TYPE CODE 1`, ignore.case=TRUE)]))
   twoWheeled<-getType(twoWheeled,"scoo")
   twoWheeled<-getType(twoWheeled,"cyc")

לטובת הנוחות, ניצור Data Frame המכיל את עמודות ’ מס’ סידורי של התאונה’, ‘סוג כלי i’. באמצעות שימוש בפורמט long ניצור רשימה הכוללת את כלל הכלים הדו-גלגליים אשר היו מעורבים בתאונות דרכים.

פעולות שבוצעו

  1. נבחרו כלל העמודות ששמן מכיל ’ VEHICLE TYPE CODE’

  2. יצירת עמודה חדשה המכילה את שמות העמודות המקוריים

  3. יצירת עמודה חדשה המכילה את הערכים

  4. שינוי ערכי העמודות לאותיות קטנות (Lower Case)


crashes_df_selected<-crashes_df[, c("COLLISION_ID", "VEHICLE TYPE CODE 1", "VEHICLE TYPE CODE 2", "VEHICLE TYPE CODE 3", "VEHICLE TYPE CODE 4", "VEHICLE TYPE CODE 5")]

df_long <- crashes_df_selected %>%
  pivot_longer(
    cols = starts_with("VEHICLE TYPE CODE"), 
    names_to = "Vehicle_Column", 
    values_to = "Vehicle_Type_Code" 
  )

יצירת מסגרת נתונים על כלים דו-גלגליים

שלפנו את הנתונים המתקבלים על סמך השמות שחילצנו לתוך twoWheeled, תוך בחירת העמודות הרלוונטיות – מס’ סידורי של התאונה, אמצעי בטיחות, ומצב קוגניטיבי של המעורב


crashes_df_finel <- df_long %>% 
  filter(Vehicle_Type_Code %in% twoWheeled$Name) %>%
  select("COLLISION_ID","Vehicle_Type_Code")
person_df <- person_df[c("COLLISION_ID", "SAFETY_EQUIPMENT", "EMOTIONAL_STATUS")]

סינון רשומות עם ערכים ריקים ותיקון ערכים המופיעים כ- “helmet.” ל- “helmet”.


person_df<- person_df%>%filter_all(all_vars(. != "Unknown")) %>%
  filter(EMOTIONAL_STATUS != "Does Not Apply") %>%na.omit() 

person_df$SAFETY_EQUIPMENT <- gsub(".*helmet.*", "helmet", person_df$SAFETY_EQUIPMENT, ignore.case = TRUE)

מיזוג הנתונים

נבצע Inner Join בין הנתונים שנערכו על סמך עמודת ‘מס’ סידורי של התאונה’. מיזוג זה יאפשר לנו לאתר את נתוניו של אדם אשר היה מעורב בתאונה עם כלי דו-גלגלי. הסינון שבוצע מותיר אותנו עם תאונות בהן היה שימוש בקסדה (או לא) בלבד.


combined_data <- inner_join(Vehicle_df, crashes_df_finel, by = "COLLISION_ID") %>%
  inner_join(person_df, by = "COLLISION_ID") %>%
  filter(SAFETY_EQUIPMENT == "helmet" | SAFETY_EQUIPMENT == "None")

סטטיסטיקה תיאורית

בטבלה הבאה נוכל להבחין בפילוח של כמות המעורבים בתאונות בקרב כלים דו-גלגליים בהתבסס על השימוש בקסדה ומצבם הקוגניטיבי לאחר התאונה:


Apparent Death Conscious Incoherent Semiconscious Shock Unconscious
helmet 160 13906 74 202 429 162
None 132 21849 148 329 716 299

בגרף העמודות הבא נבחין ביחסים שבין השימוש לקסדה ואי שימוש בה בהתאם לסוגי המצבים הקוגניטיביים השונים:


בגרף העוגה הבא נבחין בהתפלגות המצבים הקוגניטיביים בקרב המעורבים. מאחר וכפי שהוצג לפני כן – ערך המצב הקוגניטיבי ‘בהכרה’ גבוה באופן ניכר מיתר המצבים – הוא הוסר מגרף זה על מנת לאפשר קריאה נוחה של כלל הנתונים.


ניתוח סטטיסטי

מבחן חי-בריבוע

נבצע מבחן חי-בריבוע לאי תלות בין ‘אמצעי הבטיחות’ (משתנה קטגורי) לבין ה’מצב הקוגניטיבי’ של המעורב (משתנה קטגורי). באמצעות מבחן זה נוכל להעריך האם קיים קשר מובהק בין שני המשתנים.


chi_square_result <- chisq.test(table(combined_data$SAFETY_EQUIPMENT, combined_data$EMOTIONAL_STATUS))
print(chi_square_result)
## 
##  Pearson's Chi-squared test
## 
## data:  table(combined_data$SAFETY_EQUIPMENT, combined_data$EMOTIONAL_STATUS)
## X-squared = 37.828, df = 5, p-value = 4.085e-07

נדחה את השערת האפס ונסיק כי קיים קשר מובהק בין ‘אמצעי הבטיחות’ לבין ‘המצב הקוגניטיבי’.

מבחן פירסון לקורלאציה

המרת המשתנים ‘אמצעי הבטיחות’ ו- ‘המצב הקוגניטיבי’ לפקטורים ושימושם כערכים מספריים על מנת שנוכל לבצע את בדיקת הקורלאציה:


combined_data$SAFETY_EQUIPMENT_numeric <- as.numeric(factor(combined_data$SAFETY_EQUIPMENT))
combined_data$EMOTIONAL_STATUS_numeric <- as.numeric(factor(combined_data$EMOTIONAL_STATUS))

pearson_correlation <- cor.test(combined_data$SAFETY_EQUIPMENT_numeric, combined_data$EMOTIONAL_STATUS_numeric)
print(pearson_correlation)
## 
##  Pearson's product-moment correlation
## 
## data:  combined_data$SAFETY_EQUIPMENT_numeric and combined_data$EMOTIONAL_STATUS_numeric
## t = 2.7365, df = 38404, p-value = 0.006213
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.003961876 0.023960328
## sample estimates:
##       cor 
## 0.0139625

נדחה את השערת האפס ונסיק כי קיימת קורלאציה בין ‘אמצעי הבטיחות’ לבין ‘המצב הקוגניטיבי’.


לסיכום, המבחנים הסטטיסטיים והוויזואליזציה מספקים תובנות אודות הקשר שבין שימוש בקסדה לבין המצב הקוגניטיבי של אנשים המעורבים בתאונות עם כלים דו-גלגליים. ניתוחים אלה מסייעים בהבנת דפוסים ומתאמים פוטנציאליים בתוך הנתונים, אשר יכולים לספק מידע נוסף לשיפור הבטיחות בדרכים.

האם קיים מתאם והבדל מובהק במספר התאונות בקרב כלים דו גלגליים בבריטניה בין השנים 2015 ל-2016 בסידור ע”פ חודשי השנה

לירז עזריה - 207167347

שאלת מחקר מתבססת על בסיסי הנתונים של העיר לידס בבריטניה בין השנים 2015, 2016:

אופן תהליך השימוש בבסיסי הנתונים

ראשית, בשתי הדאטות התמקדתי בתאונות דרכים של כלים דו גלגליים – אופנוע ואופניים.

# Focusing only in motorcycle / bike accidents:
is_motor_vehicle1 <- 
  as.integer(grepl("motor|motorcycle|MOTOR|bike|BIKE|cc|CC", data_2015$Type.of.Vehicle, ignore.case = TRUE))

data_2015 <- data.frame(data_2015, is_motor_vehicle1)
data_2015_relevant <- subset(data_2015, is_motor_vehicle1 == 1)


is_motor_vehicle1 <- 
  as.integer(grepl("motor|motorcycle|MOTOR|bike|BIKE|cc|CC", data_2016$Type.of.Vehicle, ignore.case = TRUE))

data_2016 <- data.frame(data_2016, is_motor_vehicle1)
data_2016_relevant <- subset(data_2016, is_motor_vehicle1 == 1)

שנית, חילצתי את מספרי חודשים ושמות חודשים עבור המיפוי של נתוני החודשים לגרף.

# Extracting Months for 2015 and 2016
data_2015_relevant$Month <- substr(data_2015_relevant$Accident.Date, 4, 6) 
data_2016_relevant$Month <- data_2016_relevant$Accident.Date %>% month()

מיזוג הטבלאות של הכמויות של תאונות בכל חודש עבור 2015 ו- 2016 ע”י פונקציית MERGE.

# Generating the main data table:
month_numbers <- seq(1,12)
month_names <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug",
                 "Sep", "Oct", "Nov", "Dec")

month_index <- data.frame(month_numbers, month_names)

table_2015_by_month <- data_2015_relevant %>%
  group_by(Month) %>% summarise("count_2015" = n()) %>% data.frame()

table_by_month <- merge(x = month_index, y = table_2015_by_month, 
      by.x = "month_names", by.y = "Month")


table_2016_by_month <- data_2016_relevant %>%
  group_by(Month) %>% summarise("count_2016" = n()) %>% data.frame()

table_by_month <- merge(x = table_by_month, y = table_2016_by_month, 
                        by.x = "month_numbers", by.y = "Month")

יצרתי טבלה שמציגה את כמות תאונות בסידור חודשי בשנת 2015.

month_numbers month_names count_2015 count_2016
1 Jan 10 13
2 Feb 7 14
3 Mar 17 12
4 Apr 10 19
5 May 22 17
6 Jun 16 23
7 Jul 16 20
8 Aug 30 19
9 Sep 28 12
10 Oct 24 15
11 Nov 14 16
12 Dec 9 9

תרשים BAR לכמות תאונות הדרכים בקרב כלים דו גלגליים בין השנים 2015 ל- 2016 בעיר לידס בבריטניה

לבדיקה

גרף צפיפות המתאר את כמות תאונות הדרכים בבריטניה בין השנים 2015, 2016 בכלים דו גלגליים

גרף תלת מימד המסכם את נתוני תאונות הדרכים בין השנים 2015 ו- 2016 בקרב כלים דו גלגליים

לבדיקה


ניתוחים סטטיסטים

בדיקת קורלציה (מתאם) בין מספר התאונות בשנת 2015 לשנת 2016 כלומר האם מס’ תאונות הדרכים בשנת 2015 מרמזת ומעניקה מידע מדויק עבור מס’ תאונות הדרכים בקרב כלים דו גלגליים לשנת 2016.

# Correlation test between count 2015 and count 2016:
print("Correlation between month of 2015 and 2016", quote = F)
## [1] Correlation between month of 2015 and 2016
cor(table_by_month$count_2015, table_by_month$count_2016)
## [1] 0.1713941
# Result: the correlation between the variables is positive but 
# quite low.


מתוצאות המבחן ניתן לראות כי הקורלציה חיובית אבל מעט נמוכה, מה שאומר שהקשר בין המשתנים חלש - אין מתאם בכמות תאונות הדרכים בכלים דו גלגליים בין השנים 2015 ל-2016.


יצרתי את הגרף שמציג את הקורלציה בין מספר התאונות בשנת 2015 לשנת 2016:


הקורלציה יצאה חיובית ונמוכה, זה אומר שישנה התאמה חיובית חלשה בין הנתונים במילים אחרות, קיים קשר חיובי בין המשתנים, אך הקשר לא חזק מה שמצביע על כך שאין התאמה בין כמות תאונות הדרכים בין השנים 2015 ל-2016 בבריטניה.


בהמשך, רציתי לבדוק האם קיים הבדל מובהק בין כמות תאונות הדרכים בקרב כלים דו גלגליים בין השנים 2015 ל- 2016 בעיר לידס שבבריטניה ולכן ביצעתי מבחן T מזווג.

הנחות המודל

• בדיקת התפלגות נורמלית של הנתונים ביצעתי ע”י מבחן שפירו.


print("Shapiro Wilks test for 2015:", quote=F)
## [1] Shapiro Wilks test for 2015:
shapiro.test(table_by_month$count_2015)
## 
##  Shapiro-Wilk normality test
## 
## data:  table_by_month$count_2015
## W = 0.933, p-value = 0.413
print("Shapiro Wilks test for 2016:", quote=F)
## [1] Shapiro Wilks test for 2016:
shapiro.test(table_by_month$count_2016)
## 
##  Shapiro-Wilk normality test
## 
## data:  table_by_month$count_2016
## W = 0.98154, p-value = 0.9891


תוצאת המבחן היא שאכן ישנה התפלגות נורמלית (קבלת 0H ודחיית 1H בר”מ 5%) של הנתונים גם בשנת 2015 וגם ב- 2016.
• לא ידועה תלות בין הנתונים
• גודל המדגם: קיבצנו ל-12 חודשים עבור כל אחת מהשנים, אז 12 חודשים בשנת 2015 ו- 12 חודשים ב- 2016.
בבסיסי הנתונים (המשתנים המקוריים) ישנם: 203 תצפיות עבור שנת 2015, ועבור 2016, 189 תצפיות.
• סוג המשתנים: מספר תאונות בחלוקה לפי חודש - משתנה רציף מסוג יחס-מנה.
• לא היה צורך לבדוק שוויון שונויות מכיוון שמדובר במבחן T מזווג (המבחן מדגיש את ההבדל בין הערכים באותו זמן או באותה קבוצה של נתונים, ולכן אין צורך לבדוק את השיוון של ההפרשים בין המדגמים)

מבחן T מזווג

t.test(x = table_by_month$count_2015, y = table_by_month$count_2016, paired = T)
## 
##  Paired t-test
## 
## data:  table_by_month$count_2015 and table_by_month$count_2016
## t = 0.50674, df = 11, p-value = 0.6223
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -3.900623  6.233956
## sample estimates:
## mean difference 
##        1.166667


לאחר ביצוע מבחן T מזווג ניתן היה להסיק באופן ברור כי אין הבדל בין כמות תאונות הדרכים בקרב כלים דו גלגליים בין השנים 2015 ל- 2016 כלומר אנו תומכים בהשערת האפס ונדחה את 1H בר”מ 5%.

האם יש הבדל בכמות כלי הרכב המעורבים התאונות בעיר לידס בבריטניה בין השנים 2015-2016

אופן ייצוא הדאטה

לקחתי את דאטה בריטניה 2015, דאטה בריטניה 2016 ודאטה ניו-יורק ואיחדתי אותה באמצעות JOIN עם השדה HOUR. אין שדות ריקים כלל בגלל שיש בכל הבסיסי נתונים מידע על HOUR. המרתי את השעות משעה 2312 בבריטניה ל23, ובניו יורק מהשעה “04:00” ל-4. כך יש לי שדות מ0-23.


פילטור של בסיסי הנתונים של UK2015 UK2016 לפי כלי רכב דו גלגליים.

is_motor_vehicle1 <- 
  as.integer(grepl("motor|motorcycle|MOTOR|bike|BIKE|cc|CC", data_2015$Type.of.Vehicle, ignore.case = TRUE))

data_2015 <- data.frame(data_2015, is_motor_vehicle1)
data_2015_relevant <- subset(data_2015, is_motor_vehicle1 == 1)


is_motor_vehicle1 <- 
  as.integer(grepl("motor|motorcycle|MOTOR|bike|BIKE|cc|CC", data_2016$Type.of.Vehicle, ignore.case = TRUE))

data_2016 <- data.frame(data_2016, is_motor_vehicle1)
data_2016_relevant <- subset(data_2016, is_motor_vehicle1 == 1)


פילטור של בסיס הנתונים של NY לפי כלי רכב דו גלגליים.

is_motor_vehicle <- 
  as.integer(grepl("motor", data_ny$VEHICLE.TYPE.CODE.1, ignore.case = TRUE))

is_motor_vehicle1 <- 
  as.integer(grepl("motor|motorcycle|MOTOR|bike|BIKE", data_ny$VEHICLE.TYPE.CODE.1, ignore.case = TRUE))


data_vector <- data.frame(is_motor_vehicle, is_motor_vehicle1)

data_ny <- data.frame(data_ny, is_motor_vehicle1) 

data_ny_relevant <- subset(data_ny, is_motor_vehicle1 == 1)


תחילה, סידרתי בסיס נתונים חדש לפי השעות ביום וערכתי רגרסיה לינארית, על השכיחויות של התאונות דרכים לפי שעה. מכיוון שסדרי הגודל של הבסיסי נתונים לא שווים בגודלם אז החלטתי להשאיר את בסיס הנתונים הזה ולעשות איתו סטטיסטיקה תיאורית עם נתונים מנורמלים על ידי scale.
בשאלת המחקר שלי בחרתי לבדוק את ההבדל בכמות כלי הרכב המעורבים בתאונות בבריטניה, כאשר הטווח הוא עד 4 רכבים, יותר מ-4 רכבים ייכנסו ל-4. גם בטבלה זו אין לי ערכים ריקים.


גרף קורלציה ראשון בין מספר התאונות לפי שעה בין UK לNY

לבדיקה


• בין count_ny ו־ uk_2015 יש קורלציה חיובית חזקה בערך של 0.794. משמע ככל שמספר התאונות בניו יורק עולה, כך מספר התאונות בעיר לידס בבריטניה עולה בשנת 2015.
• בין count_ny ו־uk_2016 יש קורלציה חיובית חזקה בערך של 0.743. זו משמע ככל שמספר התאונות בניו יורק עולה, כך מספר התאונות בעיר לידס בבריטניה עולה בשנת 2016.
• בין uk_2015 ל־ uk_2016 יש קורלציה חיובית חזקה בערך של 0.863. משמע ככל שמספר התאונות בעיר לידס בבריטניה בשנת 2015, כך מספר התאונות בעיר לידס בבריטניה עולה בשנת 2016.


טבלה חדשה עם מניפולציה על השעות ונרמול נתונים

Hour Scaled_NY Scaled_UK_2015 Scaled_UK_2016
0 -0.1317 -1.0976 -1.2122
1 -1.0377 -0.9678 -1.2122
2 -1.2604 -1.0976 -1.0583
3 -1.3552 -1.0976 -1.0583
4 -1.4061 -1.0976 -1.2122
5 -1.4529 -1.0976 -1.0583
6 -1.1174 -0.5785 -0.4426
7 -0.7902 0.8489 0.6350
8 -0.2774 0.0703 0.6350
9 -0.2499 -0.1892 -0.5965
10 -0.3365 -0.1892 -0.4426
11 -0.0162 -0.0595 0.0192
12 0.2436 0.4596 0.4810
13 0.4155 1.1084 2.1743
14 0.9351 -0.1892 1.0968
15 1.0781 0.3298 -0.1347
16 1.4479 1.7572 0.4810
17 1.7779 2.6655 1.8665
18 1.5799 1.3679 1.4047
19 1.1743 0.4596 0.9429
20 0.7605 0.2000 0.0192
21 0.3192 -0.4488 -0.2886
22 0.1144 -0.3190 -0.1347
23 -0.4149 -0.8380 -0.9044

גרף 3D המציג את השכיחויות של תאונות ברכבים דו גלגליים לפי השעה ביום

לבדיקה

גרף שכיחות תאונות דרכים ברכבים דו גלגליים בUK- לפי השעה ביום

לבדיקה


ניתן לראות בגרף זה כי מרבית התאונות בשנת 2015 וגם בשנת 2016 מתרחשות בשעות הצהריים בהתאמה.


גרף קורלציה שני בין מספר כלי הרכב המעורבים בתאונה לבין מספר התאונות


נראה כי יש קשר חזק מאוד.

גרף שכיחויות בכמות כלי הרכב המעורבים בתאונה בUK



ניתן לראות בגרף שכיחויות בכמות כלי הרכב המעורבים בתאונה כי אין הבדל משמעותי בין שנים 2015 ו-2016.


השערות


H0: אין הבדל משמעותי במספר הרכבים המעורבים בתאונות בין שנת 2015 לשנת 2016.
H1: קיים הבדל משמעותי במספר הרכבים המעורבים בתאונות בין שנת 2015 לשנת 2016.


הנחות מודל T-TEST

נורמליות

shapiro.test(vec_2015)
## 
##  Shapiro-Wilk normality test
## 
## data:  vec_2015
## W = 0.67541, p-value < 2.2e-16
shapiro.test(vec_2016)
## 
##  Shapiro-Wilk normality test
## 
## data:  vec_2016
## W = 0.72839, p-value < 2.2e-16


מכיוון ששני הבדיקות הביאו לתוצאות משמעותיות, אנו יכולים להניח כי הנתונים עבור שני הווקטורים אינם מתפלגים בדרך נורמלית. עם זאת, מאחר וגודל התצפיות עבור שני הווקטורים גדול מאוד - גדול מ־30 , אנו יכולים להתעלם מכך ולהשתמש עדיין בבדיקת t ולהניח נורמליות של התוצאות. עם זאת, נשתמש בשתי הבדיקות ונסקור את המסקנות על פי שתיהן.

שיוויון שונויות

## Test of variances for both results for 2015 and 2016:
var.test(vec_2015, vec_2016)
## 
##  F test to compare two variances
## 
## data:  vec_2015 and vec_2016
## F = 0.71111, num df = 202, denom df = 188, p-value = 0.01747
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5361349 0.9418266
## sample estimates:
## ratio of variances 
##          0.7111085


נמצא כי השונויות שונות באופן משמעותי - לכן אנו צריכים להניח את זה בעת השימוש בניתוח ה-t-test.

T-Test

t.test(x = data_2015_relevant$Number.of.Vehicles, 
       y = data_2016_relevant$Number.of.Vehicles, var.equal = F)
## 
##  Welch Two Sample t-test
## 
## data:  data_2015_relevant$Number.of.Vehicles and data_2016_relevant$Number.of.Vehicles
## t = -0.84809, df = 368.86, p-value = 0.3969
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.15863568  0.06303269
## sample estimates:
## mean of x mean of y 
##  1.788177  1.835979


לפי ניתוח ה-t-test - לא נמצאה הבדל משמעותי בין שנת 2015 ו-2016, זאת אומרת שאין לפסול את השערת האפס. (H0)

Wilcox Test - א-פרמטרי

wilcox.test(data_2015_relevant$Number.of.Vehicles, 
            data_2016_relevant$Number.of.Vehicles, 
            alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data_2015_relevant$Number.of.Vehicles and data_2016_relevant$Number.of.Vehicles
## W = 18727, p-value = 0.6183
## alternative hypothesis: true location shift is not equal to 0


הבדיקה של וילקוקסון הביאה לתוצאה דומה: אין הבדל משמעותי במספר הרכבים המעורבים בין שנת 2015 ל-2016.
התוצאה הדומה הגיעה מכך שגודל המדגם שלנו גדול מאוד ולכן אין צורך בבדיקת הנחות במקרה זה.