Counterfactual Estimates of Car-free Friedrichstrasse Using Fect
Yogada Joshi
2025-04-18
# libraries
library(fect)## Warning: package 'fect' was built under R version 4.3.3## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(panelView)## ## See bit.ly/panelview4r for more info.
## ## Report bugs -> yiqingxu@stanford.edu.library(abind)## Warning: package 'abind' was built under R version 4.3.3library(doParallel)## Warning: package 'doParallel' was built under R version 4.3.3## Loading required package: foreach## Loading required package: iterators## Loading required package: parallellibrary(doRNG) ## Warning: package 'doRNG' was built under R version 4.3.3## Loading required package: rngtools## Warning: package 'rngtools' was built under R version 4.3.3library(fixest)## Warning: package 'fixest' was built under R version 4.3.3library(foreach)
library(future)
library(GGally) ## Warning: package 'GGally' was built under R version 4.3.3## Loading required package: ggplot2library(ggplot2)
library(grid)
library(gridExtra)
library(Rcpp)
library(MASS)
# for code language:
library(dplyr)## Warning: package 'dplyr' was built under R version 4.3.2##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# for data prep:
library(lubridate)## Warning: package 'lubridate' was built under R version 4.3.3##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(gsynth)## Warning: package 'gsynth' was built under R version 4.3.3## Registered S3 method overwritten by 'gsynth':
## method from
## print.interFE fect## ## Syntax has been updated since v.1.2.0.## ## Comments and suggestions -> yiqingxu@stanford.edu.##
## Attaching package: 'gsynth'## The following object is masked from 'package:fect':
##
## interFElibrary(plm)##
## Attaching package: 'plm'## The following objects are masked from 'package:dplyr':
##
## between, lag, lead# To plot from plm and fixest
library(modelsummary)## Warning: package 'modelsummary' was built under R version 4.3.2library(ggplot2)
library(broom) ## Warning: package 'broom' was built under R version 4.3.3## Registered S3 method overwritten by 'broom':
## method from
## nobs.felm lfeData preparation to suit Fect and PanelView
# load normalized panel dataframe
N_accidents <- read.csv("streetwise_accident_panel_normalized_2018_2023.csv")
colnames(N_accidents)## [1] "Street_Group" "Year"
## [3] "Month" "YearMonth"
## [5] "D1" "D2"
## [7] "Count_Total" "Count_Rad"
## [9] "Count_Fuss" "Count_PKW"
## [11] "Count_Krad" "Count_Gkfz"
## [13] "Count_Sonstige" "Count_PKW_vs_Rad"
## [15] "Count_PKW_vs_Fuss" "Count_PKW_vs_VRU"
## [17] "Street_Length_km" "Rate_Total_per_km"
## [19] "Rate_Rad_per_km" "Rate_Fuss_per_km"
## [21] "Rate_PKW_per_km" "Rate_Krad_per_km"
## [23] "Rate_Gkfz_per_km" "Rate_Sonstige_per_km"
## [25] "Rate_PKW_vs_Rad_per_km" "Rate_PKW_vs_Fuss_per_km"
## [27] "Rate_PKW_vs_VRU_per_km"# data prep for panelview and fect:
# creates single D variable, pooled D1 and D2
N_accidents <- N_accidents %>%
mutate(D = ifelse(D1 == 1 | D2 == 1, 1, 0))
# time variable in count of 72 periods:
class(N_accidents$YearMonth)## [1] "character"N_accidents <- N_accidents %>%
# 1. Convert character string to Date objects
mutate(Date = ymd(YearMonth)) %>%
mutate(TimeIndex = 12 * (year(Date) - 2018) + month(Date)) %>%
arrange(Street_Group, TimeIndex)Visualized Treatment Periods and Units
panelview(Rate_Total_per_km ~ D, data = N_accidents, index = c("Street_Group","TimeIndex"),
axis.lab = "time", xlab = "Time", ylab = "Unit",
gridOff = TRUE, by.timing = TRUE,
background = "white", main = "Flaniermeile Friedrichstrasse and Controls: Treatment Status")
Visualising the outcome variable Normalized Total Accidents (Rate_Total_per_km)
panelview(Rate_Total_per_km ~ D, data = N_accidents, index = c("Street_Group","TimeIndex"),
axis.lab = "time", xlab = "Time", ylab = "Unit",
theme.bw = TRUE, type = "outcome",
main = "Normalized Total Accidents (Rate_Total_per_km)")## Treatment has reversals.
Model-1: Effect on total accidents (normalized by length (km))
out.mc <- fect(formula = Rate_Total_per_km ~ D, # Y ~ D
data = N_accidents,
index = c("Street_Group", "TimeIndex"), # Specify unit and time IDs
method = "mc", # Matrix Completion method
se = TRUE, # Calculate standard errors (using bootstrap/jackknife)
vartype = "jackknife", # Or "jackknife", recommended if few treated units [cite: 175]
nboots = 200, # Number of bootstrap iterations (adjust as needed, 200 is often ok, more is slower)
parallel = TRUE,
CV = TRUE, # Use cross-validation to choose lambda (tuning parameter for MC) [cite: 172]
nlambda = 10, # Optionally specify how many lambdas CV should try [cite: 1]
criterion = "moment",
force = "two-way") # Assuming two-way fixed effects structure as baseline## Parallel computing ...## Cross-validating ...## Criterion: Moment Conditions## Matrix completion method...##
## lambda.norm = 1.00000; MSPE = 1.13156; MSPTATT = 2.48177; MSE = 3.55987##
## lambda.norm = 0.42170; MSPE = 1.09900; MSPTATT = 0.45365; MSE = 0.66141##
## lambda.norm = 0.17783; MSPE = 1.11703; MSPTATT = 0.08058; MSE = 0.11774## *##
## lambda.norm = 0.07499; MSPE = 1.11834; MSPTATT = 0.01434; MSE = 0.02093##
## lambda.norm = 0.03162; MSPE = 1.11889; MSPTATT = 0.00255; MSE = 0.00371##
## lambda.norm = 0.01334; MSPE = 1.12513; MSPTATT = 0.00045; MSE = 0.00066##
##
## lambda.norm* = 0.177827941003892## ## Jackknife estimates ...## 10 runs## Cannot use full pre-treatment periods in the F test. The first period is removed.## ## The estimated covariance matrix is irreversible.## # Print summary of the results (ATT etc.)
summary(out.mc)## Length Class Mode
## Y.dat 792 -none- numeric
## D.dat 792 -none- numeric
## I.dat 792 -none- numeric
## Y 1 -none- character
## D 1 -none- character
## X 0 -none- NULL
## T.on 792 -none- numeric
## G 0 -none- NULL
## balance.period 0 -none- NULL
## hasRevs 1 -none- numeric
## T.off 792 -none- numeric
## T.on.balance 792 -none- logical
## index 2 -none- character
## id 11 -none- character
## rawtime 72 -none- numeric
## binary 1 -none- logical
## loo 1 -none- logical
## proportion 1 -none- numeric
## pre.periods 32 -none- numeric
## tost.threshold 1 -none- numeric
## placeboTest 1 -none- logical
## placebo.period 0 -none- NULL
## carryoverTest 1 -none- logical
## carryover.period 0 -none- NULL
## unit.type 11 -none- numeric
## obs.missing 792 -none- numeric
## obs.missing.balance 1 -none- logical
## sigma2 1 -none- numeric
## sigma2.fect 1 -none- numeric
## T.on 792 -none- numeric
## Y.ct 792 -none- numeric
## eff 792 -none- numeric
## att.avg 1 -none- numeric
## att.avg.unit 1 -none- numeric
## force 1 -none- numeric
## T 1 -none- numeric
## N 1 -none- numeric
## Ntr 1 -none- numeric
## Nco 1 -none- numeric
## p 1 -none- numeric
## D 792 -none- numeric
## Y 792 -none- numeric
## X 0 -none- numeric
## I 792 -none- numeric
## II 792 -none- numeric
## tr 1 -none- numeric
## co 10 -none- numeric
## est 9 -none- list
## method 1 -none- character
## mu 1 -none- numeric
## beta 1 -none- logical
## validX 1 -none- numeric
## validF 1 -none- numeric
## niter 1 -none- numeric
## time 46 -none- numeric
## att 46 -none- numeric
## count 46 -none- numeric
## eff.calendar 144 -none- numeric
## N.calendar 72 -none- numeric
## eff.calendar.fit 144 -none- numeric
## calendar.enp 1 -none- numeric
## eff.pre 141 -none- numeric
## eff.pre.equiv 141 -none- numeric
## pre.sd 96 -none- numeric
## rmse 1 -none- numeric
## rmCV 10 -none- list
## estCV 10 -none- list
## res 792 -none- numeric
## time.off 29 -none- numeric
## att.off 29 -none- numeric
## count.off 29 -none- numeric
## eff.off 63 -none- numeric
## eff.off.equiv 63 -none- numeric
## off.sd 45 -none- numeric
## alpha 11 -none- numeric
## xi 72 -none- numeric
## alpha.tr 1 -none- numeric
## alpha.co 10 -none- numeric
## lambda.cv 1 -none- numeric
## lambda.seq 10 -none- numeric
## lambda.norm 1 -none- numeric
## eigen.all 11 -none- numeric
## CV.out.mc 100 -none- numeric
## est.avg 5 -none- numeric
## att.avg.boot 10 -none- numeric
## est.avg.unit 5 -none- numeric
## att.avg.unit.boot 10 -none- numeric
## est.eff.calendar 432 -none- numeric
## est.eff.calendar.fit 432 -none- numeric
## est.att 276 -none- numeric
## att.bound 92 -none- numeric
## att.boot 460 -none- numeric
## att.count.boot 460 -none- numeric
## est.att.off 174 -none- numeric
## att.off.boot 290 -none- numeric
## att.off.bound 58 -none- numeric
## att.off.count.boot 290 -none- numeric
## test.out 9 -none- list
## call 13 -none- call
## formula 3 formula callplot(out.mc, main = "Estimated ATT (MC)")
::: Warning of unreliable estimates: The estimated covariance matrix is
irreversible. Avg ATT is significant at 5% level but cannot trust these
estimates. The plot shows that many pre-treatment periods do not hover
around zero but are closer to zero than post treatment.
The post-treatment estimates do not follow any particular
pattern.
:::
Gap plot of dynamic treatment effects
# We can use the plot function to visualize the estimation results. By default, the plot function produces a “gap” plot – as if we type plot(out.fect, type = "gap") — which visualizes the estimated period-wise ATT (dynamic treatment effects).
plot(out.mc, main = "Estimated ATT (MC)", ylab = "Effect of D on Y",
cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
#The plot() function can visualize the estimated period-wise ATTs as well as their uncertainty estimates. stats = "F.p" shows the p-value for the F test of no-pretrend.
plot(out.mc, main = "Estimated ATT (MC)", ylab = "Effect of D on Y",
cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p")
::: Cannot use full pre-treatment periods in the F test. The first period is removed. Showed as a warning during the model estimation. Hence, I cannot compute this statistic. :::
# Result summary : est.avg and est.avg.unit show the ATT averaged over all periods – the former weights each treated observation equally while the latter weights each treated unit equally. est.beta reports the coefficients of the time-varying covariates. est.att reports the average treatment effect on the treated (ATT) by period. Treatment effect estimates from each bootstrap run is stored in eff.boot, an array whose dimension = (#time periods * #treated * #bootstrap runs).
plot(out.mc, main = "Estimated ATT (MC)", ylab = "Effect of D on Y",
cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p")
print(out.mc)## Call:
## fect.formula(formula = Rate_Total_per_km ~ D, data = N_accidents,
## index = c("Street_Group", "TimeIndex"), force = "two-way",
## nlambda = 10, CV = TRUE, criterion = "moment", method = "mc",
## se = TRUE, vartype = "jackknife", nboots = 200, parallel = TRUE)
##
## ATT:
## ATT S.E. CI.lower CI.upper p.value
## Tr obs equally weighted -0.5305 0.0538 -0.6359 -0.425 6.218e-23
## Tr units equally weighted -0.5305 0.0538 -0.6359 -0.425 6.218e-23out.mc$est.att## ATT S.E. CI.lower CI.upper p.value count
## -31 0.13290292 0.02404192 0.085781629 0.18002422 3.239635e-08 1
## -30 0.03834833 0.04362907 -0.047163082 0.12385975 3.794215e-01 1
## -29 -0.28642098 0.04746583 -0.379452300 -0.19338967 1.596963e-09 1
## -28 0.06463350 0.02833430 0.009099282 0.12016771 2.254229e-02 1
## -27 0.64590069 0.03256267 0.582079039 0.70972235 0.000000e+00 1
## -26 0.37402019 0.03127708 0.312718243 0.43532214 0.000000e+00 1
## -25 0.05007957 0.04249763 -0.033214252 0.13337340 2.386337e-01 1
## -24 0.02613011 0.04712243 -0.066228148 0.11848837 5.792262e-01 1
## -23 -0.01297782 0.02956922 -0.070932434 0.04497679 6.607366e-01 1
## -22 0.43347844 0.05041409 0.334668643 0.53228823 0.000000e+00 1
## -21 0.76234634 0.03327085 0.697136667 0.82755602 0.000000e+00 1
## -20 0.40659027 0.02901385 0.349724175 0.46345637 0.000000e+00 1
## -19 -0.27182255 0.03602618 -0.342432564 -0.20121254 4.518028e-14 1
## -18 0.07282018 0.02243490 0.028848585 0.11679177 1.171030e-03 1
## -17 0.75142152 0.03343640 0.685887386 0.81695565 0.000000e+00 1
## -16 -0.30567345 0.03674180 -0.377686058 -0.23366084 8.833445e-17 1
## -15 -0.30284499 0.03679533 -0.374962508 -0.23072747 1.863895e-16 1
## -14 -0.28412319 0.02919665 -0.341347572 -0.22689882 2.216017e-22 2
## -13 0.38357585 0.02403295 0.336472129 0.43067956 0.000000e+00 2
## -12 0.03150712 0.02662451 -0.020675967 0.08369020 2.366556e-01 2
## -11 -0.10491511 0.01858780 -0.141346525 -0.06848370 1.658541e-08 2
## -10 -0.11203303 0.02177481 -0.154710871 -0.06935518 2.674141e-07 2
## -9 0.04666754 0.03315376 -0.018312633 0.11164771 1.592468e-01 2
## -8 -0.29079100 0.01612405 -0.322393545 -0.25918845 1.042131e-72 2
## -7 -0.28495436 0.01485938 -0.314078218 -0.25583051 5.791124e-82 2
## -6 -0.29524891 0.03006193 -0.354169210 -0.23632861 9.110974e-23 2
## -5 -0.25975179 0.02803713 -0.314703543 -0.20480003 1.958731e-20 2
## -4 -0.11275999 0.03054627 -0.172629588 -0.05289040 2.229805e-04 2
## -3 -0.09273214 0.02888108 -0.149338017 -0.03612626 1.323537e-03 2
## -2 -0.12106445 0.02251499 -0.165193007 -0.07693589 7.571159e-08 2
## -1 -0.25672869 0.02229158 -0.300419385 -0.21303800 1.085101e-30 2
## 0 0.04653192 0.01623495 0.014711993 0.07835185 4.154889e-03 2
## 1 -0.91904269 0.20052686 -1.312068120 -0.52601726 4.580450e-06 2
## 2 -0.40421262 0.14524636 -0.688890263 -0.11953498 5.386796e-03 2
## 3 -1.38722222 0.08647169 -1.556703612 -1.21774083 6.450210e-58 2
## 4 0.32347930 0.10131669 0.124902230 0.52205636 1.409228e-03 2
## 5 1.57898155 0.09469659 1.393379635 1.76458347 0.000000e+00 2
## 6 -0.95589906 0.11153924 -1.174511958 -0.73728615 1.034228e-17 1
## 7 -0.87794304 0.12404083 -1.121058600 -0.63482747 1.464027e-12 1
## 8 -1.05069728 0.13921075 -1.323545335 -0.77784923 4.435907e-14 1
## 9 0.69374435 0.17159421 0.357425880 1.03006282 5.278610e-05 1
## 10 0.57494121 0.18294303 0.216379455 0.93350296 1.673781e-03 1
## 11 -1.64995609 0.18605685 -2.014620825 -1.28529136 7.445729e-19 1
## 12 -1.65126164 0.18815542 -2.020039487 -1.28248379 1.693146e-18 1
## 13 -1.93481025 0.13540613 -2.200201392 -1.66941910 2.564804e-46 1
## 14 -1.61094983 0.15749712 -1.919638514 -1.30226115 1.478824e-24 1out.mc$est.avg## ATT.avg S.E. CI.lower CI.upper p.value
## [1,] -0.5304666 0.05380099 -0.6359146 -0.4250186 6.217792e-23out.mc$est.beta## NULLModel:2 Effects on Log-transformed Normalised Car vs VRU Accidents using Matrix Completion Estimator
NL_accidents <- read.csv("streetwise_accident_panel_log_transformed_2018_2023.csv")#data preparation
NL_accidents <- NL_accidents %>%
mutate(D = ifelse(D1 == 1 | D2 == 1, 1, 0))
class(NL_accidents$YearMonth)## [1] "character"NL_accidents <- NL_accidents %>%
mutate(Date = ymd(YearMonth)) %>%
mutate(TimeIndex = 12 * (year(Date) - 2018) + month(Date)) %>%
arrange(Street_Group, TimeIndex)
colnames(NL_accidents)## [1] "Street_Group" "Year"
## [3] "Month" "YearMonth"
## [5] "D1" "D2"
## [7] "Count_Total" "Count_Rad"
## [9] "Count_Fuss" "Count_PKW"
## [11] "Count_Krad" "Count_Gkfz"
## [13] "Count_Sonstige" "Count_PKW_vs_Rad"
## [15] "Count_PKW_vs_Fuss" "Count_PKW_vs_VRU"
## [17] "Street_Length_km" "Rate_Total_per_km"
## [19] "Rate_Rad_per_km" "Rate_Fuss_per_km"
## [21] "Rate_PKW_per_km" "Rate_Krad_per_km"
## [23] "Rate_Gkfz_per_km" "Rate_Sonstige_per_km"
## [25] "Rate_PKW_vs_Rad_per_km" "Rate_PKW_vs_Fuss_per_km"
## [27] "Rate_PKW_vs_VRU_per_km" "log_Count_Total"
## [29] "log_Count_Rad" "log_Count_Fuss"
## [31] "log_Count_PKW" "log_Count_Krad"
## [33] "log_Count_Gkfz" "log_Count_Sonstige"
## [35] "log_Count_PKW_vs_Rad" "log_Count_PKW_vs_Fuss"
## [37] "log_Count_PKW_vs_VRU" "log_Rate_Total_per_km"
## [39] "log_Rate_Rad_per_km" "log_Rate_Fuss_per_km"
## [41] "log_Rate_PKW_per_km" "log_Rate_Krad_per_km"
## [43] "log_Rate_Gkfz_per_km" "log_Rate_Sonstige_per_km"
## [45] "log_Rate_PKW_vs_Rad_per_km" "log_Rate_PKW_vs_Fuss_per_km"
## [47] "log_Rate_PKW_vs_VRU_per_km" "D"
## [49] "Date" "TimeIndex"Visualising the outcome variable
panelview(log_Count_PKW_vs_VRU ~ D, data = NL_accidents, index = c("Street_Group","TimeIndex"),
axis.lab = "time", xlab = "Time", ylab = "Unit",
theme.bw = TRUE, type = "outcome",
main = "Log-transformed Normalised Car vs VRU Accidents (log_Count_PKW_vs_VRU)")## Treatment has reversals.
out.lnmc <- fect(formula = log_Count_PKW_vs_VRU ~ D, # Y ~ D
data = NL_accidents,
index = c("Street_Group", "TimeIndex"), # Specify unit and time IDs
method = "mc", # Matrix Completion method
se = TRUE, # Calculate standard errors (using bootstrap/jackknife)
vartype = "jackknife", # Or "jackknife", recommended if few treated units [cite: 175]
nboots = 200, # Number of bootstrap iterations (adjust as needed, 200 is often ok, more is slower)
parallel = TRUE,
CV = TRUE, # Use cross-validation to choose lambda (tuning parameter for MC) [cite: 172]
nlambda = 10, # Optionally specify how many lambdas CV should try [cite: 1]
criterion = "moment",
force = "two-way")## Parallel computing ...## Cross-validating ...## Criterion: Moment Conditions## Matrix completion method...##
## lambda.norm = 1.00000; MSPE = 0.24654; MSPTATT = 0.08113; MSE = 0.11131## *##
## lambda.norm = 0.42170; MSPE = 0.25623; MSPTATT = 0.05070; MSE = 0.07534##
## lambda.norm = 0.17783; MSPE = 0.25594; MSPTATT = 0.00912; MSE = 0.01354##
## lambda.norm = 0.07499; MSPE = 0.25216; MSPTATT = 0.00163; MSE = 0.00240##
##
## lambda.norm* = 1## ## Jackknife estimates ...## 10 runs## Cannot use full pre-treatment periods in the F test. The first period is removed.## ## The estimated covariance matrix is irreversible.## # Print summary of the results (ATT etc.)
summary(out.lnmc)## Length Class Mode
## Y.dat 792 -none- numeric
## D.dat 792 -none- numeric
## I.dat 792 -none- numeric
## Y 1 -none- character
## D 1 -none- character
## X 0 -none- NULL
## T.on 792 -none- numeric
## G 0 -none- NULL
## balance.period 0 -none- NULL
## hasRevs 1 -none- numeric
## T.off 792 -none- numeric
## T.on.balance 792 -none- logical
## index 2 -none- character
## id 11 -none- character
## rawtime 72 -none- numeric
## binary 1 -none- logical
## loo 1 -none- logical
## proportion 1 -none- numeric
## pre.periods 32 -none- numeric
## tost.threshold 1 -none- numeric
## placeboTest 1 -none- logical
## placebo.period 0 -none- NULL
## carryoverTest 1 -none- logical
## carryover.period 0 -none- NULL
## unit.type 11 -none- numeric
## obs.missing 792 -none- numeric
## obs.missing.balance 1 -none- logical
## sigma2 1 -none- numeric
## sigma2.fect 1 -none- numeric
## T.on 792 -none- numeric
## Y.ct 792 -none- numeric
## eff 792 -none- numeric
## att.avg 1 -none- numeric
## att.avg.unit 1 -none- numeric
## force 1 -none- numeric
## T 1 -none- numeric
## N 1 -none- numeric
## Ntr 1 -none- numeric
## Nco 1 -none- numeric
## p 1 -none- numeric
## D 792 -none- numeric
## Y 792 -none- numeric
## X 0 -none- numeric
## I 792 -none- numeric
## II 792 -none- numeric
## tr 1 -none- numeric
## co 10 -none- numeric
## est 9 -none- list
## method 1 -none- character
## mu 1 -none- numeric
## beta 1 -none- logical
## validX 1 -none- numeric
## validF 1 -none- numeric
## niter 1 -none- numeric
## time 46 -none- numeric
## att 46 -none- numeric
## count 46 -none- numeric
## eff.calendar 144 -none- numeric
## N.calendar 72 -none- numeric
## eff.calendar.fit 144 -none- numeric
## calendar.enp 1 -none- numeric
## eff.pre 141 -none- numeric
## eff.pre.equiv 141 -none- numeric
## pre.sd 96 -none- numeric
## rmse 1 -none- numeric
## rmCV 10 -none- list
## estCV 10 -none- list
## res 792 -none- numeric
## time.off 29 -none- numeric
## att.off 29 -none- numeric
## count.off 29 -none- numeric
## eff.off 63 -none- numeric
## eff.off.equiv 63 -none- numeric
## off.sd 45 -none- numeric
## alpha 11 -none- numeric
## xi 72 -none- numeric
## alpha.tr 1 -none- numeric
## alpha.co 10 -none- numeric
## lambda.cv 1 -none- numeric
## lambda.seq 10 -none- numeric
## lambda.norm 1 -none- numeric
## eigen.all 11 -none- numeric
## CV.out.mc 100 -none- numeric
## est.avg 5 -none- numeric
## att.avg.boot 10 -none- numeric
## est.avg.unit 5 -none- numeric
## att.avg.unit.boot 10 -none- numeric
## est.eff.calendar 432 -none- numeric
## est.eff.calendar.fit 432 -none- numeric
## est.att 276 -none- numeric
## att.bound 92 -none- numeric
## att.boot 460 -none- numeric
## att.count.boot 460 -none- numeric
## est.att.off 174 -none- numeric
## att.off.boot 290 -none- numeric
## att.off.bound 58 -none- numeric
## att.off.count.boot 290 -none- numeric
## test.out 9 -none- list
## call 13 -none- call
## formula 3 formula callplot(out.lnmc, main = "Estimated ATT (MC)")
::: This model also shows a gap plot where the counterfactual estimates
are not matching the real estimates on pre-treatment periods.
Due to the model warnings and poor pre-trending outcome matches + noise in outcome plots, I did not check the placebo tests. :::
Model-3: Using 7 control streets to reduce pre-trend noise
Outcome variable is the same as Model-2 (Log-transformed Normalised Car vs VRU Accidents )
Matrix completion estimator
::: Removing Boxhagnerstr., Wilmersdorfer Straße, Kastanienallee and Unter den Linden and retaining more commercial through-ways similar to Friedrichstrasse. :::
# data preparation
streets_to_exclude <- c("Boxhagenerstraße",
"Kastanienallee",
"Unter den Linden",
"Wilmersdorfer Straße")
# Create the new dataframe 'FNL_Accidents' by filtering out the specified streets
FNL_Accidents <- NL_accidents %>%
filter(!Street_Group %in% streets_to_exclude)panelview(log_Count_PKW_vs_VRU ~ D, data = FNL_Accidents, index = c("Street_Group","TimeIndex"),
axis.lab = "time", xlab = "Time", ylab = "Unit",
theme.bw = TRUE, type = "outcome",
main = "Log-transformed Normalised Car vs VRU Accidents (log_Count_PKW_vs_VRU)")## Treatment has reversals.
out.flnmc <- fect(formula = log_Count_PKW_vs_VRU ~ D, # Y ~ D
data = FNL_Accidents,
index = c("Street_Group", "TimeIndex"), # Specify unit and time IDs
method = "mc", # Matrix Completion method
se = TRUE, # Calculate standard errors (using bootstrap/jackknife)
vartype = "jackknife", # Or "jackknife", recommended if few treated units [cite: 175]
nboots = 200, # Number of bootstrap iterations (adjust as needed, 200 is often ok, more is slower)
parallel = TRUE,
CV = TRUE, # Use cross-validation to choose lambda (tuning parameter for MC) [cite: 172]
nlambda = 10, # Optionally specify how many lambdas CV should try [cite: 1]
criterion = "moment",
force = "two-way") # Assuming two-way fixed effects structure as baseline## Parallel computing ...## Cross-validating ...## Criterion: Moment Conditions## Matrix completion method...##
## lambda.norm = 1.00000; MSPE = 0.31255; MSPTATT = 0.06209; MSE = 0.09518## *##
## lambda.norm = 0.42170; MSPE = 0.32505; MSPTATT = 0.04430; MSE = 0.06964##
## lambda.norm = 0.17783; MSPE = 0.32354; MSPTATT = 0.00788; MSE = 0.01237##
## lambda.norm = 0.07499; MSPE = 0.31938; MSPTATT = 0.00140; MSE = 0.00219##
##
## lambda.norm* = 1## ## Jackknife estimates ...## 6 runs## Cannot use full pre-treatment periods in the F test. The first period is removed.## ## The estimated covariance matrix is irreversible.## # Print summary of the results (ATT etc.)
summary(out.flnmc)## Length Class Mode
## Y.dat 504 -none- numeric
## D.dat 504 -none- numeric
## I.dat 504 -none- numeric
## Y 1 -none- character
## D 1 -none- character
## X 0 -none- NULL
## T.on 504 -none- numeric
## G 0 -none- NULL
## balance.period 0 -none- NULL
## hasRevs 1 -none- numeric
## T.off 504 -none- numeric
## T.on.balance 504 -none- logical
## index 2 -none- character
## id 7 -none- character
## rawtime 72 -none- numeric
## binary 1 -none- logical
## loo 1 -none- logical
## proportion 1 -none- numeric
## pre.periods 32 -none- numeric
## tost.threshold 1 -none- numeric
## placeboTest 1 -none- logical
## placebo.period 0 -none- NULL
## carryoverTest 1 -none- logical
## carryover.period 0 -none- NULL
## unit.type 7 -none- numeric
## obs.missing 504 -none- numeric
## obs.missing.balance 1 -none- logical
## sigma2 1 -none- numeric
## sigma2.fect 1 -none- numeric
## T.on 504 -none- numeric
## Y.ct 504 -none- numeric
## eff 504 -none- numeric
## att.avg 1 -none- numeric
## att.avg.unit 1 -none- numeric
## force 1 -none- numeric
## T 1 -none- numeric
## N 1 -none- numeric
## Ntr 1 -none- numeric
## Nco 1 -none- numeric
## p 1 -none- numeric
## D 504 -none- numeric
## Y 504 -none- numeric
## X 0 -none- numeric
## I 504 -none- numeric
## II 504 -none- numeric
## tr 1 -none- numeric
## co 6 -none- numeric
## est 9 -none- list
## method 1 -none- character
## mu 1 -none- numeric
## beta 1 -none- logical
## validX 1 -none- numeric
## validF 1 -none- numeric
## niter 1 -none- numeric
## time 46 -none- numeric
## att 46 -none- numeric
## count 46 -none- numeric
## eff.calendar 144 -none- numeric
## N.calendar 72 -none- numeric
## eff.calendar.fit 144 -none- numeric
## calendar.enp 1 -none- numeric
## eff.pre 141 -none- numeric
## eff.pre.equiv 141 -none- numeric
## pre.sd 96 -none- numeric
## rmse 1 -none- numeric
## rmCV 10 -none- list
## estCV 10 -none- list
## res 504 -none- numeric
## time.off 29 -none- numeric
## att.off 29 -none- numeric
## count.off 29 -none- numeric
## eff.off 63 -none- numeric
## eff.off.equiv 63 -none- numeric
## off.sd 45 -none- numeric
## alpha 7 -none- numeric
## xi 72 -none- numeric
## alpha.tr 1 -none- numeric
## alpha.co 6 -none- numeric
## lambda.cv 1 -none- numeric
## lambda.seq 10 -none- numeric
## lambda.norm 1 -none- numeric
## eigen.all 7 -none- numeric
## CV.out.mc 100 -none- numeric
## est.avg 5 -none- numeric
## att.avg.boot 6 -none- numeric
## est.avg.unit 5 -none- numeric
## att.avg.unit.boot 6 -none- numeric
## est.eff.calendar 432 -none- numeric
## est.eff.calendar.fit 432 -none- numeric
## est.att 276 -none- numeric
## att.bound 92 -none- numeric
## att.boot 276 -none- numeric
## att.count.boot 276 -none- numeric
## est.att.off 174 -none- numeric
## att.off.boot 174 -none- numeric
## att.off.bound 58 -none- numeric
## att.off.count.boot 174 -none- numeric
## test.out 9 -none- list
## call 13 -none- call
## formula 3 formula callplot(out.flnmc, main = "Estimated ATT (MC)")
::: Noisier estimates :::
plot(out.flnmc, main = "Estimated ATT (MC)", ylab = "Effect of D on Y",
cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p")
print(out.mc)## Call:
## fect.formula(formula = Rate_Total_per_km ~ D, data = N_accidents,
## index = c("Street_Group", "TimeIndex"), force = "two-way",
## nlambda = 10, CV = TRUE, criterion = "moment", method = "mc",
## se = TRUE, vartype = "jackknife", nboots = 200, parallel = TRUE)
##
## ATT:
## ATT S.E. CI.lower CI.upper p.value
## Tr obs equally weighted -0.5305 0.0538 -0.6359 -0.425 6.218e-23
## Tr units equally weighted -0.5305 0.0538 -0.6359 -0.425 6.218e-23out.mc$est.att## ATT S.E. CI.lower CI.upper p.value count
## -31 0.13290292 0.02404192 0.085781629 0.18002422 3.239635e-08 1
## -30 0.03834833 0.04362907 -0.047163082 0.12385975 3.794215e-01 1
## -29 -0.28642098 0.04746583 -0.379452300 -0.19338967 1.596963e-09 1
## -28 0.06463350 0.02833430 0.009099282 0.12016771 2.254229e-02 1
## -27 0.64590069 0.03256267 0.582079039 0.70972235 0.000000e+00 1
## -26 0.37402019 0.03127708 0.312718243 0.43532214 0.000000e+00 1
## -25 0.05007957 0.04249763 -0.033214252 0.13337340 2.386337e-01 1
## -24 0.02613011 0.04712243 -0.066228148 0.11848837 5.792262e-01 1
## -23 -0.01297782 0.02956922 -0.070932434 0.04497679 6.607366e-01 1
## -22 0.43347844 0.05041409 0.334668643 0.53228823 0.000000e+00 1
## -21 0.76234634 0.03327085 0.697136667 0.82755602 0.000000e+00 1
## -20 0.40659027 0.02901385 0.349724175 0.46345637 0.000000e+00 1
## -19 -0.27182255 0.03602618 -0.342432564 -0.20121254 4.518028e-14 1
## -18 0.07282018 0.02243490 0.028848585 0.11679177 1.171030e-03 1
## -17 0.75142152 0.03343640 0.685887386 0.81695565 0.000000e+00 1
## -16 -0.30567345 0.03674180 -0.377686058 -0.23366084 8.833445e-17 1
## -15 -0.30284499 0.03679533 -0.374962508 -0.23072747 1.863895e-16 1
## -14 -0.28412319 0.02919665 -0.341347572 -0.22689882 2.216017e-22 2
## -13 0.38357585 0.02403295 0.336472129 0.43067956 0.000000e+00 2
## -12 0.03150712 0.02662451 -0.020675967 0.08369020 2.366556e-01 2
## -11 -0.10491511 0.01858780 -0.141346525 -0.06848370 1.658541e-08 2
## -10 -0.11203303 0.02177481 -0.154710871 -0.06935518 2.674141e-07 2
## -9 0.04666754 0.03315376 -0.018312633 0.11164771 1.592468e-01 2
## -8 -0.29079100 0.01612405 -0.322393545 -0.25918845 1.042131e-72 2
## -7 -0.28495436 0.01485938 -0.314078218 -0.25583051 5.791124e-82 2
## -6 -0.29524891 0.03006193 -0.354169210 -0.23632861 9.110974e-23 2
## -5 -0.25975179 0.02803713 -0.314703543 -0.20480003 1.958731e-20 2
## -4 -0.11275999 0.03054627 -0.172629588 -0.05289040 2.229805e-04 2
## -3 -0.09273214 0.02888108 -0.149338017 -0.03612626 1.323537e-03 2
## -2 -0.12106445 0.02251499 -0.165193007 -0.07693589 7.571159e-08 2
## -1 -0.25672869 0.02229158 -0.300419385 -0.21303800 1.085101e-30 2
## 0 0.04653192 0.01623495 0.014711993 0.07835185 4.154889e-03 2
## 1 -0.91904269 0.20052686 -1.312068120 -0.52601726 4.580450e-06 2
## 2 -0.40421262 0.14524636 -0.688890263 -0.11953498 5.386796e-03 2
## 3 -1.38722222 0.08647169 -1.556703612 -1.21774083 6.450210e-58 2
## 4 0.32347930 0.10131669 0.124902230 0.52205636 1.409228e-03 2
## 5 1.57898155 0.09469659 1.393379635 1.76458347 0.000000e+00 2
## 6 -0.95589906 0.11153924 -1.174511958 -0.73728615 1.034228e-17 1
## 7 -0.87794304 0.12404083 -1.121058600 -0.63482747 1.464027e-12 1
## 8 -1.05069728 0.13921075 -1.323545335 -0.77784923 4.435907e-14 1
## 9 0.69374435 0.17159421 0.357425880 1.03006282 5.278610e-05 1
## 10 0.57494121 0.18294303 0.216379455 0.93350296 1.673781e-03 1
## 11 -1.64995609 0.18605685 -2.014620825 -1.28529136 7.445729e-19 1
## 12 -1.65126164 0.18815542 -2.020039487 -1.28248379 1.693146e-18 1
## 13 -1.93481025 0.13540613 -2.200201392 -1.66941910 2.564804e-46 1
## 14 -1.61094983 0.15749712 -1.919638514 -1.30226115 1.478824e-24 1out.mc$est.avg## ATT.avg S.E. CI.lower CI.upper p.value
## [1,] -0.5304666 0.05380099 -0.6359146 -0.4250186 6.217792e-23out.mc$est.beta## NULL::: The model cannot run F-statistics due to the same error as before. :::
Visualizing ATT by treatment timing
::: (No idea what the red line is.) I find that most estimates are below zero, showing a negative effect on accidents involving cars and vulnerable road users for most of the treatment periods but some are positive, which is known that D2 had a spike in accidents.)
plot(out.flnmc, type = "calendar", xlim = c(1, 72))
Placebo tests
::: Test 1: If t-test p-value is smaller than a pre-specified threshold (e.g. 5%), we reject the null of no-differences. Hence, the placebo test is deemed failed.(larger placebo p preferred) :::
out.flnmc <- fect(log_Count_PKW_vs_VRU ~ D, data = FNL_Accidents, index = c("Street_Group", "TimeIndex"),
force = "two-way", method = "mc", lambda = out.mc$lambda.cv,
CV = 0, parallel = TRUE, se = TRUE,
nboots = 200, placeboTest = TRUE, placebo.period = c(-2, 0))## Parallel computing ...## Bootstrapping for uncertainties ...## 200 runs::: Test:2 - TOST. The TOST checks whether the 90% confidence intervals for estimated ATTs in the placebo period exceed a pre-specified range (defined by a threshold), or the equivalence range. A TOST p-value smaller than a pre-specified threshold suggests that the null of difference bigger than the threshold is rejected; hence, the placebo test is passed.(smaller placebo TOST p-value are preferred.) :::
::: I find desirable p-values (larger for test-1 and smaller for test-2) :::
plot(out.flnmc, cex.text = 0.8, stats = c("placebo.p","equiv.p"),
main = "Estimated ATT (MC)")
### Pre-trend (or not) tests:
::: Test-1 F test. We offer a goodness-of-fit test (a variant of the test) and to gauge the presence of pretreatment (differential) trends. A larger F-test p-value suggests a better pre-trend fitting. , which I find.
Test-2 TOST. The TOST checks whether the 90% confidence intervals for estimated ATTs in the pretreatment periods (again, subject to the proportion option) exceed a pre-specified range, or the equivalence range. A smaller TOST p-value suggests a better pre-trend fitting, which I find. :::
plot(out.flnmc, type = "equiv", ylim = c(-2,2),
cex.legend = 0.6, main = "Testing Pre-Trend (MC)", cex.text = 0.8)
Exit plot
::: fect allows the treatment to switch back and forth and provides diagnostic tools for this setting. After the estimation, we can visualize the period-wise ATTs relative to the exit of treatments by setting type = “exit” (one can still draw the classic gap plot by setting type = “gap”). The x-axis is then realigned based on the timing of the treatment’s exit, not onset, e.g., 1 represents 1 period after the treatment ends.
It shows that the treatment effect hovers around zero post-treatment. :::
plot(out.flnmc, type = "exit", ylim = c(-2.5,4.5), main = "What Happens after the Treatment Switches Off?")
Tests for (no) carryover effects
::: Explanation from source: “If carryover effects do not exist, we would expect the average prediction error in those periods to be close to zero. We can treat a range of exit-treatment periods in option carryover.period to remove observations in the specified range for model fitting, and then test whether the estimated ATT in this range is significantly different from zero.
We can treat a range of exit-treatment periods in option carryover.period to remove observations in the specified range for model fitting, and then test whether the estimated ATT in this range is significantly different from zero.
Below, we set carryover.period = c(1, 3). As we deduct the treatment effect from the outcome in simdata, we expect the average prediction error for these removed periods to be close to zero.” :::
out.flnmc.c <- fect(log_Count_PKW_vs_VRU ~ D, data = FNL_Accidents, index = c("Street_Group", "TimeIndex"),
force = "two-way", method = "mc", lambda = out.mc$lambda.cv,
CV = 0, parallel = TRUE, se = TRUE,
nboots = 200, carryoverTest = TRUE, carryover.period = c(1, 15))## Parallel computing ...## Bootstrapping for uncertainties ...## 200 runs::: Like the placebo test, the plot will display the p-value of the carryover effect test (stats = “carryover.p”). Users can also add the p-value of a corresponding TOST test by setting stats = c(“carryover.p”,“equiv.p”).
I do not find substantial carry over effects. Which is perhaps good for the model fit which rests on using control periods data to estimate counterfactual outcomes. :::
plot(out.flnmc.c, type = "exit", ylim = c(-2.5,4.5),
cex.text = 0.8, main = "Carryover Effects (MC)")
::: I find that the carry over effect is around zero. So to the best of
my understanding this inference is that there are no effects carried
over after treatment ends. :::
Model-4 Effects on Log-transformed Car vs VRU Accidents using IFE or Interactive Fixed Effects Estimator
out.life <- fect(formula = log_Count_PKW_vs_VRU ~ D, # Y ~ D
data = FNL_Accidents,
index = c("Street_Group", "TimeIndex"),
method = "ife",
se = TRUE,
vartype = "jackknife",
nboots = 200,
parallel = TRUE,
CV = TRUE,
r = c(0, 5),
criterion = "moment",
force = "two-way") ## Parallel computing ...## Cross-validating ...## Criterion: Moment Conditions## Interactive fixed effects model...##
## r = 0; sigma2 = 0.20770; IC = -0.57710; PC = 0.17430; MSPE = 0.29261## *##
## r = 1; sigma2 = 0.18447; IC = 0.27338; PC = 0.60318; MSPE = 0.59456##
## r = 2; sigma2 = 0.16336; IC = 1.09537; PC = 0.93186; MSPE = 0.79695##
## r = 3; sigma2 = 0.14534; IC = 1.89659; PC = 1.18358; MSPE = 1.67937##
## r = 4; sigma2 = 0.12082; IC = 2.60439; PC = 1.27907; MSPE = 2.54517## *##
## r = 5; sigma2 = 0.10254; IC = 3.30736; PC = 1.33644; MSPE = 1.48302##
##
## r* = 4## ## Jackknife estimates ...## 6 runs## Cannot use full pre-treatment periods in the F test. The first period is removed.## ## The estimated covariance matrix is irreversible.## # Print summary of the results (ATT etc.)
summary(out.life)## Length Class Mode
## Y.dat 504 -none- numeric
## D.dat 504 -none- numeric
## I.dat 504 -none- numeric
## Y 1 -none- character
## D 1 -none- character
## X 0 -none- NULL
## T.on 504 -none- numeric
## G 0 -none- NULL
## balance.period 0 -none- NULL
## hasRevs 1 -none- numeric
## T.off 504 -none- numeric
## T.on.balance 504 -none- logical
## index 2 -none- character
## id 7 -none- character
## rawtime 72 -none- numeric
## binary 1 -none- logical
## loo 1 -none- logical
## proportion 1 -none- numeric
## pre.periods 32 -none- numeric
## tost.threshold 1 -none- numeric
## placeboTest 1 -none- logical
## placebo.period 0 -none- NULL
## carryoverTest 1 -none- logical
## carryover.period 0 -none- NULL
## unit.type 7 -none- numeric
## obs.missing 504 -none- numeric
## obs.missing.balance 1 -none- logical
## sigma2 1 -none- numeric
## sigma2.fect 1 -none- numeric
## T.on 504 -none- numeric
## Y.ct 504 -none- numeric
## eff 504 -none- numeric
## att.avg 1 -none- numeric
## att.avg.unit 1 -none- numeric
## force 1 -none- numeric
## T 1 -none- numeric
## N 1 -none- numeric
## Ntr 1 -none- numeric
## Nco 1 -none- numeric
## p 1 -none- numeric
## D 504 -none- numeric
## Y 504 -none- numeric
## X 0 -none- numeric
## I 504 -none- numeric
## II 504 -none- numeric
## tr 1 -none- numeric
## co 6 -none- numeric
## est 13 -none- list
## method 1 -none- character
## mu 1 -none- numeric
## beta 1 -none- logical
## validX 1 -none- numeric
## validF 1 -none- numeric
## niter 1 -none- numeric
## time 46 -none- numeric
## att 46 -none- numeric
## count 46 -none- numeric
## eff.calendar 144 -none- numeric
## N.calendar 72 -none- numeric
## eff.calendar.fit 144 -none- numeric
## calendar.enp 1 -none- numeric
## eff.pre 141 -none- numeric
## eff.pre.equiv 141 -none- numeric
## pre.sd 96 -none- numeric
## rmse 1 -none- numeric
## rmCV 10 -none- list
## estCV 10 -none- list
## res 504 -none- numeric
## time.off 29 -none- numeric
## att.off 29 -none- numeric
## count.off 29 -none- numeric
## eff.off 63 -none- numeric
## eff.off.equiv 63 -none- numeric
## off.sd 45 -none- numeric
## alpha 7 -none- numeric
## xi 72 -none- numeric
## alpha.tr 1 -none- numeric
## alpha.co 6 -none- numeric
## r.cv 1 -none- numeric
## IC 1 -none- numeric
## PC 1 -none- numeric
## factor 288 -none- numeric
## lambda 28 -none- numeric
## lambda.tr 4 -none- numeric
## lambda.co 24 -none- numeric
## CV.out.ife 78 -none- numeric
## est.avg 5 -none- numeric
## att.avg.boot 6 -none- numeric
## est.avg.unit 5 -none- numeric
## att.avg.unit.boot 6 -none- numeric
## est.eff.calendar 432 -none- numeric
## est.eff.calendar.fit 432 -none- numeric
## est.att 276 -none- numeric
## att.bound 92 -none- numeric
## att.boot 276 -none- numeric
## att.count.boot 276 -none- numeric
## est.att.off 174 -none- numeric
## att.off.boot 174 -none- numeric
## att.off.bound 58 -none- numeric
## att.off.count.boot 174 -none- numeric
## test.out 9 -none- list
## call 13 -none- call
## formula 3 formula callplot(out.life, main = "Estimated ATT (IFE)")
::: Pre-trend estimates are far away from zero. The model again warns
“The estimated covariance matrix is irreversible.” :::
plot(out.life, type = "calendar", xlim = c(1, 72))
# heterogenous effects box plot
plot(out.life, type = "box", xlim = c(-30, 10))
::: Unsure on the interpretation of the this heterogeneity effects box
plot. :::
out.life.p <- fect(log_Count_PKW_vs_VRU ~ D , data = FNL_Accidents, index = c("Street_Group", "TimeIndex"),
force = "two-way", method = "ife", r = 2, CV = 0,
parallel = TRUE, se = TRUE,
nboots = 200, placeboTest = TRUE, placebo.period = c(-2, 0))## Parallel computing ...## Bootstrapping for uncertainties ...## 200 runsplot(out.life.p, cex.text = 0.8, stats = c("placebo.p","equiv.p"),
main = "Estimated ATT (TWFE)")
::: This model passes the placebo tests. (A larger placebo p-value from
a t-test and a smaller placebo TOST p-value are preferred.) :::
Pre-trend (or not) tests:
::: I find that pre-trends are absent and the test shows it with the incredibly high p-value :::
plot(out.life, type = "equiv", ylim = c(-4,4),
cex.legend = 0.6, main = "Testing Pre-Trend (IFEct)", cex.text = 0.8)
Exit plot
plot(out.life, type = "exit", ylim = c(-2.5,4.5), main = "Exit Plot (IFE)")
Tests for (no) carryover effects
::: If carryover effects do not exist, we would expect the average prediction error in those periods to be close to zero.
We can treat a range of exit-treatment periods in option carryover.period to remove observations in the specified range for model fitting, and then test whether the estimated ATT in this range is significantly different from zero.
I find that carryover effects do not exist. :::
out.life.c <- fect(log_Count_PKW_vs_VRU ~ D, data = FNL_Accidents, index = c("Street_Group", "TimeIndex"),
force = "two-way", method = "ife", r = 2, CV = 0,
parallel = TRUE, se = TRUE,
nboots = 200, carryoverTest = TRUE, carryover.period = c(1, 15))## Parallel computing ...## Bootstrapping for uncertainties ...## 200 runsplot(out.life.c, type = "exit", ylim = c(-2.5,4.5),
cex.text = 0.8, main = "Carryover Effects (IFE)")
Model - 5 Gsynth, log_Count_PKW_vs_VRU
out <- gsynth(formula = NULL,
data = FNL_Accidents, # Explicitly name the data argument
Y = "log_Count_PKW_vs_VRU", # Provide outcome name as a STRING
D = "D", # Provide treatment name as a STRING
X = NULL,
na.rm = FALSE,
index = c("Street_Group", "TimeIndex"),
weight = NULL,
force = "two-way",
cl = NULL,
r = 0,
lambda = NULL,
nlambda = 10,
CV = TRUE,
criterion = "mspe",
k = 5,
EM = FALSE,
estimator = "ife",
se = FALSE,
nboots = 200,
inference = "parametric",
cov.ar = 1,
parallel = TRUE,
cores = NULL,
tol = 0.001,
seed = NULL,
min.T0 = 5,
alpha = 0.05,
normalize = FALSE)## Cross-validating ...
## r = 0; sigma2 = 0.22188; IC = -0.42400; PC = 0.18233; MSPE = 0.15437*
## r = 1; sigma2 = 0.19798; IC = 0.51560; PC = 0.74020; MSPE = 0.16466
## r = 2; sigma2 = 0.17579; IC = 1.42218; PC = 1.17084; MSPE = 0.17314
## r = 3; sigma2 = 0.15715; IC = 2.30744; PC = 1.50653; MSPE = 0.18380
## r = 4; sigma2 = 0.12505; IC = 3.04824; PC = 1.56536; MSPE = 0.19022
## r = 5; sigma2 = Inf; IC = Inf; PC = Inf; MSPE = 0.19834
##
## r* = 0out.se <- gsynth(formula = NULL,
data = FNL_Accidents, # Explicitly name the data argument
Y = "log_Count_PKW_vs_VRU", # Provide outcome name as a STRING
D = "D", # Provide treatment name as a STRING
X = NULL,
na.rm = FALSE,
index = c("Street_Group", "TimeIndex"),
weight = NULL,
force = "two-way",
cl = NULL,
r = 0,
lambda = NULL,
nlambda = 10,
CV = TRUE,
criterion = "mspe",
k = 5,
EM = FALSE,
estimator = "ife",
se = TRUE,
nboots = 200,
inference = "parametric",
cov.ar = 1,
parallel = TRUE,
cores = NULL,
tol = 0.001,
seed = NULL,
min.T0 = 5,
alpha = 0.05,
normalize = FALSE)## Parallel computing ...
## Cross-validating ...
## r = 0; sigma2 = 0.22188; IC = -0.42400; PC = 0.18233; MSPE = 0.15437*
## r = 1; sigma2 = 0.19798; IC = 0.51560; PC = 0.74020; MSPE = 0.16466
## r = 2; sigma2 = 0.17579; IC = 1.42218; PC = 1.17084; MSPE = 0.17314
## r = 3; sigma2 = 0.15715; IC = 2.30744; PC = 1.50653; MSPE = 0.18380
## r = 4; sigma2 = 0.12505; IC = 3.04824; PC = 1.56536; MSPE = 0.19022
## r = 5; sigma2 = Inf; IC = Inf; PC = Inf; MSPE = 0.19834
##
## r* = 0
##
## Simulating errors ...Bootstrapping ...
## print(out)## Call:
## gsynth.default(formula = NULL, data = FNL_Accidents, Y = "log_Count_PKW_vs_VRU",
## D = "D", X = NULL, na.rm = FALSE, index = c("Street_Group",
## "TimeIndex"), weight = NULL, force = "two-way", cl = NULL,
## r = 0, lambda = NULL, nlambda = 10, CV = TRUE, criterion = "mspe",
## k = 5, EM = FALSE, estimator = "ife", se = FALSE, nboots = 200,
## inference = "parametric", cov.ar = 1, parallel = TRUE, cores = NULL,
## tol = 0.001, seed = NULL, min.T0 = 5, alpha = 0.05, normalize = FALSE)
##
## Average Treatment Effect on the Treated:
## [1] -0.1368
##
## ~ by Period (including Pre-treatment Periods):
## 1 2 3 4 5 6 7 8
## -0.261868 -0.078765 -0.030818 0.363702 0.115093 0.200231 0.363702 -0.397023
## 9 10 11 12 13 14 15 16
## -0.146343 0.290063 0.836745 -0.020062 -0.030818 0.185295 0.464479 -0.165973
## 17 18 19 20 21 22 23 24
## -0.482160 -0.309815 0.701590 0.405588 -0.848364 -0.414583 0.673273 -0.347005
## 25 26 27 28 29 30 31 32
## -0.213921 -0.213921 0.152284 0.084706 -0.146343 0.113023 -0.560494 -0.281498
## 33 34 35 36 37 38 39 40
## -0.806479 -0.560494 -0.347005 -0.194290 -0.078765 0.084706 0.084706 -0.030818
## 41 42 43 44 45 46 47 48
## -0.146343 -0.231481 -0.384196 -0.377392 -0.444970 -0.261868 -0.030818 0.152284
## 49 50 51 52 53 54 55 56
## -0.030818 0.267808 -0.098396 -0.329445 -0.540864 -0.444970 -0.146343 -0.213921
## 57 58 59 60 61 62 63 64
## -0.008563 -0.859121 -0.030818 0.036759 0.152284 0.777853 -0.146343 -0.329445
## 65 66 67 68 69 70 71 72
## -0.530107 -0.020062 -0.068009 -0.030818 -0.329445 0.546804 0.431280 0.045445
##
## Uncertainty estimates not available.out$est.att## NULLout$est.avg## NULLout$est.beta## NULLprint(out.se)## Call:
## gsynth.default(formula = NULL, data = FNL_Accidents, Y = "log_Count_PKW_vs_VRU",
## D = "D", X = NULL, na.rm = FALSE, index = c("Street_Group",
## "TimeIndex"), weight = NULL, force = "two-way", cl = NULL,
## r = 0, lambda = NULL, nlambda = 10, CV = TRUE, criterion = "mspe",
## k = 5, EM = FALSE, estimator = "ife", se = TRUE, nboots = 200,
## inference = "parametric", cov.ar = 1, parallel = TRUE, cores = NULL,
## tol = 0.001, seed = NULL, min.T0 = 5, alpha = 0.05, normalize = FALSE)
##
## Average Treatment Effect on the Treated:
## Estimate S.E. CI.lower CI.upper p.value
## ATT.avg -0.1368 0.1727 -0.4753 0.2017 0.4282
##
## ~ by Period (including Pre-treatment Periods):
## ATT S.E. CI.lower CI.upper p.value n.Treated
## -31 -0.261868 0.7702 -1.77145 1.24772 0.73386 0
## -30 -0.078765 0.5749 -1.20559 1.04806 0.89103 0
## -29 -0.030818 0.7207 -1.44339 1.38175 0.96589 0
## -28 0.363702 0.3636 -0.34886 1.07626 0.31712 0
## -27 0.115093 0.5422 -0.94768 1.17787 0.83191 0
## -26 0.200231 0.4529 -0.68749 1.08795 0.65843 0
## -25 0.363702 0.3636 -0.34886 1.07626 0.31712 0
## -24 -0.397023 0.7867 -1.93889 1.14485 0.61378 0
## -23 -0.146343 0.7714 -1.65821 1.36552 0.84953 0
## -22 0.290063 0.4906 -0.67152 1.25165 0.55437 0
## -21 0.836745 0.3577 0.13564 1.53785 0.01933 0
## -20 -0.020062 0.4449 -0.89210 0.85198 0.96404 0
## -19 -0.030818 0.4895 -0.99030 0.92867 0.94980 0
## -18 0.185295 0.7179 -1.22182 1.59241 0.79633 0
## -17 0.464479 0.4100 -0.33907 1.26803 0.25725 0
## -16 -0.165973 0.3515 -0.85495 0.52301 0.63682 0
## -15 -0.482160 0.4686 -1.40059 0.43627 0.30350 0
## -14 -0.309815 0.4032 -1.10009 0.48046 0.44227 0
## -13 0.701590 0.4440 -0.16859 1.57177 0.11405 0
## -12 0.405588 0.8292 -1.21956 2.03073 0.62474 0
## -11 -0.848364 0.4638 -1.75744 0.06071 0.06739 0
## -10 -0.414583 0.5581 -1.50848 0.67931 0.45759 0
## -9 0.673273 0.5943 -0.49149 1.83803 0.25724 0
## -8 -0.347005 0.5581 -1.44090 0.74689 0.53411 0
## -7 -0.213921 0.4952 -1.18446 0.75662 0.66574 0
## -6 -0.213921 0.8857 -1.94994 1.52210 0.80916 0
## -5 0.152284 0.2910 -0.41814 0.72271 0.60080 0
## -4 0.084706 0.4830 -0.86190 1.03131 0.86078 0
## -3 -0.146343 0.5203 -1.16606 0.87337 0.77850 0
## -2 0.113023 0.7468 -1.35074 1.57679 0.87971 0
## -1 -0.560494 0.7606 -2.05124 0.93025 0.46118 0
## 0 -0.281498 0.4734 -1.20936 0.64636 0.55210 0
## 1 -0.806479 0.7376 -2.25210 0.63914 0.27421 1
## 2 -0.560494 0.5362 -1.61138 0.49039 0.29586 1
## 3 -0.347005 0.6562 -1.63305 0.93904 0.59691 1
## 4 -0.194290 0.3021 -0.78642 0.39784 0.52016 1
## 5 -0.078765 0.5749 -1.20559 1.04806 0.89103 1
## 6 0.084706 0.7058 -1.29869 1.46810 0.90448 1
## 7 0.084706 0.5373 -0.96845 1.13786 0.87474 1
## 8 -0.030818 0.5853 -1.17790 1.11626 0.95800 1
## 9 -0.146343 0.6294 -1.37997 1.08728 0.81614 1
## 10 -0.231481 0.4759 -1.16430 0.70134 0.62670 1
## 11 -0.384196 0.6046 -1.56926 0.80087 0.52516 1
## 12 -0.377392 0.3550 -1.07312 0.31834 0.28771 1
## 13 -0.444970 0.3955 -1.22005 0.33011 0.26050 1
## 14 -0.261868 0.4720 -1.18693 0.66319 0.57901 1
## 15 -0.030818 0.5182 -1.04643 0.98480 0.95257 1
## 16 0.152284 0.2910 -0.41814 0.72271 0.60080 1
## 17 -0.030818 0.4528 -0.91834 0.85670 0.94574 1
## 18 0.267808 0.4914 -0.69537 1.23099 0.58578 1
## 19 -0.098396 0.6721 -1.41570 1.21891 0.88361 1
## 20 -0.329445 0.4071 -1.12741 0.46852 0.41841 1
## 21 -0.540864 0.6254 -1.76656 0.68483 0.38711 1
## 22 -0.444970 0.5518 -1.52653 0.63659 0.42004 1
## 23 -0.146343 0.7150 -1.54778 1.25509 0.83783 1
## 24 -0.213921 0.4952 -1.18446 0.75662 0.66574 1
## 25 -0.008563 0.8504 -1.67528 1.65816 0.99197 1
## 26 -0.859121 0.7222 -2.27460 0.55636 0.23421 1
## 27 -0.030818 0.8113 -1.62093 1.55929 0.96970 1
## 28 0.036759 0.6481 -1.23352 1.30704 0.95477 1
## 29 0.152284 0.6146 -1.05232 1.35688 0.80431 1
## 30 0.777853 0.5039 -0.20974 1.76545 0.12266 1
## 31 -0.146343 0.4882 -1.10316 0.81047 0.76435 1
## 32 -0.329445 0.5853 -1.47662 0.81773 0.57353 1
## 33 -0.530107 0.4623 -1.43614 0.37593 0.25149 1
## 34 -0.020062 0.4344 -0.87138 0.83125 0.96316 1
## 35 -0.068009 0.5635 -1.17238 1.03636 0.90393 1
## 36 -0.030818 0.4528 -0.91834 0.85670 0.94574 1
## 37 -0.329445 0.5225 -1.35358 0.69469 0.52838 1
## 38 0.546804 0.2539 0.04922 1.04439 0.03125 1
## 39 0.431280 0.5769 -0.69945 1.56201 0.45473 1
## 40 0.045445 0.3358 -0.61264 0.70353 0.89234 1out$est.att## NULLout$est.avg## NULLout$est.beta## NULLa <- plot(out) # by default, type = "gap"## Uncertainty estimates not available.
print(a)
a.se <- plot(out.se) # by default, type = "gap"
print(a.se)
plot(out, theme.bw = FALSE) ## Uncertainty estimates not available.
plot(out.se, theme.bw = FALSE) 
plot(out, show.points = FALSE) ## Uncertainty estimates not available.
plot(out.se, show.points = FALSE)
plot(out, type = "gap", ylim = c(-2,2), xlab = "Period",
main = "Estimated ATT (FEct)")## Uncertainty estimates not available.
plot(out.se, type = "gap", ylim = c(-2,2), xlab = "Period",
main = "Estimated ATT (FEct)")
Model 6 - two-way fixed effects coutnerfacutal estimator (FEct)
Outcome variable: Count_PKW_vs_VRU
Model 6.1/No uncertainty estimates
Visualising the outcome variable
colnames(NL_accidents)## [1] "Street_Group" "Year"
## [3] "Month" "YearMonth"
## [5] "D1" "D2"
## [7] "Count_Total" "Count_Rad"
## [9] "Count_Fuss" "Count_PKW"
## [11] "Count_Krad" "Count_Gkfz"
## [13] "Count_Sonstige" "Count_PKW_vs_Rad"
## [15] "Count_PKW_vs_Fuss" "Count_PKW_vs_VRU"
## [17] "Street_Length_km" "Rate_Total_per_km"
## [19] "Rate_Rad_per_km" "Rate_Fuss_per_km"
## [21] "Rate_PKW_per_km" "Rate_Krad_per_km"
## [23] "Rate_Gkfz_per_km" "Rate_Sonstige_per_km"
## [25] "Rate_PKW_vs_Rad_per_km" "Rate_PKW_vs_Fuss_per_km"
## [27] "Rate_PKW_vs_VRU_per_km" "log_Count_Total"
## [29] "log_Count_Rad" "log_Count_Fuss"
## [31] "log_Count_PKW" "log_Count_Krad"
## [33] "log_Count_Gkfz" "log_Count_Sonstige"
## [35] "log_Count_PKW_vs_Rad" "log_Count_PKW_vs_Fuss"
## [37] "log_Count_PKW_vs_VRU" "log_Rate_Total_per_km"
## [39] "log_Rate_Rad_per_km" "log_Rate_Fuss_per_km"
## [41] "log_Rate_PKW_per_km" "log_Rate_Krad_per_km"
## [43] "log_Rate_Gkfz_per_km" "log_Rate_Sonstige_per_km"
## [45] "log_Rate_PKW_vs_Rad_per_km" "log_Rate_PKW_vs_Fuss_per_km"
## [47] "log_Rate_PKW_vs_VRU_per_km" "D"
## [49] "Date" "TimeIndex"panelview(Count_PKW_vs_VRU ~ D, data = NL_accidents, index = c("Street_Group","TimeIndex"),
axis.lab = "time", xlab = "Time", ylab = "Unit",
theme.bw = TRUE, type = "outcome",
main = "Car vs VRU Accidents (Count_PKW_vs_VRU)")## Treatment has reversals.
out.fect <- fect(Count_PKW_vs_VRU ~ D, data = NL_accidents, index = c("Street_Group","TimeIndex"),
method = "fe", force = "two-way")
summary(out.fect)## Length Class Mode
## Y.dat 792 -none- numeric
## D.dat 792 -none- numeric
## I.dat 792 -none- numeric
## Y 1 -none- character
## D 1 -none- character
## X 0 -none- NULL
## T.on 792 -none- numeric
## G 0 -none- NULL
## balance.period 0 -none- NULL
## hasRevs 1 -none- numeric
## T.off 792 -none- numeric
## T.on.balance 792 -none- logical
## index 2 -none- character
## id 11 -none- character
## rawtime 72 -none- numeric
## binary 1 -none- logical
## loo 1 -none- logical
## proportion 1 -none- numeric
## pre.periods 32 -none- numeric
## tost.threshold 1 -none- numeric
## placeboTest 1 -none- logical
## placebo.period 0 -none- NULL
## carryoverTest 1 -none- logical
## carryover.period 0 -none- NULL
## unit.type 11 -none- numeric
## obs.missing 792 -none- numeric
## obs.missing.balance 1 -none- logical
## method 1 -none- character
## Y.ct 792 -none- numeric
## Y.ct.full 792 -none- numeric
## D 792 -none- numeric
## Y 792 -none- numeric
## X 0 -none- numeric
## eff 792 -none- numeric
## I 792 -none- numeric
## II 792 -none- numeric
## att.avg 1 -none- numeric
## att.avg.unit 1 -none- numeric
## force 1 -none- numeric
## T 1 -none- numeric
## N 1 -none- numeric
## Ntr 1 -none- numeric
## Nco 1 -none- numeric
## tr 1 -none- numeric
## co 10 -none- numeric
## p 1 -none- numeric
## r.cv 1 -none- numeric
## IC 1 -none- numeric
## beta 1 -none- logical
## est 10 -none- list
## mu 1 -none- numeric
## niter 1 -none- numeric
## validX 1 -none- numeric
## validF 1 -none- numeric
## time 46 -none- numeric
## att 46 -none- numeric
## count 46 -none- numeric
## eff.calendar 144 -none- numeric
## N.calendar 72 -none- numeric
## eff.calendar.fit 144 -none- numeric
## calendar.enp 1 -none- numeric
## eff.pre 141 -none- numeric
## eff.pre.equiv 141 -none- numeric
## pre.sd 96 -none- numeric
## PC 1 -none- numeric
## sigma2 1 -none- numeric
## sigma2.fect 1 -none- numeric
## res 792 -none- numeric
## res.full 792 -none- numeric
## rmse 1 -none- numeric
## time.off 29 -none- numeric
## att.off 29 -none- numeric
## count.off 29 -none- numeric
## eff.off 63 -none- numeric
## eff.off.equiv 63 -none- numeric
## off.sd 45 -none- numeric
## alpha 11 -none- numeric
## xi 72 -none- numeric
## alpha.tr 1 -none- numeric
## alpha.co 10 -none- numeric
## call 6 -none- call
## formula 3 formula callout.fect## Call:
## fect.formula(formula = Count_PKW_vs_VRU ~ D, data = NL_accidents,
## index = c("Street_Group", "TimeIndex"), force = "two-way",
## method = "fe")
##
## ATT:
## ATT
## Tr obs. equally weighted -0.2816
## Tr units equally weighted -0.2816
##
## Uncertainty estimates not available.plot(out.fect, main = "Estimated ATT (FEct)", ylab = "Effect of D on Count_PKW_vs_VRU",
cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)## Uncertainty estimates not available.
#> Uncertainty estimates not available.Model 6.1/With uncertainty estimates
out.fect.se <- fect(formula = Count_PKW_vs_VRU ~ D, data = NL_accidents, index = c("Street_Group", "TimeIndex"), method = "fe", force = "two-way", se = TRUE,
vartype = "jackknife", # try the other one next
nboots = 200, # check if more would help
parallel = TRUE,
CV = TRUE, # Use cross-validation to choose lambda (tuning parameter for MC) [cite: 172]
nlambda = 10, # Optionally specify how many lambdas CV should try [cite: 1]
criterion = "moment",)## Parallel computing ...## Jackknife estimates ...## 10 runs## Cannot use full pre-treatment periods in the F test. The first period is removed.## ## The estimated covariance matrix is irreversible.## out.fect.se## Call:
## fect.formula(formula = Count_PKW_vs_VRU ~ D, data = NL_accidents,
## index = c("Street_Group", "TimeIndex"), force = "two-way",
## nlambda = 10, CV = TRUE, criterion = "moment", method = "fe",
## se = TRUE, vartype = "jackknife", nboots = 200, parallel = TRUE)
##
## ATT:
## ATT S.E. CI.lower CI.upper p.value
## Tr obs equally weighted -0.2816 0.1043 -0.4861 -0.07719 0.006933
## Tr units equally weighted -0.2816 0.1043 -0.4861 -0.07719 0.006933plot(out.fect.se, main = "Estimated ATT (FEct)", ylab = "Effect of D on Y",
cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p")
## Trying without CV
out.fect.se1 <- fect(formula = Count_PKW_vs_VRU ~ D, data = NL_accidents, index = c("Street_Group", "TimeIndex"), method = "fe", force = "two-way", se = TRUE,
vartype = "jackknife", # try the other one next
nboots = 200, # check if more would help
parallel = TRUE)## Parallel computing ...## Jackknife estimates ...## 10 runs## Cannot use full pre-treatment periods in the F test. The first period is removed.## ## The estimated covariance matrix is irreversible.## summary(out.fect.se1)## Length Class Mode
## Y.dat 792 -none- numeric
## D.dat 792 -none- numeric
## I.dat 792 -none- numeric
## Y 1 -none- character
## D 1 -none- character
## X 0 -none- NULL
## T.on 792 -none- numeric
## G 0 -none- NULL
## balance.period 0 -none- NULL
## hasRevs 1 -none- numeric
## T.off 792 -none- numeric
## T.on.balance 792 -none- logical
## index 2 -none- character
## id 11 -none- character
## rawtime 72 -none- numeric
## binary 1 -none- logical
## loo 1 -none- logical
## proportion 1 -none- numeric
## pre.periods 32 -none- numeric
## tost.threshold 1 -none- numeric
## placeboTest 1 -none- logical
## placebo.period 0 -none- NULL
## carryoverTest 1 -none- logical
## carryover.period 0 -none- NULL
## unit.type 11 -none- numeric
## obs.missing 792 -none- numeric
## obs.missing.balance 1 -none- logical
## method 1 -none- character
## Y.ct 792 -none- numeric
## Y.ct.full 792 -none- numeric
## D 792 -none- numeric
## Y 792 -none- numeric
## X 0 -none- numeric
## eff 792 -none- numeric
## I 792 -none- numeric
## II 792 -none- numeric
## att.avg 1 -none- numeric
## att.avg.unit 1 -none- numeric
## force 1 -none- numeric
## T 1 -none- numeric
## N 1 -none- numeric
## Ntr 1 -none- numeric
## Nco 1 -none- numeric
## tr 1 -none- numeric
## co 10 -none- numeric
## p 1 -none- numeric
## r.cv 1 -none- numeric
## IC 1 -none- numeric
## beta 1 -none- logical
## est 10 -none- list
## mu 1 -none- numeric
## niter 1 -none- numeric
## validX 1 -none- numeric
## validF 1 -none- numeric
## time 46 -none- numeric
## att 46 -none- numeric
## count 46 -none- numeric
## eff.calendar 144 -none- numeric
## N.calendar 72 -none- numeric
## eff.calendar.fit 144 -none- numeric
## calendar.enp 1 -none- numeric
## eff.pre 141 -none- numeric
## eff.pre.equiv 141 -none- numeric
## pre.sd 96 -none- numeric
## PC 1 -none- numeric
## sigma2 1 -none- numeric
## sigma2.fect 1 -none- numeric
## res 792 -none- numeric
## res.full 792 -none- numeric
## rmse 1 -none- numeric
## time.off 29 -none- numeric
## att.off 29 -none- numeric
## count.off 29 -none- numeric
## eff.off 63 -none- numeric
## eff.off.equiv 63 -none- numeric
## off.sd 45 -none- numeric
## alpha 11 -none- numeric
## xi 72 -none- numeric
## alpha.tr 1 -none- numeric
## alpha.co 10 -none- numeric
## est.avg 5 -none- numeric
## att.avg.boot 10 -none- numeric
## est.avg.unit 5 -none- numeric
## att.avg.unit.boot 10 -none- numeric
## est.eff.calendar 432 -none- numeric
## est.eff.calendar.fit 432 -none- numeric
## est.att 276 -none- numeric
## att.bound 92 -none- numeric
## att.boot 460 -none- numeric
## att.count.boot 460 -none- numeric
## est.att.off 174 -none- numeric
## att.off.boot 290 -none- numeric
## att.off.bound 58 -none- numeric
## att.off.count.boot 290 -none- numeric
## test.out 9 -none- list
## call 10 -none- call
## formula 3 formula callout.fect.se1$est.att # why is count 1 and 2? why are there 31+14 periods only## ATT S.E. CI.lower CI.upper p.value count
## -31 -0.28473413 0.3188376 -0.9096443 0.340175992 3.718367e-01 1
## -30 -0.01200686 0.2918368 -0.5839964 0.559982722 9.671824e-01 1
## -29 -0.10291595 0.3767448 -0.8413222 0.635490256 7.847214e-01 1
## -28 0.53344768 0.1480911 0.2431945 0.823700893 3.155851e-04 1
## -27 0.07890223 0.2886482 -0.4868378 0.644642241 7.845835e-01 1
## -26 0.26072041 0.2804121 -0.2888773 0.810318120 3.524872e-01 1
## -25 0.53344768 0.1480911 0.2431945 0.823700893 3.155851e-04 1
## -24 -1.01200686 0.5879713 -2.1644094 0.140395657 8.521740e-02 1
## -23 -0.19382504 0.3643019 -0.9078437 0.520193569 5.946948e-01 1
## -22 -0.01200686 0.3543780 -0.7065751 0.682561353 9.729716e-01 1
## -21 1.62435678 0.1716700 1.2878897 1.960823825 0.000000e+00 1
## -20 0.35162950 0.3090198 -0.2540382 0.957297179 2.551678e-01 1
## -19 0.16981132 0.1937875 -0.2100053 0.549627934 3.808802e-01 1
## -18 0.35162950 0.4461817 -0.5228706 1.226129636 4.306465e-01 1
## -17 1.26072041 0.3718678 0.5318728 1.989567981 6.983201e-04 1
## -16 -0.55746141 0.2820124 -1.1101955 -0.004727335 4.807251e-02 1
## -15 -0.83018868 0.3262087 -1.4695459 -0.190831412 1.092896e-02 1
## -14 -0.33018868 0.1391190 -0.6028568 -0.057520528 1.762388e-02 2
## -13 0.85162950 0.1285233 0.5997285 1.103530522 3.442846e-11 2
## -12 0.12435678 0.2832516 -0.4308061 0.679519622 6.606376e-01 2
## -11 -0.14837050 0.1754491 -0.4922444 0.195503421 3.977418e-01 2
## -10 -0.23927959 0.1542827 -0.5416681 0.063108885 1.209216e-01 2
## -9 0.57890223 0.1713854 0.2429930 0.914811428 7.307340e-04 2
## -8 -0.51200686 0.1988157 -0.9016784 -0.122335311 1.001577e-02 2
## -7 -0.55746141 0.2405827 -1.0289948 -0.085928008 2.049664e-02 2
## -6 -0.10291595 0.2608998 -0.6142702 0.408438267 6.932374e-01 2
## -5 0.16981132 0.1651667 -0.1539094 0.493532043 3.038929e-01 2
## -4 -0.14837050 0.2777425 -0.6927357 0.395994733 5.932020e-01 2
## -3 -0.73927959 0.1850467 -1.1019645 -0.376594660 6.466785e-05 2
## -2 -0.14837050 0.2628209 -0.6634900 0.366749054 5.723930e-01 2
## -1 -0.23927959 0.2960611 -0.8195488 0.340989582 4.189697e-01 2
## 0 -0.28473413 0.1553658 -0.5892454 0.019777158 6.685169e-02 2
## 1 -0.26320755 0.2883751 -0.8284123 0.301997199 3.613864e-01 2
## 2 -0.36320755 0.2228903 -0.8000646 0.073649457 1.031998e-01 2
## 3 -0.36320755 0.2555358 -0.8640485 0.137633415 1.552131e-01 2
## 4 -0.46320755 0.1251197 -0.7084376 -0.217977448 2.138094e-04 2
## 5 0.03679245 0.1849551 -0.3257129 0.399297833 8.423203e-01 2
## 6 0.38679245 0.2639599 -0.1305595 0.904144446 1.428268e-01 1
## 7 0.18679245 0.2856267 -0.3730256 0.746610497 5.131285e-01 1
## 8 0.18679245 0.2143436 -0.2333134 0.606898279 3.835017e-01 1
## 9 -0.51320755 0.3075054 -1.1159070 0.089491897 9.512957e-02 1
## 10 -0.21320755 0.2505964 -0.7043674 0.277952327 3.948801e-01 1
## 11 -1.01320755 0.3661658 -1.7308793 -0.295535828 5.656202e-03 1
## 12 -0.31320755 0.1869142 -0.6795526 0.053137506 9.380173e-02 1
## 13 -0.71320755 0.3733943 -1.4450469 0.018631778 5.612477e-02 1
## 14 -0.51320755 0.3208509 -1.1420637 0.115648646 1.097050e-01 1out.fect.se1$est.avg## ATT.avg S.E. CI.lower CI.upper p.value
## [1,] -0.2816286 0.1043052 -0.4860631 -0.07719412 0.006933049out.fect.se1$est.beta## NULLplot(out.fect.se1, main = "Estimated ATT (FEct)", ylab = "Effect of D on Y",
cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p")
Model - 7 (TWFE simple)
PLM
Using “D”
# Convert data to panel data structure
pdata <- pdata.frame(NL_accidents, index = c("Street_Group", "TimeIndex"))
# Estimate TWFE with no additional covariates
model_plm <- plm(Count_PKW_vs_VRU ~ D, data = pdata, model = "within", effect = "twoways")
summary(model_plm)## Twoways effects Within Model
##
## Call:
## plm(formula = Count_PKW_vs_VRU ~ D, data = pdata, effect = "twoways",
## model = "within")
##
## Balanced Panel: n = 11, T = 72, N = 792
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.512183 -0.603442 -0.089553 0.466090 5.477716
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## D -0.28163 0.30077 -0.9364 0.3494
##
## Total Sum of Squares: 816.49
## Residual Sum of Squares: 815.49
## R-Squared: 0.0012351
## Adj. R-Squared: -0.11428
## F-statistic: 0.876773 on 1 and 709 DF, p-value: 0.34941# Estimate TWFE model
model_fixest <- feols(Count_PKW_vs_VRU ~ D | Street_Group + TimeIndex, data = NL_accidents)
summary(model_fixest)## OLS estimation, Dep. Var.: Count_PKW_vs_VRU
## Observations: 792
## Fixed-effects: Street_Group: 11, TimeIndex: 72
## Standard-errors: Clustered (Street_Group)
## Estimate Std. Error t value Pr(>|t|)
## D -0.281629 0.108855 -2.5872 0.02708 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.01472 Adj. R2: 0.347108
## Within R2: 0.001235modelsummary(list(model_fixest,model_plm))| (1) | (2) | |
|---|---|---|
| D | −0.282 | −0.282 |
| (0.109) | (0.301) | |
| Num.Obs. | 792 | 792 |
| R2 | 0.415 | 0.001 |
| R2 Adj. | 0.347 | −0.114 |
| R2 Within | 0.001 | |
| R2 Within Adj. | 0.000 | |
| AIC | 2436.7 | 2274.7 |
| BIC | 2824.7 | 2284.1 |
| RMSE | 1.01 | 1.01 |
| Std.Errors | by: Street_Group | |
| FE: Street_Group | X | |
| FE: TimeIndex | X |
plotting both Models-7
# Step 1: Tidy up both model outputs
tidy_plm <- broom::tidy(model_plm) %>%
mutate(model = "plm")
tidy_fixest <- broom::tidy(model_fixest, conf.int = TRUE) %>%
mutate(model = "fixest")
# Step 2: Combine into one dataframe
coef_df <- bind_rows(tidy_plm, tidy_fixest) %>%
filter(term == "D") # We're only plotting coefficient on D
# Step 3: Plot
ggplot(coef_df, aes(x = model, y = estimate)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = estimate - std.error * 1.96,
ymax = estimate + std.error * 1.96),
width = 0.2) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Comparison of TWFE Estimates from plm and fixest",
y = "Coefficient Estimate (D)",
x = "Model") +
theme_minimal()
### observed vs predicted values (fixest)
# Step 1: Add predicted values to the data
NL_accidents$predicted_fixest <- predict(model_fixest)
# Step 2: Aggregate observed and predicted values by TimeIndex
agg_df <- NL_accidents %>%
group_by(TimeIndex) %>%
summarise(observed = mean(Count_PKW_vs_VRU),
predicted = mean(predicted_fixest))
# Step 3: Plot
ggplot(agg_df, aes(x = TimeIndex)) +
geom_line(aes(y = observed), color = "black", linewidth = 1.2, linetype = "solid") +
geom_line(aes(y = predicted), color = "blue", linewidth = 3, linetype = "dashed") +
labs(title = "Observed vs Predicted Values (fixest model)",
y = "Count_PKW_vs_VRU",
x = "TimeIndex",
caption = "Dashed = Predicted, Solid = Observed") +
theme_minimal()
### time-index plot/event study / Estimates Over Time (TWFE)
# Step 1: Run interaction model (event study style)
event_model <- feols(Count_PKW_vs_VRU ~ i(TimeIndex, D, ref = 1) | Street_Group + TimeIndex,
data = NL_accidents)## The variables 'TimeIndex::2:D', 'TimeIndex::3:D', 'TimeIndex::4:D', 'TimeIndex::5:D', 'TimeIndex::6:D', 'TimeIndex::7:D' and 46 others have been removed because of collinearity (see $collin.var).# Step 2: Plot using built-in fixest plotting
iplot(event_model,
main = "Event Study: Effect of D Over TimeIndex",
xlab = "TimeIndex (relative to ref = 1)",
ylab = "Estimate (D × TimeIndex)")
event_model$collin.var## [1] "TimeIndex::2:D" "TimeIndex::3:D" "TimeIndex::4:D" "TimeIndex::5:D"
## [5] "TimeIndex::6:D" "TimeIndex::7:D" "TimeIndex::8:D" "TimeIndex::9:D"
## [9] "TimeIndex::10:D" "TimeIndex::11:D" "TimeIndex::12:D" "TimeIndex::13:D"
## [13] "TimeIndex::14:D" "TimeIndex::15:D" "TimeIndex::16:D" "TimeIndex::17:D"
## [17] "TimeIndex::18:D" "TimeIndex::19:D" "TimeIndex::20:D" "TimeIndex::21:D"
## [21] "TimeIndex::22:D" "TimeIndex::23:D" "TimeIndex::24:D" "TimeIndex::25:D"
## [25] "TimeIndex::26:D" "TimeIndex::27:D" "TimeIndex::28:D" "TimeIndex::29:D"
## [29] "TimeIndex::30:D" "TimeIndex::31:D" "TimeIndex::32:D" "TimeIndex::47:D"
## [33] "TimeIndex::48:D" "TimeIndex::49:D" "TimeIndex::50:D" "TimeIndex::51:D"
## [37] "TimeIndex::52:D" "TimeIndex::53:D" "TimeIndex::54:D" "TimeIndex::55:D"
## [41] "TimeIndex::56:D" "TimeIndex::57:D" "TimeIndex::58:D" "TimeIndex::59:D"
## [45] "TimeIndex::60:D" "TimeIndex::61:D" "TimeIndex::67:D" "TimeIndex::68:D"
## [49] "TimeIndex::69:D" "TimeIndex::70:D" "TimeIndex::71:D" "TimeIndex::72:D"Model-8 TWFE with D1 and D2
# Estimate TWFE model
fixest_2D <- feols(Count_PKW_vs_VRU ~ D1 + D2 | Street_Group + TimeIndex, data = pdata)
summary(fixest_2D)## OLS estimation, Dep. Var.: Count_PKW_vs_VRU
## Observations: 792
## Fixed-effects: Street_Group: 11, TimeIndex: 72
## Standard-errors: Clustered (Street_Group)
## Estimate Std. Error t value Pr(>|t|)
## D1 -0.370350 0.123487 -2.999106 0.013364 *
## D2 -0.033208 0.136507 -0.243266 0.812717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.01448 Adj. R2: 0.346491
## Within R2: 0.001701# Step 1: Name your models
models <- list(
"PLM TWFE" = model_plm,
"Fixest TWFE (D)" = model_fixest,
"Fixest TWFE (D1 & D2)" = fixest_2D
)
# Step 2: Custom coefficient labels
coef_labels <- c(
"D" = "Treatment (D)",
"D1" = "Treatment Period 1",
"D2" = "Treatment Period 2"
)
# Step 3: Custom goodness-of-fit (GOF) statistics
gof_map_custom <- tribble(
~raw, ~clean, ~fmt,
"r.squared", "R²", 3,
"adj.r.squared", "Adjusted R²", 3,
"nobs", "Observations", 0,
"FE: Street_Group", "Street FE", 0,
"FE: TimeIndex", "Time FE", 0
)
# Step 4: Modelsummary table with all extras
modelsummary(
models,
coef_map = coef_labels,
gof_map = gof_map_custom,
stars = TRUE,
notes = list(
"PLM uses default SEs.",
"Fixest models use clustered SEs by Street_Group.",
"All models include two-way fixed effects (Street_Group, TimeIndex)."
),
title = "Comparison of TWFE Models: Treatment Effects"
)| PLM TWFE | Fixest TWFE (D) | Fixest TWFE (D1 & D2) | |
|---|---|---|---|
| Treatment (D) | −0.282 | −0.282* | |
| (0.301) | (0.109) | ||
| Treatment Period 1 | −0.370* | ||
| (0.123) | |||
| Treatment Period 2 | −0.033 | ||
| (0.137) | |||
| R² | 0.001 | 0.415 | 0.415 |
| Adjusted R² | −0.114 | 0.347 | 0.346 |
| Observations | 792 | 792 | 792 |
| Street FE | X | X | |
| Time FE | X | X | |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| PLM uses default SEs. | |||
| Fixest models use clustered SEs by Street_Group. | |||
| All models include two-way fixed effects (Street_Group, TimeIndex). |