Replication R Code for
Crime under-reporting in Bogotá: a spatial panel model with fixed effects
Chanci, Kumbhakar, and Sandoval, 2023. Crime under-reporting in Bogotá: a spatial panel model with fixed effects. Empirical Economics . Doi: s00181-023-02517-4
- Version : August 2023 (using R 4.3.1)
- By : Luis Chancí
- Email : luischanci@santotomas.cl
1. Setup
Working directory and packages. See footnote 17 on the paper for more details.
Pkg <- c("spdep" , "spatialreg" , "tmap", "sf", "haven","tidyverse",
"rgdal", "splm" , "plm" , "frontier" , "texreg")
invisible(lapply(Pkg, library, character.only = TRUE))
setwd("D:/OneDrive - Corporación Santo Tomas/2_Research/ProyectoInterno/Estimation")1.1. Spatial Weights Matrix: the SWM created after editing shapefiles in ArcGIS
# The SWM (created using Arcgis)
W <- './gis/ArcGISpro_project/Wn3_for_R.txt' %>%
read.table %>%
as.matrix
WQ <- './gis/ArcGISpro_project/WnQ_for_R.txt' %>%
read.table %>%
as.matrix
nobs <- dim(W)[1]1.2. Dataset: Loading the dataset prepared from crime reports processed in Stata
# Dataset (previously cleaned using Stata):
Data <- './data/regression_data_cuadrante.dta' %>%
read_dta %>%
as.data.frame
Crimes <- c('TheftHo', # Residential Burglary
'PersInj', # Personal Injuries
'Homicid', # Homicides
'TheftPe', # Theft and Robbery
'Extorsi', # Extortion
'Sexualc' # Sexual Assault
)
aCr <- paste0( "asinh_", Crimes) 2. Descriptive Statistics (Table 1)
2.1. Summary Statistics
Data %>%
select(all_of(Crimes)) %>%
na.omit %>%
summary(., digits=1)
## TheftHo PersInj Homicid TheftPe Extorsi
## Min. : 0 Min. : 0 Min. : 0.0 Min. : 0 Min. : 0.0
## 1st Qu.: 1 1st Qu.: 4 1st Qu.: 0.0 1st Qu.: 13 1st Qu.: 0.0
## Median : 3 Median : 8 Median : 0.5 Median : 29 Median : 0.1
## Mean : 4 Mean : 11 Mean : 1.0 Mean : 45 Mean : 0.4
## 3rd Qu.: 5 3rd Qu.: 14 3rd Qu.: 1.3 3rd Qu.: 57 3rd Qu.: 0.5
## Max. :136 Max. :516 Max. :16.0 Max. :622 Max. :11.0
## Sexualc
## Min. : 0.0
## 1st Qu.: 0.4
## Median : 1.0
## Mean : 1.8
## 3rd Qu.: 2.2
## Max. :183.02.2. Explorign spatial correlations (Moran’s I)
dta.m <- Data %>%
select(uniqueid, year, Crimes) %>%
drop_na %>%
group_by(uniqueid) %>%
summarise_at(vars(Crimes), list(m = mean))
# Using Theoretical Values
lapply(dta.m[,2:7], moran.test, mat2listw(W), alternative="greater")
## $TheftHo_m
##
## Moran I test under randomisation
##
## data: X[[i]]
## weights: mat2listw(W)
##
## Moran I statistic standard deviate = 35.068, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 2.387576e-01 -9.541985e-04 4.672689e-05
##
##
## $PersInj_m
##
## Moran I test under randomisation
##
## data: X[[i]]
## weights: mat2listw(W)
##
## Moran I statistic standard deviate = 10.148, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 6.716000e-02 -9.541985e-04 4.505075e-05
##
##
## $Homicid_m
##
## Moran I test under randomisation
##
## data: X[[i]]
## weights: mat2listw(W)
##
## Moran I statistic standard deviate = 46.555, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 3.211545e-01 -9.541985e-04 4.787028e-05
##
##
## $TheftPe_m
##
## Moran I test under randomisation
##
## data: X[[i]]
## weights: mat2listw(W)
##
## Moran I statistic standard deviate = 28.99, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 1.993324e-01 -9.541985e-04 4.773196e-05
##
##
## $Extorsi_m
##
## Moran I test under randomisation
##
## data: X[[i]]
## weights: mat2listw(W)
##
## Moran I statistic standard deviate = 18.954, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 1.295184e-01 -9.541985e-04 4.738243e-05
##
##
## $Sexualc_m
##
## Moran I test under randomisation
##
## data: X[[i]]
## weights: mat2listw(W)
##
## Moran I statistic standard deviate = 3.0741, p-value = 0.001056
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 1.144663e-02 -9.541985e-04 1.627279e-05
# Using Simulated Values
lapply(dta.m[,2:7], moran.mc , mat2listw(W), nsim=999, alternative="greater")
## $TheftHo_m
##
## Monte-Carlo simulation of Moran I
##
## data: X[[i]]
## weights: mat2listw(W)
## number of simulations + 1: 1000
##
## statistic = 0.23876, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
##
##
## $PersInj_m
##
## Monte-Carlo simulation of Moran I
##
## data: X[[i]]
## weights: mat2listw(W)
## number of simulations + 1: 1000
##
## statistic = 0.06716, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
##
##
## $Homicid_m
##
## Monte-Carlo simulation of Moran I
##
## data: X[[i]]
## weights: mat2listw(W)
## number of simulations + 1: 1000
##
## statistic = 0.32115, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
##
##
## $TheftPe_m
##
## Monte-Carlo simulation of Moran I
##
## data: X[[i]]
## weights: mat2listw(W)
## number of simulations + 1: 1000
##
## statistic = 0.19933, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
##
##
## $Extorsi_m
##
## Monte-Carlo simulation of Moran I
##
## data: X[[i]]
## weights: mat2listw(W)
## number of simulations + 1: 1000
##
## statistic = 0.12952, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
##
##
## $Sexualc_m
##
## Monte-Carlo simulation of Moran I
##
## data: X[[i]]
## weights: mat2listw(W)
## number of simulations + 1: 1000
##
## statistic = 0.011447, observed rank = 987, p-value = 0.013
## alternative hypothesis: greater3. Spatial Stochastic Panel Model with F.E.
The main function has as inputs a criminal offense in the panel dataset and a spatial weights matrix. Then, the function reports estimates of: - the parameters, \(\rho\), \(\xi\), \(\sigma_{\dot{u}}\), and \(\sigma_{\dot{v}}\), as described in equations (3), (4), and (5) in the paper,
\[\begin{eqnarray} y_{it} &=& \alpha_i + \rho\sum_{j=1}^N w_{ij}y_{jt} + \boldsymbol{x}_{it}'\boldsymbol{\beta} +\varepsilon_{it} \\ \varepsilon_{it}&=& \nu_{it} - u_{it} \end{eqnarray}\]
where,
\[\begin{eqnarray} u_{it} = \xi\sum_{j=1}^N w_{ij}\dot{u}_{jt} + \dot{u}_{it} \hspace{0.5cm};\hspace{0.5cm} \nu_{it} = \xi\sum_{j=1}^N w_{ij}\dot{\nu}_{jt} + \dot{\nu}_{it} \end{eqnarray}\]
It also reports the (average) shares of Direct and Indirect (Spillover) Effects, as described in equation (11).
##############################
# (Function to conduct the 2-steps estimation)
main.fn <- function(Offense,W){
df <- Data %>%
select(uniqueid, year, all_of(Offense)) %>%
drop_na
nt <- nrow(df)
tobs <- nt/nobs
df <- df[order(df$year,df$uniqueid),] # sort by t and i (spatial panel)
Y <- matrix(df[[Offense]],nt,1)
df$Y <- Y
SWMp <- kronecker(diag(tobs) , W)
df$Wy<- (SWMp%*%Y)
X <- model.matrix(~factor(year)-1,df)
df <- cbind(df,X[,-1])
# 1. First stage:
# 1.1. Transformation Approach for FE (Lee and Yu, 2010)
df <- df[order(df$uniqueid,df$year),] # sort by i and t (panel)
Q <- diag(tobs)-(1/tobs)*matrix(1,tobs,tobs)
F <- kronecker(diag(nobs) , eigen(Q)$vectors[,1:(tobs-1)])
Fy <- t(F)%*%as.matrix(df$Y)
FWy <- t(F)%*%as.matrix(df$Wy)
FX <- t(F)%*%as.matrix(df%>%select(starts_with("factor")))
df.2 <- df|>select(uniqueid,year)
df.2 <- df.2[order(df.2$uniqueid,df.2$year),] # sort by i and t (panel)
Tyr <- df.2$year[tobs] # latest year (t=T)
df.2 <- subset(df.2, year!=Tyr) # Drop obs t=T
df.2 <- cbind(df.2,Fy,FWy,FX)
df.2 <- df.2[order(df.2$year,df.2$uniqueid),] # sort by t and i (spatial panel)
# 1.2. Estimation (first-stage)
SWMp.F<- kronecker(diag(tobs-1), W)
Stage1<-summary(
stsls(Fy ~ . - 1,
data = df.2%>%select(-c('uniqueid','year','FWy')),
mat2listw(SWMp.F)))
Rho <- Stage1$coefficients[1]
Beta <- Stage1$coefficients[2:length(Stage1$coefficients)]
#Note:
# Alternatively, one could use 2SLS using as instruments: Z=cbind(SWMp%*%FX,SWMp%*%(SWMp%*%FX))
# and the param are solve(t(cbind(FWy,FX))%*%PH%*%(cbind(FWy,FX)))%*%(t(cbind(FWy,FX))%*%PH%*%Fy)
# where PH <- Z%*%solve(t(Z)%*%Z)%*%t(Z) , but the command runds faster
# 2. Second stage:
name.b <- gsub("`","",names(Beta))
df <- df[order(df$year,df$uniqueid),]
X.used <- as.matrix(as.data.frame(df)%>%select(all_of(name.b)))
Xb <- (X.used)%*%as.matrix(Beta)
rit <- df$Y - (Rho*df$Wy) - Xb
df.s2 <- cbind(df%>%select(uniqueid,year,Y),rit,Xb,X.used) # Dataframe Stage 2
df.s2 <- df.s2[order(df.s2$uniqueid,df.s2$year),] # sort by i and t (panel)
df.s2 <- df.s2 |> group_by(uniqueid) |> mutate(alphaStar=mean(rit)) # alpha star
df.s2 <- df.s2 |> group_by(uniqueid) |> mutate(ri=rit - mean(rit)) # rtilde
E.eta.1<- mean(df.s2$rit) # E(alpha_i+beta_0+v-u)
# 2.1. Estimation (second stage)
empty <- array(numeric(), c(nobs,nobs,tobs))
df.s2 <- df.s2[order(df.s2$year,df.s2$uniqueid),] # sort by t and i (spanel)
temp.1 <- split(df.s2$ri, ceiling(seq_along(df.s2$ri)/nobs))
for (tt in 1:tobs) {
empty[,,tt] <- matrix(temp.1[[tt]])%*%t(as.matrix(temp.1[[tt]]))
}
E.eta.2 <- rowSums(empty, dims = 2)/tobs
obj.fn <- function(parm){
.xi <- exp(parm[1])
.sv <- exp(parm[2])
.su <- exp(parm[3])
norm(E.eta.2-((.sv^2+.su^2)*(1-2/pi))*(diag(nobs)+.xi*W)%*%t(diag(nobs)+.xi*W),"F")
}
theta <- optim(c(-3,0.01,0.01), obj.fn) # guess for xi based on SAC
xi <- exp(theta$par[1])
sv <- exp(theta$par[2])
su <- exp(theta$par[3])
sv2 <- sv^2
su2 <- su^2
hat.E.u <- su*sqrt(2/pi)*(1+xi)
df.s2$FE <- df.s2$alphaStar + hat.E.u
df.s2 <- df.s2[order(df.s2$year,df.s2$uniqueid),] # sort by t and i (spanel)
df.s2$vu<- (df.s2$ri - hat.E.u) # v-u
#2.2. s.e. (bootstrap)
set.seed(123)
B <- 50
.twop <- c(-(sqrt(5)-1)/2, (sqrt(5)+1)/2)
.prob <- rev(abs(.twop)/sqrt(5))
Srho <- (solve(diag(nobs) - Rho*W))
A <- (diag(nobs) + xi*W)
SA <- kronecker(diag(tobs),Srho%*%A)
SA.F <- kronecker(diag(tobs-1),Srho%*%A)
df.2.b <- df.2 |> select(c(uniqueid,year,all_of(name.b)))
df.2.b <- df.2.b[order(df.2.b$year, df.2.b$uniqueid),] # sort t,i (spanel)
F.e.tild<- kronecker(diag(tobs-1),solve(A))%*%as.matrix(Stage1$residuals)
FX.b <- as.matrix(as.data.frame(df.2.b)|>select(all_of(name.b)))
Fy.temp <- (kronecker(diag(tobs-1),Srho))%*%( FX.b%*%as.matrix(Beta) )
df.2.b$Fytemp <- Fy.temp
df.s2.b <- df.s2 |> select(c(uniqueid,year,FE,vu,Xb,all_of(name.b)))
df.s2.b <- df.s2.b[order(df.s2.b$year,df.s2.b$uniqueid),] # sort t,i (spanel)
e.tild <- kronecker(diag(tobs),solve(A))%*%as.matrix(df.s2.b$vu)
df.s2.b$vuT<- e.tild
y.tempNT <- (kronecker(diag(tobs),Srho))%*%(as.matrix(df.s2.b$FE+df.s2.b$Xb))
df.s2.b$yte<- y.tempNT
out.b <- matrix(0,4,B)
for (bb in 1:B) {
e.b <- e.tild*sample(.twop, size=length(e.tild), replace=T, prob=.prob)
df.s2.b <- df.s2.b[order(df.s2.b$year,df.s2.b$uniqueid),] # sort t,i (spanel)
df.s2.b$eB<- e.b
df.s2.b$Yb<- df.s2.b$yte + SA%*%as.matrix(df.s2.b$eB)
df.s2.b <-df.s2.b[order(df.s2.b$uniqueid,df.s2.b$year),] # sort i,t (panel)
FYb <- t(F)%*%as.matrix(df.s2.b$Yb)
df.2.b <- df.2.b[order(df.2.b$uniqueid,df.2.b$year),] # sort i,t (panel)
df.2.b$FYb<-FYb
df.2.b <- df.2.b[order(df.2.b$year,df.2.b$uniqueid),] # sort t,i (spanel)
Stage1.b<-summary(
stsls(FYb ~ . - 1,
data = df.2.b%>%select(c(FYb,all_of(name.b))),
mat2listw(SWMp.F,style='C')))
Rho.b <- Stage1.b$coefficients[1]
Beta.b <- Stage1.b$coefficients[2:length(Stage1.b$coefficients)]
df.s2.b <- df.s2.b[order(df.s2.b$year, df.s2.b$uniqueid),] # sort t,i (spanel)
df.s2.b$ritb <- df.s2.b$Yb -
(Rho.b*(SWMp%*%as.matrix(df.s2.b$Yb))) -
(as.matrix(as.data.frame(df.s2.b)%>%select(all_of(name.b))))%*%as.matrix(Beta.b)
df.s2.b <- df.s2.b[order(df.s2.b$uniqueid,df.s2.b$year),] # sort i,t (panel)
df.s2.b <- df.s2.b |> group_by(uniqueid) |> mutate(rib = ritb - mean(ritb))
E.eta.1.b <- mean(df.s2.b$ritb)
empty <- array(numeric(), c(nobs,nobs,tobs))
df.s2.b <- df.s2.b[order(df.s2.b$year, df.s2.b$uniqueid),] # sort t,i (spanel)
temp.1 <- split(df.s2.b$rib, ceiling(seq_along(df.s2.b$rib)/nobs))
for (tt in 1:tobs) {
empty[,,tt] <- matrix(temp.1[[tt]])%*%t(as.matrix(temp.1[[tt]]))
}
E.eta.2.b <- rowSums(empty, dims = 2)/tobs
obj.fn.b <- function(parm){
.xi <- exp(parm[1])
.sv <- exp(parm[2])
.su <- exp(parm[3])
norm(E.eta.2.b-((.sv^2+.su^2)*(1-2/pi))*(diag(nobs)+.xi*W)%*%t(diag(nobs)+.xi*W),"F")
}
theta.b <- optim(c(-3,0.01,0.01), obj.fn.b)
xi.b <- exp(theta.b$par[1])
sv.b <- exp(theta.b$par[2])
su.b <- exp(theta.b$par[3])
sv2.b <- (sv.b)^2
su2.b <- (su.b)^2
out.b[,bb] <- c(Rho.b,xi.b,su2.b,sv2.b)
}
se.b <- c(sd(out.b[1,]), sd(out.b[2,]), sd(out.b[3,]), sd(out.b[4,]))
p.moran <- moran.test(e.tild, mat2listw(SWMp,style='C'), alternative="greater")$p.val
#2.3. Underreporting
mu.start <- -(su2)/(sv2+su2)*e.tild
s2.start <- (su2*sv2)/(sv2+su2)
ait <- -(mu.start/sqrt(s2.start))
phi <- dnorm(ait)
Phi <- (1-pnorm(ait))
E.u.tild <- mu.start + sqrt(s2.start)*(phi/Phi)
uhat <- kronecker(diag(tobs),A)%*%E.u.tild
df.s2 <- df.s2[order(df.s2$year,df.s2$uniqueid),] # sort t,i (spanel)
df <- df[order(df$year, df$uniqueid),] # sort t,i (spanel)
Diag <- diag(diag(Srho))
T.eff <- kronecker(diag(tobs) , Srho)%*%uhat
D.eff <- (kronecker(diag(tobs), Diag)%*%uhat)/T.eff
I.eff <- (kronecker(diag(tobs), (Srho - Diag))%*%uhat)/T.eff
#3.Output
.name <- c("Rho", "Xi")
.para <- c( Rho , xi)
.se <- c(se.b[1], se.b[2])
.pval <- 2*(1-pnorm(abs(.para/.se)))
Re <- createTexreg(
coef.names = .name,
coef = .para,
se = .se,
pvalues = .pval,
gof = c(sv2, su2, p.moran, nt),
gof.names = c("Sigma2_v", "Sigma2_u", "Res. Moran's I (pval)", "Num. obs."),
gof.decimal = c(TRUE,TRUE,TRUE,FALSE),
model.name = c(substring(Offense, first=7))
)
return(list(Reg=Re, Shares=matrix(c(mean(D.eff),sd(D.eff),mean(I.eff),sd(I.eff)),2,2)))
}
##############################4. Main Results (Tables 2)
Results in Panel A Table 2:
########################
# Coeff in Panel A - T2
out.W <- lapply(aCr, main.fn, W)
lapply(1:length(aCr),
function(x) screenreg(out.W[[x]]$Reg, stars = c(0.01, 0.05, 0.1), digits = 3, dcolumn = TRUE))
## [[1]]
##
## ===================================
## TheftHo
## -----------------------------------
## Rho 0.915 ***
## (0.158)
## Xi 0.185
## (0.635)
## -----------------------------------
## Sigma2_v 0.291
## Sigma2_u 0.091
## Res. Moran's I (pval) 0.205
## Num. obs. 7343
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[2]]
##
## ===================================
## PersInj
## -----------------------------------
## Rho 0.562 **
## (0.249)
## Xi 0.582 ***
## (0.125)
## -----------------------------------
## Sigma2_v 0.187
## Sigma2_u 0.193
## Res. Moran's I (pval) 1.000
## Num. obs. 8392
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[3]]
##
## ===================================
## Homicid
## -----------------------------------
## Rho 0.792 ***
## (0.138)
## Xi 0.152 ***
## (0.028)
## -----------------------------------
## Sigma2_v 0.223
## Sigma2_u 0.086
## Res. Moran's I (pval) 0.217
## Num. obs. 9441
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[4]]
##
## ==================================
## TheftPe
## ----------------------------------
## Rho 0.241
## (0.438)
## Xi 1.220 **
## (0.498)
## ----------------------------------
## Sigma2_v 0.223
## Sigma2_u 0.009
## Res. Moran's I (pval) 1.000
## Num. obs. 6294
## ==================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[5]]
##
## ===================================
## Extorsi
## -----------------------------------
## Rho 0.933 ***
## (0.100)
## Xi 0.000
## (0.206)
## -----------------------------------
## Sigma2_v 0.083
## Sigma2_u 0.115
## Res. Moran's I (pval) 1.000
## Num. obs. 6294
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[6]]
##
## ===================================
## Sexualc
## -----------------------------------
## Rho 0.819 ***
## (0.146)
## Xi 0.061
## (0.070)
## -----------------------------------
## Sigma2_v 0.194
## Sigma2_u 0.203
## Res. Moran's I (pval) 0.535
## Num. obs. 6294
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.16. Robustness check (Table 5)
Robustness checks using a Queen SWM.
# Table 5:
# Robustness check: Queen
out.WQ <- lapply(aCr, main.fn, WQ)
lapply(1:length(aCr),
function(x) screenreg(out.WQ[[x]]$Reg, stars = c(0.01, 0.05, 0.1), digits = 3, dcolumn = TRUE))
## [[1]]
##
## ===================================
## TheftHo
## -----------------------------------
## Rho 0.831 ***
## (0.078)
## Xi 0.000
## (0.000)
## -----------------------------------
## Sigma2_v 0.332
## Sigma2_u 0.000
## Res. Moran's I (pval) 1.000
## Num. obs. 7343
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[2]]
##
## ===================================
## PersInj
## -----------------------------------
## Rho 0.993 ***
## (0.013)
## Xi 0.000
## (0.037)
## -----------------------------------
## Sigma2_v 0.041
## Sigma2_u 0.301
## Res. Moran's I (pval) 1.000
## Num. obs. 8392
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[3]]
##
## ===================================
## Homicid
## -----------------------------------
## Rho 0.869 ***
## (0.084)
## Xi 0.000
## (0.029)
## -----------------------------------
## Sigma2_v 0.000
## Sigma2_u 0.267
## Res. Moran's I (pval) 1.000
## Num. obs. 9441
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[4]]
##
## ===================================
## TheftPe
## -----------------------------------
## Rho 0.696 ***
## (0.211)
## Xi 0.000
## (0.148)
## -----------------------------------
## Sigma2_v 0.010
## Sigma2_u 0.192
## Res. Moran's I (pval) 1.000
## Num. obs. 6294
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[5]]
##
## ===================================
## Extorsi
## -----------------------------------
## Rho 0.254
## (0.242)
## Xi 0.281 ***
## (0.019)
## -----------------------------------
## Sigma2_v 0.153
## Sigma2_u 0.039
## Res. Moran's I (pval) 0.932
## Num. obs. 6294
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[6]]
##
## ===================================
## Sexualc
## -----------------------------------
## Rho 1.063 ***
## (0.003)
## Xi 0.000
## (0.160)
## -----------------------------------
## Sigma2_v 0.091
## Sigma2_u 0.265
## Res. Moran's I (pval) 1.000
## Num. obs. 6294
## ===================================
## *** p < 0.01; ** p < 0.05; * p < 0.17. Other results
Other models in Table 1 (panels B and C for the SAC and SAR, respectively) using standard packages (see footnote 17 in the paper)
# Panels B and C T2:
# Other models using standard packages
# SAC
panel.sac.gmm <- function(Offense){
Data$y <- Data[[Offense]]
df <- Data %>%
select(uniqueid, year, y) %>%
drop_na
out <- summary(spgm(
y ~ factor(year),
data = df,
index = c("uniqueid", "year"),
lag = TRUE,
listw = mat2listw(W),
model = "within",
spatial.error = TRUE ))
Re <- createTexreg(
coef.names = c("Rho"),
coef = c(out$CoefTable[1,1]),
se = c(out$CoefTable[1,2]),
pvalues = c(out$CoefTable[1,4]),
gof = c(out$rho[1], out$rho[2]),
gof.names = c("Lambda", "sigma^2"),
gof.decimal = c(TRUE,TRUE),
model.name = c(substring(Offense, first=7))
)
return(screenreg(Re, stars = c(0.01, 0.05, 0.1), digits = 3, dcolumn = TRUE))
}
# SAR:
panel.sar.gmm <- function(Offense){
Data$y <- Data[[Offense]]
df <- Data %>%
select(uniqueid, year, y) %>%
drop_na
out <- summary(spgm(
y ~ factor(year),
data = df,
index = c("uniqueid", "year"),
lag = TRUE,
listw = mat2listw(W),
model = "within",
spatial.error = FALSE ))
Re <- createTexreg(
coef.names = c("Rho"),
coef = c(out$CoefTable[1,1]),
se = c(out$CoefTable[1,2]),
pvalues = c(out$CoefTable[1,4]),
gof = c(out$sigmav),
gof.names = c("sigma2"),
gof.decimal = c(TRUE),
model.name = c(substring(Offense, first=7))
)
return(screenreg(Re, stars = c(0.01, 0.05, 0.1), digits = 3, dcolumn = TRUE))
}
# Other results in Table 2:
# Coeff in Panel B: SAC model with FE
lapply(aCr, panel.sac.gmm)
## [[1]]
##
## ===================
## TheftHo
## -------------------
## Rho 0.855 ***
## (0.143)
## -------------------
## Lambda 0.181
## sigma^2 0.161
## ===================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[2]]
##
## ================
## PersInj
## ----------------
## Rho 0.186
## (0.133)
## ----------------
## Lambda 0.386
## sigma^2 0.157
## ================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[3]]
##
## ===================
## Homicid
## -------------------
## Rho 0.754 ***
## (0.202)
## -------------------
## Lambda 0.135
## sigma^2 0.126
## ===================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[4]]
##
## ===================
## TheftPe
## -------------------
## Rho -0.524 ***
## (0.169)
## -------------------
## Lambda 0.586
## sigma^2 0.100
## ===================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[5]]
##
## ===================
## Extorsi
## -------------------
## Rho 1.027 ***
## (0.124)
## -------------------
## Lambda -0.417
## sigma^2 0.086
## ===================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[6]]
##
## ===================
## Sexualc
## -------------------
## Rho 0.811 ***
## (0.148)
## -------------------
## Lambda 0.029
## sigma^2 0.173
## ===================
## *** p < 0.01; ** p < 0.05; * p < 0.1
# Other results in Table 2:
# Coeff in Panel C: SAR model with FE
lapply(aCr, panel.sar.gmm)
## [[1]]
##
## ==================
## TheftHo
## ------------------
## Rho 0.915 ***
## (0.112)
## ------------------
## sigma2 0.162
## ==================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[2]]
##
## ==================
## PersInj
## ------------------
## Rho 0.562 ***
## (0.102)
## ------------------
## sigma2 0.160
## ==================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[3]]
##
## ==================
## Homicid
## ------------------
## Rho 0.792 ***
## (0.164)
## ------------------
## sigma2 0.127
## ==================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[4]]
##
## ================
## TheftPe
## ----------------
## Rho 0.241 *
## (0.131)
## ----------------
## sigma2 0.106
## ================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[5]]
##
## ==================
## Extorsi
## ------------------
## Rho 0.933 ***
## (0.150)
## ------------------
## sigma2 0.086
## ==================
## *** p < 0.01; ** p < 0.05; * p < 0.1
##
## [[6]]
##
## ==================
## Sexualc
## ------------------
## Rho 0.819 ***
## (0.132)
## ------------------
## sigma2 0.173
## ==================
## *** p < 0.01; ** p < 0.05; * p < 0.1