Econometrics: Panel Data Analysis
Read-only preview of the saved R notebook. Code is collapsed, and figures and results are from the original run. Nothing is executed here.
Panel Data Analysis – Guns Dataset¶
בחלק זה נבצע ניתוח נתוני פאנל באמצעות מודלי Fixed Effects (FE) ו-Random Effects (RE).
מטרת הניתוח היא לבחון את הקשר בין שיעור הפשיעה האלימה במדינות ארה"ב לבין משתנים דמוגרפיים, כלכליים ומשתני מדיניות לאורך זמן, ולבחור בין שני המודלים באמצעות מבחן Hausman.
לאחר בחירת המודל נבחן את שאריותיו ונבדוק האם קיימות בעיות במפרט המודל.
Show code, cell 1
# Packages
library(AER)
library(plm)
library(corrplot)
library(lmtest)
library(tseries)
library(nlme)
Show code, cell 2
# Load the dataset
data("Guns", package = "AER")
guns <- Guns
1. Data Description¶
הניתוח מבוסס על מאגר הנתונים Guns מספריית AER ב-R.
זהו מאגר נתוני פאנל המתאר את 50 מדינות ארה"ב ואת מחוז קולומביה לאורך השנים 1977–1999. כל תצפית מייצגת מדינה מסוימת בשנה מסוימת.
המשתנה המוסבר בניתוח הוא violent, המייצג את שיעור מקרי הפשיעה האלימה ל-100,000 תושבים. מטרת הניתוח היא לבחון כיצד שיעור זה קשור למאפיינים דמוגרפיים וכלכליים של המדינה ולמשתנה law, המציין האם חוק Shall-Carry היה בתוקף במדינה באותה שנה.
המשתנים במאגר הם:
- state – המדינה.
- year – שנת התצפית.
- violent – שיעור מקרי פשיעה אלימה ל-100,000 תושבים.
- murder – שיעור מקרי רצח ל-100,000 תושבים.
- robbery – שיעור מקרי שוד ל-100,000 תושבים.
- prisoners – שיעור האסירים שנידונו למאסר ל-100,000 תושבים, בשנה הקודמת.
- afam – אחוז האוכלוסייה האפרו-אמריקאית בגילאי 10–64.
- cauc – אחוז האוכלוסייה הקווקזית בגילאי 10–64.
- male – אחוז הגברים בגילאי 10–29.
- population – אוכלוסיית המדינה, במיליוני תושבים.
- income – הכנסה אישית ריאלית לנפש.
- density – צפיפות אוכלוסייה לשטח.
- law – האם חוק Shall-Carry היה בתוקף במדינה ובשנה הנתונה.
Show code, cell 3
# Inspect the data
head(guns)
str(guns)
summary(guns)
dim(guns)
2. Data Preparation and Panel Structure¶
בשלב זה בחנו את מבנה מאגר הנתונים ואת הסטטיסטיקה התיאורית של המשתנים, ובדקנו האם קיימים ערכים חסרים.
בנוסף, הגדרנו את state כממד הרוחבי ואת year כממד הזמן.
במקרה שלנו, כל מדינה היא יחידה נפרדת וכל שנה היא תקופת זמן.
המשתנה year כבר מוגדר במאגר כמשתנה מסוג factor בעל 23 רמות. לכן כאשר הוא נכלל במודל, R יוצר באופן אוטומטי משתני דמה לשנים, כאשר שנת 1977 משמשת כקטגוריית הבסיס.
Show code, cell 4
# Check for missing values
anyNA(guns)
colSums(is.na(guns))
# Create a panel data frame
guns_panel <- pdata.frame(
guns,
index = c("state", "year")
)
# Check the panel structure
pdim(guns_panel)
המאגר כולל 1,173 תצפיות ו-13 משתנים, המתקבלים מ-51 מדינות הנצפות במשך 23 שנים.
לא נמצאו ערכים חסרים במאגר.
בנוסף, נמצא כי מדובר ב- Balanced Panel, כלומר לכל אחת מ-51 המדינות קיימת תצפית בכל אחת מ-23 השנים.
לכן ניתן להמשיך לניתוח ללא צורך בהסרת תצפיות עקב ערכים חסרים או חוסר איזון במבנה הפאנל.
3. Exploratory Data Analysis¶
לפני הגדרת מודלי Fixed Effects ו-Random Effects, בחנו את ההתפלגויות של המשתנה המוסבר ושל המשתנים הרציפים שעשויים להיכלל במודל.
מטרת הבדיקה היא לזהות משתנים בעלי התפלגות מוטה מאוד או פערים גדולים בסדרי הגודל, ולבחון האם טרנספורמציית log עשויה להתאים להם.
בשלב זה איננו דורשים שהמשתנים המסבירים עצמם יתפלגו נורמלית. ההיסטוגרמות משמשות בעיקר לבחינת צורת ההתפלגות ולבחינת הצורך בטרנספורמציה.
Show code, cell 5
options(repr.plot.width = 16, repr.plot.height = 12)
par(mfrow = c(2, 2))
hist(guns$violent,
main = "Violent Crime Rate",
xlab = "violent")
hist(guns$prisoners,
main = "Prisoners",
xlab = "prisoners")
hist(guns$population,
main = "Population",
xlab = "population")
hist(guns$income,
main = "Income",
xlab = "income")
par(mfrow = c(1, 1))
Show code, cell 6
options(repr.plot.width = 16, repr.plot.height = 12)
par(mfrow = c(2, 2))
hist(guns$density,
main = "Population Density",
xlab = "density")
hist(guns$afam,
main = "African-American Population",
xlab = "afam")
hist(guns$cauc,
main = "Caucasian Population",
xlab = "cauc")
hist(guns$male,
main = "Young Male Population",
xlab = "male")
par(mfrow = c(1, 1))
4. Log Transformation¶
מההיסטוגרמות ניתן לראות כי violent, prisoners, population ו-density מאופיינים בהטיה משמעותית ימינה.
לכן נבצע עבור משתנים אלה טרנספורמציית log.
לבחירה זו יש שתי מטרות. מבחינה סטטיסטית, הלוגריתם מצמצם את השפעתם של ערכים גבוהים מאוד ומקרב את המשתנים לסקאלה מאוזנת יותר. מבחינה כלכלית, הוא גם מאפשר לפרש קשרים יחסיים, כאשר גם המשתנה המוסבר וגם מסביר מסוים נמצאים בלוגריתם, המקדם של אותו מסביר ניתן לפרש כגמישות בקירוב - אחוז השינוי ב-violent הקשור לעלייה של 1% במסביר, בהינתן יתר המשתנים.
לאחר הטרנספורמציה נבחן מחדש את ההתפלגויות כדי לוודא שהשינוי אכן שיפר את צורתן.
המשתנים income, afam, cauc ו-male יישארו בשלב זה בסקאלה המקורית, משום שלא נצפתה עבורם הטיה המצדיקה טרנספורמציה דומה.
Show code, cell 7
# Log transformations for right-skewed variables
guns$log_violent <- log(guns$violent)
guns$log_prisoners <- log(guns$prisoners)
guns$log_population <- log(guns$population)
guns$log_density <- log(guns$density)
Show code, cell 8
options(repr.plot.width = 16, repr.plot.height = 12)
par(mfrow = c(2, 2))
hist(guns$log_violent,
main = "Log Violent Crime Rate",
xlab = "log(violent)")
hist(guns$log_prisoners,
main = "Log Prisoners",
xlab = "log(prisoners)")
hist(guns$log_population,
main = "Log Population",
xlab = "log(population)")
hist(guns$log_density,
main = "Log Population Density",
xlab = "log(density)")
par(mfrow = c(1, 1))
לאחר ביצוע טרנספורמציית log ניתן לראות שיפור ברור בהתפלגויות של ארבעת המשתנים.
ההטיה החזקה ימינה של violent ו-prisoners הצטמצמה משמעותית, ושתי ההתפלגויות הפכו סימטריות יותר.
גם עבור population ו-density הטרנספורמציה צמצמה באופן משמעותי את הפערים בין הערכים ואת ההטיה שנצפתה בנתונים המקוריים.
לכן בהמשך הניתוח נשתמש ב- log(violent) כמשתנה המוסבר, וב- log(prisoners), log(population) ו-log(density) כגרסאות של המשתנים המסבירים המתאימים.
5. Correlation and Multicollinearity¶
לפני הגדרת מודלי הפאנל נבחן את הקשרים בין המשתנים המסבירים הרציפים.
מטרת הבדיקה היא לזהות זוגות של משתנים בעלי מתאם גבוה מאוד, אשר עלולים ליצור בעיית multicollinearity ולהקשות על הפרדת ההשפעה של כל אחד מהם במודל.
המשתנים murder ו-robbery אינם נכללים כמסבירים, מכיוון שהם מתארים סוגים של פשיעה הנכללים במדד הפשיעה האלימה שאותו אנו מנסים להסביר.
Show code, cell 9
# Correlation matrix
cor_data <- guns[, c(
"log_prisoners",
"afam",
"cauc",
"male",
"log_population",
"income",
"log_density"
)]
cor_matrix <- cor(cor_data)
round(cor_matrix, 2)
Show code, cell 10
options(repr.plot.width = 14, repr.plot.height = 12)
corrplot(
cor_matrix,
method = "number",
type = "upper",
tl.cex = 1.1,
number.cex = 1,
number.digits = 2
)
מטריצת המתאמים מצביעה על מתאם שלילי גבוה מאוד בין afam ל-cauc, כאשר מקדם המתאם הוא -0.98.
כדי לשמור על מפרט חסכוני ולהימנע מ- multicollinearity נשאיר את afam כמשתנה המייצג היבט זה של ההרכב הדמוגרפי ונשמיט את cauc. הבחירה אינה טענה שלפיה אחד המשתנים חשוב כלכלית יותר מהשני, היא נועדה למנוע ייצוג כפול של מידע כמעט זהה.
שאר המתאמים בין המשתנים המסבירים הם מתונים יותר, ולכן בשלב זה אין הצדקה להסיר משתנים נוספים על סמך מטריצת המתאמים בלבד.
6. Model Specification¶
לאחר שלבי הכנת הנתונים, הטרנספורמציות ובדיקת המתאמים, נגדיר את המפרט הראשוני של מודלי הפאנל.
המשתנה המוסבר הוא log_violent.
המשתנים המסבירים כוללים את משתנה המדיניות law, את המשתנים הדמוגרפיים afam ו-male, את המשתנים הכלכליים והמבניים log_prisoners, log_population, income ו-log_density, וכן משתני שנה לצורך שליטה בשינויים משותפים לכל המדינות לאורך זמן.
המשתנים murder ו-robbery אינם נכללים במודל משום שהם מתארים סוגי פשיעה הנכללים במדד violent שאותו אנו מנסים להסביר.
המשתנה cauc אינו נכלל בעקבות המתאם הגבוה עם afam.
Initial Model Equation¶
בהתאם לבחירת המשתנים ולטרנספורמציות שבוצעו, המפרט הראשוני של המודל הוא:
כאשר:
- $i$ מייצג מדינה.
- $t$ מייצג שנה.
- $violent_{it}$ הוא שיעור הפשיעה האלימה במדינה $i$ בשנה $t$.
- $D_t$ הם משתני דמה לשנים, כאשר שנת 1977 משמשת כשנת הבסיס.
- $\beta_1,\ldots,\beta_7$ הם המקדמים של המשתנים המסבירים.
- $\gamma_t$ מייצגים את השפעות השנים.
- $\varepsilon_{it}$ הוא רכיב השגיאה.
Show code, cell 11
# Initial panel model specification
formula_initial <- log_violent ~
law +
log_prisoners +
afam +
male +
log_population +
income +
log_density +
year
7. Fixed Effects and Random Effects Models¶
כעת נאמוד את שני מודלי הפאנל על בסיס אותו מפרט ראשוני.
מודל Fixed Effects (FE) מאפשר לאפקטים הבלתי נצפים והקבועים של כל מדינה להיות מתואמים עם המשתנים המסבירים. האמידה מבוססת על השינויים המתרחשים בתוך אותה מדינה לאורך זמן.
מודל Random Effects (RE) מניח שהאפקט הייחודי של כל מדינה אינו מתואם עם המשתנים המסבירים.
בשלב זה נאמוד את שני המודלים על אותו מפרט ונבחן את תוצאותיהם.
אם תזוהה בעיה מבנית במפרט המשותף לשני המודלים, כגון multicollinearity חזקה, נטפל בה לפני ההשוואה.
לאחר קביעת מפרט משותף תקין נשתמש במבחן Hausman לבחירה בין FE ל-RE. סינון המשתנים לפי רמת מובהקות יתבצע לאחר מכן במודל שנבחר.
Fixed Effects Model¶
$$ \ln(violent_{it}) = \alpha_i + X_{it}'\beta + \sum_{t=1978}^{1999}\gamma_t D_t + \varepsilon_{it} $$במודל Fixed Effects לכל מדינה קיים אפקט קבוע משלה, $\alpha_i$. כך ניתן לשלוט במאפיינים בלתי נצפים של המדינה שאינם משתנים לאורך זמן, גם כאשר הם מתואמים עם המשתנים המסבירים.
Random Effects Model¶
$$ \ln(violent_{it}) = \alpha + X_{it}'\beta + \sum_{t=1978}^{1999}\gamma_t D_t + \mu_i + \varepsilon_{it} $$במודל Random Effects ההבדל הייחודי בין המדינות מיוצג באמצעות $\mu_i$, שהוא רכיב אקראי של המדינה.
ההנחה המרכזית של מודל זה היא שהאפקט הייחודי למדינה אינו מתואם עם המשתנים המסבירים:
מבחן Hausman ישמש בהמשך לבחינת התאמת הנחה זו ולבחירה בין שני המודלים.
Show code, cell 12
# Fixed Effects model
fe_initial <- plm(
formula_initial,
data = guns,
index = c("state", "year"),
model = "within"
)
summary(fe_initial)
Show code, cell 13
# Random Effects model
re_initial <- plm(
formula_initial,
data = guns,
index = c("state", "year"),
model = "random"
)
summary(re_initial)
8. Specification Refinement Before Hausman¶
תוצאות המודלים הראשוניים הראו Standard Errors גדולים במיוחד עבור log_population ו-log_density במודל Fixed Effects.
מאחר שאמידת FE מבוססת על השינויים בתוך כל מדינה לאורך זמן, נבדוק האם שני המשתנים כמעט נעים יחד בתוך אותה מדינה.
מטרת שלב זה אינה לסנן משתנים לפי מובהקות, אלא לוודא שהמפרט המשותף שישמש להשוואת FE ו-RE אינו סובל מבעיית multicollinearity .
Show code, cell 14
# Check within-state correlation between population and density
log_population_within <- guns$log_population -
ave(guns$log_population, guns$state)
log_density_within <- guns$log_density -
ave(guns$log_density, guns$state)
cor(log_population_within, log_density_within)
המתאם בין השינויים בתוך המדינה של log_population ושל log_density עומד על 0.9992, ומצביע על כך ששני המשתנים מכילים כמעט אותו מידע מבחינת השונות לאורך זמן.
שטחה של כל מדינה כמעט קבוע לאורך תקופת המדגם. לכן, בתוך אותה מדינה, שינוי בצפיפות משקף כמעט באופן מלא שינוי בגודל האוכלוסייה.
מאחר שמודל Fixed Effects מזהה את המקדמים מתוך השינויים בתוך כל מדינה לאורך זמן, הכללת שני המשתנים יחד יוצרת multicollinearity כמעט מושלמת. תוצאה זו מסבירה את ה- Standard Errors הגדולים שהתקבלו עבורם באמידה הראשונית.
לפיכך, נסיר את log_density ונשאיר את log_population. בחירה זו מונעת כפילות מידע ומאפשרת לשמור במודל מדד ישיר וברור לגודל האוכלוסייה של המדינה.
לאחר הסרת log_density ניתן לאמוד מחדש את מודלי Fixed Effects ו-Random Effects על בסיס מפרט משותף ויציב יותר, ולאחר מכן להשוות ביניהם באמצעות מבחן Hausman.
Show code, cell 15
# Reduced model specification after removing within multicollinearity
formula_reduced <- log_violent ~
law +
log_prisoners +
afam +
male +
log_population +
income +
year
Show code, cell 16
# Fixed Effects reduced model
fe_reduced <- plm(
formula_reduced,
data = guns,
index = c("state", "year"),
model = "within"
)
summary(fe_reduced)
Show code, cell 17
# Random Effects reduced model
re_reduced <- plm(
formula_reduced,
data = guns,
index = c("state", "year"),
model = "random"
)
summary(re_reduced)
לאחר הסרת log_density הרצנו מחדש את מודלי Fixed Effects ו-Random Effects על אותו מפרט מצומצם.
ה- Standard Errors החריגים שנצפו קודם עבור משתני האוכלוסייה והצפיפות נעלמו, ולכן בעיית ה- multicollinearity טופלה.
מפרט משותף זה ישמש כעת לביצוע מבחן Hausman. בשלב זה איננו מבצעים עדיין סינון לפי מובהקות. סינון כזה יתבצע לאחר בחירת המודל.
Reduced Model Equation¶
לאחר הסרת log_density עקב ה- multicollinearity החזקה עם log_population בתוך המדינות, המשוואה ששימשה להשוואת מודלי FE ו-RE היא:
9. Hausman Test and Model Selection¶
לאחר קביעת מפרט משותף וסופי לשני המודלים, נשתמש במבחן Hausman כדי לבחור בין Fixed Effects ל-Random Effects.
השערת האפס היא שמודל Random Effects מתאים, כלומר האפקטים הייחודיים של המדינות אינם מתואמים עם המשתנים המסבירים.
אם ערך ה- p-value קטן מ-0.05, נדחה את השערת האפס ונעדיף את מודל Fixed Effects. אחרת, נעדיף את מודל Random Effects.
Show code, cell 18
# Hausman test
hausman_test <- phtest(fe_reduced, re_reduced)
hausman_test
במבחן Hausman התקבל ערך p-value = 0.0892.
מכיוון שערך זה גדול מרמת המובהקות של 5%, איננו דוחים את השערת האפס של המבחן.
לכן אין עדות מספקת לכך שהאפקטים הייחודיים של המדינות מתואמים עם המשתנים המסבירים, ובהתאם נבחר במודל Random Effects כמודל המתאים להמשך הניתוח.
מכאן והלאה בדיקות השאריות ותיקוני המודל, במידת הצורך, יבוצעו על מודל Random Effects שנבחר.
10. Variable Selection in the Selected Model¶
לאחר שמבחן Hausman הוביל לבחירת מודל Random Effects, בחנו את מובהקות המשתנים המסבירים במודל הנבחר.
המשתנה log_prisoners לא נמצא מובהק, עם p-value = 0.672, ולכן הוסר מהמפרט.
המשתנים afam, male, log_population ו-income נמצאו מובהקים ברמת 5%.
המשתנה law התקבל עם p-value = 0.084. מושאר במכוון במפרט למרות מובהקות גבולית, משום שהוא מהווה את מוקד שאלת המחקר המרכזית, ובלעדיו המודל מאבד את משמעותו.
השארת המשתנה מאפשרת לבחון באופן ישיר האם קיומו של חוק Shall-Carry קשור לשיעור הפשיעה האלימה, גם אם בסופו של דבר ההשפעה אינה מובהקת ברמת 5%.
משתני השנה נשארים במודל כמשתני בקרה להשפעות משותפות לכל המדינות לאורך זמן.
הסינון יתבצע באופן הדרגתי: בכל שלב יוסר המסביר שאינו מובהק ובעל ערך ה- p-value הגבוה ביותר, ולאחר מכן המודל ייאמד מחדש.
Show code, cell 19
# Random Effects model after removing log_prisoners
formula_step1 <- log_violent ~
law +
afam +
male +
log_population +
income +
year
re_step1 <- plm(
formula_step1,
data = guns,
index = c("state", "year"),
model = "random"
)
summary(re_step1)
במודל המעודכן, income הוא המשתנה המסביר בעל ערך ה- p-value הגבוה ביותר (0.191).
לכן, במסגרת תהליך סינון המשתנים, נסיר בשלב הבא את income ונאמוד את המודל מחדש.
Show code, cell 20
# Random Effects model after removing income
formula_step2 <- log_violent ~
law +
afam +
male +
log_population +
year
re_step2 <- plm(
formula_step2,
data = guns,
index = c("state", "year"),
model = "random"
)
summary(re_step2)
במודל המעודכן, log_population הוא המשתנה המסביר בעל ערך ה- p-value הגבוה ביותר (0.231), ולכן אינו מובהק ברמת 5%.
לעומתו, law נמצא קרוב מאוד לרמת מובהקות של 5%, afam מציג מובהקות חלשה ברמת 10%, ו-male מובהק מאוד.
לכן נסיר בשלב זה רק את log_population ונאמוד מחדש את המודל לפני קבלת החלטה נוספת.
Show code, cell 21
# Random Effects model after removing log_population
formula_step3 <- log_violent ~
law +
afam +
male +
year
re_step3 <- plm(
formula_step3,
data = guns,
index = c("state", "year"),
model = "random"
)
summary(re_step3)
המשתנה afam נותר לא מובהק, עם p-value = 0.132, ולכן יוסר מהמפרט.
המשתנה male נמצא מובהק מאוד, ואילו law מציג מובהקות חלשה ברמת 10% ונשאר במודל בשל תפקידו כמשתנה המדיניות המרכזי של הניתוח.
Show code, cell 22
# Final Random Effects model
formula_final <- log_violent ~
law +
male +
year
re_final <- plm(
formula_final,
data = guns,
index = c("state", "year"),
model = "random"
)
summary(re_final)
לאחר תהליך סינון המשתנים התקבל המודל הסופי.
המשתנה male נמצא מובהק מאוד (p-value < 0.001).
המשתנה law התקבל עם p-value = 0.056. לכן הוא אינו מובהק ברמת 5%, אך מציג מובהקות חלשה ברמת 10%. המשתנה נשמר במודל בשל תפקידו כמשתנה המדיניות המרכזי של הניתוח.
משתני השנה נמצאו ברובם המכריע מובהקים מאוד ונשמרו במודל כדי לשלוט בשינויים המשותפים לכל המדינות לאורך זמן.
למודל הסופי התקבל R² = 0.399.
מאחר שבשלב זה טרם נבדק מבנה השגיאה, המסקנות הסופיות לגבי מובהקות המקדמים ייקבעו לאחר ביצוע בדיקות השאריות והתיקון המתאים במידת הצורך.
מודל זה ייחשב למודל Random Effects הסופי, ועליו נבצע את בדיקות השגיאה.
Final Model Equation¶
לאחר בחירת מודל Random Effects וסיום תהליך סינון המשתנים, המפרט הסופי הוא:
במשוואה הסופית:
- $law_{it}$ מציין האם חוק Shall-Carry היה בתוקף במדינה ובשנה הנתונות.
- $male_{it}$ הוא אחוז הגברים בגילאי 10–29.
- $\gamma_t$ מייצגים השפעות שנה ביחס לשנת הבסיס 1977.
- $\mu_i$ הוא האפקט האקראי הייחודי למדינה.
- $\varepsilon_{it}$ הוא רכיב השגיאה המשתנה בין מדינה ושנה.
מודל זה הוא המודל שעליו יבוצעו בדיקות השאריות וה- misspecification.
11. Residual Diagnostics¶
לאחר קביעת מודל Random Effects הסופי, נבחן את שאריות המודל כדי לבדוק האם קיימות בעיות במבנה השגיאה.
נבחן תחילה את פיזור השאריות ביחס לערכים החזויים. לאחר מכן נבדוק heteroskedasticity ואת נורמליות השאריות.
Show code, cell 23
# Residuals vs Fitted - final Random Effects model
res_final <- as.numeric(resid(re_final))
fit_final <- as.numeric(fitted(re_final))
options(repr.plot.width = 12, repr.plot.height = 7)
plot(
fit_final,
res_final,
main = "Residuals vs Fitted - Final Random Effects Model",
xlab = "Fitted Values",
ylab = "Residuals",
pch = 20
)
abline(h = 0, col = "red")
גרף Residuals vs Fitted מראה כי השאריות מפוזרות באופן כללי סביב אפס, ללא מגמה או תבנית לא-ליניארית ברורה.
עם זאת, מידת הפיזור של השאריות אינה נראית אחידה לחלוטין לאורך כל טווח הערכים החזויים.
לכן לא ניתן להסיק מהגרף בלבד כי שונות השגיאה קבועה, ונבצע מבחן פורמלי ל- heteroskedasticity.
Heteroskedasticity Test¶
נשתמש במבחן Breusch–Pagan כדי לבדוק האם שונות השגיאות קבועה.
השערת האפס היא כי קיימת שונות קבועה בשגיאות:
לעומת ההשערה החלופית שלפיה שונות השגיאות אינה קבועה.
אם ערך ה- p-value קטן מ-0.05, נדחה את השערת האפס ונקבע שקיימת עדות ל- heteroskedasticity.
Show code, cell 24
# Breusch-Pagan heteroskedasticity test
hetero_test <- bptest(
re_final,
varformula = ~ fitted(re_final) + I(fitted(re_final)^2),
studentize = FALSE
)
hetero_test
במבחן Breusch–Pagan התקבל p-value = 0.0114.
מכיוון שערך זה קטן מרמת המובהקות של 5%, אנו דוחים את השערת האפס של שונות קבועה בשגיאות.
לכן קיימת עדות סטטיסטית ל- heteroskedasticity בשאריות מודל Random Effects הסופי.
ממצא זה מצביע על misspecification במבנה השגיאה, ולכן לאחר השלמת בדיקות השאריות ננסה לתקן את המודל.
Normality of Residuals¶
בנוסף לבדיקת שונות השגיאות, נבחן האם שאריות המודל מתפלגות בקירוב נורמלית.
הבדיקה תתבצע באופן גרפי באמצעות Q-Q Plot והיסטוגרמת השאריות, ובאופן פורמלי באמצעות מבחן Jarque–Bera.
במבחן JB השערת האפס היא שהשאריות מתפלגות נורמלית:
Show code, cell 25
# Normality diagnostics - final Random Effects model
options(repr.plot.width = 14, repr.plot.height = 6)
par(mfrow = c(1, 2))
# Q-Q plot
qqnorm(
res_final,
main = "Q-Q Plot - Final Random Effects Residuals"
)
qqline(res_final, col = "red")
# Histogram
hist(
res_final,
breaks = 30,
freq = FALSE,
main = "Histogram - Final Random Effects Residuals",
xlab = "Residuals"
)
lines(density(res_final), lwd = 2)
curve(
dnorm(
x,
mean = mean(res_final),
sd = sd(res_final)
),
col = "red",
lwd = 2,
add = TRUE
)
par(mfrow = c(1, 1))
Show code, cell 26
# Jarque-Bera normality test
jb_test <- jarque.bera.test(res_final)
jb_test
גרף Q-Q מראה כי מרבית השאריות במרכז ההתפלגות נמצאות קרוב לקו התאורטי, אך קיימות סטיות ברורות בזנבות.
גם ההיסטוגרמה מצביעה על התפלגות דמוית נורמלית במרכז, אך עם זנבות שאינם מתאימים באופן מלא להתפלגות נורמלית.
במבחן Jarque–Bera התקבל p-value = 1.506×10⁻¹⁰, ולכן אנו דוחים את השערת האפס של נורמליות השאריות.
בשילוב עם תוצאת מבחן Breusch–Pagan, המסקנה היא כי במודל הסופי קיימת heteroskedasticity וכן סטייה מנורמליות השאריות.
ננסה כעת לתקן את בעיית מבנה השגיאה.
12. Model Correction¶
בדיקות השאריות הצביעו על heteroskedasticity ועל סטייה מנורמליות.
מודל Random Effects שנבחר באמצעות מבחן Hausman נשאר מודל הפאנל המרכזי של הניתוח.
תחילה נתקן את ההסקה הסטטיסטית באמצעות robust standard errors, אשר שומרים על אומדי מודל RE אך מתקנים את אומדני אי-הוודאות תחת שונות לא קבועה.
לאחר מכן נבצע גם ניסיון FGLS כדי לבחון האם ניתן לתקן את מבנה השגיאה עצמו.
12.1 Robust Standard Errors¶
Show code, cell 27
# Robust standard errors for the final Random Effects model
re_final_robust <- coeftest(
re_final,
vcov = vcovHC(
re_final,
method = "arellano",
type = "HC1",
cluster = "group"
)
)
re_final_robust
התיקון אינו משנה את אומדי המקדמים, אלא מתקן את אומדני Standard Errors ואת מבחני המובהקות שלהם.
חישוב מטריצת השונות בוצע בשיטת Arellano HC1 עם קיבוץ לפי מדינה (cluster = "group"). כך ההסקה עמידה להטרוסקדסטיות ולתלות אפשרית בין תצפיות החוזרות של אותה מדינה לאורך זמן.
לאחר התיקון, המשתנה male נותר מובהק ברמת 1% (p-value = 0.0082).
לעומת זאת, המשתנה law אינו מובהק לאחר התיקון (p-value = 0.447).
לכן, לאחר התחשבות ב- heteroskedasticity, אין עדות סטטיסטית מספקת לכך שקיומו של חוק Shall-Carry קשור לשיעור הפשיעה האלימה, כאשר יתר המשתנים וההשפעות השנתיות מוחזקים קבועים.
12.2 FGLS Correction Attempt¶
בנוסף לתיקון השגיאות החסינות, ננסה לתקן את heteroskedasticity באמצעות FGLS.
בשלב הראשון נאמוד את שונות השגיאות באמצעות רגרסיה של ריבועי השאריות על המשתנים המסבירים. השונויות החזויות ישמשו ליצירת משקולות עבור מודל GLS.
לאחר מכן נבצע את התהליך באופן איטרטיבי: בכל איטרציה יחושבו שאריות חדשות, יאומד מחדש מודל השונות וייאמד מודל GLS חדש. התהליך ייעצר כאשר השינוי בערכים החזויים יהיה קטן מסף ההתכנסות.
מפרט הממוצע שנשמר הוא:
חשוב להדגיש כי nlme::gls אינו משחזר את מבנה ה- Random Effects של plm. לכן ה- FGLS ישמש כניסיון לתיקון מבנה השגיאה וכבדיקת עמידות, ולא כתחליף למודל RE שנבחר באמצעות Hausman.
Show code, cell 28
# Initial variance estimation for FGLS
fgls_residuals <- as.numeric(resid(re_final))
variance_model <- lm(
I(fgls_residuals^2) ~ law + male + year,
data = guns
)
fitted_variances <- fitted(variance_model)
# Initial FGLS model
fgls_model <- gls(
formula_final,
data = guns,
weights = varFixed(~ fitted_variances)
)
summary(fgls_model)
Show code, cell 29
# Iterate FGLS until convergence
tolerance <- 1e-6
max_iterations <- 10
iteration <- 1
converged <- FALSE
while (iteration <= max_iterations & !converged) {
# Update residuals
fgls_residuals <- residuals(fgls_model)
# Re-estimate the variance model
variance_model <- lm(
I(fgls_residuals^2) ~ law + male + year,
data = guns
)
new_fitted_variances <- fitted(variance_model)
# Re-estimate FGLS with updated variance estimates
new_fgls_model <- gls(
formula_final,
data = guns,
weights = varFixed(~ new_fitted_variances)
)
# Check convergence
if (sum((fitted(new_fgls_model) - fitted(fgls_model))^2) < tolerance) {
converged <- TRUE
}
# Update
fgls_model <- new_fgls_model
fitted_variances <- new_fitted_variances
iteration <- iteration + 1
}
cat("FGLS converged:", converged, "\n")
cat("Iterations:", iteration - 1, "\n")
summary(fgls_model)
Show code, cell 30
# Check the estimated variance function
range(new_fitted_variances)
sum(new_fitted_variances <= 0)
אומדני השונות שהתקבלו ממודל השונות נמצאים בטווח 0.077–0.652, ולא התקבלו אומדני שונות שליליים.
לכן פונקציית השונות ששימשה ליצירת המשקולות ב- FGLS תקינה מבחינה בסיסית, וניתן להמשיך ולבדוק האם האמידה המשוקללת אכן שיפרה את התנהגות השאריות.
12.3 FGLS Residual Diagnostics¶
כדי לבדוק האם תיקון FGLS הצליח לצמצם את ה- heteroskedasticity, נבחן מחדש את השאריות לאחר האמידה המשוקללת.
נשתמש בשאריות המנורמלות של מודל FGLS, ונבצע מחדש את גרף Residuals vs Fitted ואת מבחן Breusch–Pagan.
Show code, cell 31
# FGLS residual diagnostics
res_fgls <- as.numeric(residuals(fgls_model, type = "normalized"))
fit_fgls <- as.numeric(fitted(fgls_model))
options(repr.plot.width = 12, repr.plot.height = 7)
plot(
fit_fgls,
res_fgls,
main = "Residuals vs Fitted - FGLS Model",
xlab = "Fitted Values",
ylab = "Normalized Residuals",
pch = 20
)
abline(h = 0, col = "red")
Show code, cell 32
# Breusch-Pagan test after FGLS
hetero_test_fgls <- bptest(
res_fgls ~ fit_fgls + I(fit_fgls^2),
studentize = FALSE
)
hetero_test_fgls
מודל FGLS התכנס לאחר 4 איטרציות, אך בדיקת השאריות מראה כי תיקון ההטרוסקדסטיות לא הצליח.
במבחן Breusch–Pagan לאחר FGLS התקבל p-value = 0.0047.
מכיוון שערך זה קטן מ-0.05, אנו עדיין דוחים את השערת האפס של שונות קבועה בשגיאות.
לכן, למרות שהאלגוריתם התכנס והמשקולות התבססו על אומדני שונות חיוביים, מודל FGLS לא הצליח להסיר את ה- heteroskedasticity מהשאריות.
בהתאם לכך, תוצאות מודל Random Effects עם robust standard errors יישארו הבסיס המרכזי להסקה הסטטיסטית, בעוד שה- FGLS יוצג כניסיון נוסף לתיקון מבנה השגיאה שלא פתר את הבעיה במלואה.
12.4 Normality after FGLS¶
נבחן גם האם אמידת FGLS שיפרה את נורמליות השאריות באמצעות Q-Q Plot, היסטוגרמה ומבחן Jarque–Bera.
Show code, cell 33
# Normality diagnostics after FGLS
options(repr.plot.width = 14, repr.plot.height = 6)
par(mfrow = c(1, 2))
qqnorm(
res_fgls,
main = "Q-Q Plot - FGLS Residuals"
)
qqline(res_fgls, col = "red")
hist(
res_fgls,
breaks = 30,
freq = FALSE,
main = "Histogram - FGLS Residuals",
xlab = "Normalized Residuals"
)
lines(density(res_fgls), lwd = 2)
curve(
dnorm(
x,
mean = mean(res_fgls),
sd = sd(res_fgls)
),
col = "red",
lwd = 2,
add = TRUE
)
par(mfrow = c(1, 1))
Show code, cell 34
# Jarque-Bera test after FGLS
jb_test_fgls <- jarque.bera.test(res_fgls)
jb_test_fgls
גם לאחר אמידת FGLS מבחן Jarque–Bera דוחה את השערת הנורמליות, עם p-value = 1.836×10⁻¹⁰.
גם מבחינה גרפית, גרף Q-Q מצביע על התאמה סבירה יחסית במרכז ההתפלגות, אך על סטיות ברורות בזנבות.
בנוסף, מבחן Breusch–Pagan שבוצע לאחר FGLS התקבל עם p-value = 0.0047, ולכן גם בעיית ה- heteroskedasticity לא נפתרה.
לפיכך, אף שמודל FGLS התכנס לאחר ארבע איטרציות והשתמש באומדני שונות חיוביים, ניסיון התיקון לא הצליח להביא לקיום מלא של הנחות השגיאה.
בהתאם לכך, מודל Random Effects נשאר המודל המרכזי של הניתוח, וההסקה הסטטיסטית ממנו תתבסס על robust standard errors כדי להתמודד עם ההטרוסקדסטיות.
13. Final Results and Conclusions¶
בשלב הראשון נאמדו מודלי Fixed Effects ו-Random Effects על אותו מפרט.
לאחר טיפול ב- multicollinearity בין log_population ו-log_density, בוצע מבחן Hausman. במבחן התקבל p-value = 0.0892, ולכן ברמת מובהקות של 5% לא דחינו את השערת האפס ובחרנו במודל Random Effects.
לאחר בחירת המודל בוצע סינון הדרגתי של המשתנים המסבירים. המודל הסופי שהתקבל הוא:
במודל Random Effects הסופי התקבל עבור law מקדם של -0.0321.
אומדן נקודתי זה מתאים בקירוב לירידה של כ-3.16% בשיעור הפשיעה האלימה כאשר חוק Shall-Carry בתוקף.
עם זאת, לאחר שנמצאה heteroskedasticity וחושבו robust standard errors, התקבל עבור law ערך p-value = 0.447.
לכן אין עדות סטטיסטית מספקת לכך שקיומו של חוק Shall-Carry קשור לשינוי בשיעור הפשיעה האלימה. כלומר, אין לפרש את הירידה הנקודתית של כ-3.16% כהשפעה מובהקת.
עבור male התקבל מקדם של 0.0802, והמשתנה נשאר מובהק גם לאחר תיקון השגיאות (p-value = 0.0082).
מכיוון שהמשתנה המוסבר נמצא בסקאלה לוגריתמית, עלייה של נקודת אחוז אחת ב- male קשורה בקירוב לעלייה של כ-8% בשיעור הפשיעה האלימה, בהינתן יתר המשתנים במודל.
גם משתני השנה מצביעים על שינויים משמעותיים בשיעור הפשיעה האלימה לאורך תקופת המדגם ביחס לשנת הבסיס 1977.
למודל התקבל R² = 0.399