### R code for "GMM Estimators for Binary Spatial Models in R" ################################################### ### code chunk number 3: boston-shape ################################################### library("sf") boston.tr <- st_read(system.file("shapes/boston_tracts.shp", package ="spData")[1], quiet = TRUE) ################################################### ### code chunk number 4: Wboston ################################################### library("spdep") boston_nb <- poly2nb(boston.tr) W <- nb2listw(boston_nb, style = "W") ################################################### ### code chunk number 5: simulated-data-part-1 ################################################### library("spatialreg") set.seed(1) n <- length(W$neighbours) lambda <- 0.6 beta0 <- -0.5 beta1 <- beta2 <- beta3 <- 1 A_i <- invIrM(boston_nb, lambda) x <- rnorm(n) z <- runif(n) Wx <- lag.listw(W, x) epsilon <- rnorm(n) ################################################### ### code chunk number 6: simulated-data-part-2 ################################################### ystar <- A_i %*% (beta0 + beta1 * x + beta2 * Wx + beta3 * z + epsilon) y <- as.numeric(ystar > 0) data <- as.data.frame(cbind(y, x, z)) ################################################### ### code chunk number 7: me-true-part-1 ################################################### Sigma_u <- tcrossprod(A_i) sigma <- sqrt(diag(Sigma_u)) index <- beta0 + beta1 * x + beta2 * Wx + beta3 * z a <- as.vector(A_i %*% index / sigma) P <- diag(dnorm(a)) %*% A_i %*% diag(1 / sigma) C_x <- P %*% (diag(n) * beta1 + beta2 * listw2mat(W)) C_z <- P %*% (diag(n) * beta3) ################################################### ### code chunk number 8: me-true-part-2 ################################################### TE_x <- (sum(C_x) / n) DE_x <- (sum(diag(C_x))/ n) IE_x <- TE_x - DE_x TE_z <- (sum(C_z) / n) DE_z <- (sum(diag(C_z))/ n) IE_z <- TE_z - DE_z Teffects <- rbind(cbind(TE_x, DE_x, IE_x), cbind(TE_z, DE_z, IE_z)) colnames(Teffects) <- c("ATE", "ADE", "AIE") rownames(Teffects) <- c("x", "z") Teffects ################################################### ### code chunk number 9: Figure ################################################### library("ggplot2") library("gridExtra") library("RColorBrewer") boston.tr$ystar <- ystar boston.tr$y <- as.factor(y) boston.tr$tex <- as.vector(t(rep(1, n)) %*% C_x) boston.tr$tez <- as.vector(t(rep(1, n)) %*% C_z) plot1 <- ggplot() + geom_sf(data = boston.tr, aes(fill = ystar)) + scale_fill_gradientn( colors = c("#9DBF9E", "#FCB97D", "#A84268"), name = "" ) + theme_void() + labs(title = expression(paste("(a) Latent ", bold(y)))) plot2 <- ggplot() + geom_sf(data = boston.tr, aes(fill = factor(y))) + scale_fill_manual( values = c("#9DBF9E", "#56B4E9"), name = "" ) + theme_void() + labs(title = expression(paste("(b) Observed ", bold(y)))) grid.arrange(plot1, plot2, ncol = 2) ################################################### ### code chunk number 10: load-package ################################################### library("spldv") ################################################### ### code chunk number 11: os_R ################################################### os_R <- sbinaryGMM(y ~ x + z | x, link = "probit", listw = W, nins = 2, data = data, type = "onestep", winitial = "optimal") ################################################### ### code chunk number 12: instruments ################################################### head(os_R$H) ################################################### ### code chunk number 13: summary-os-R ################################################### summary(os_R) ################################################### ### code chunk number 14: os_I ################################################### os_I <- sbinaryGMM(y ~ x + z | x, link = "probit", listw = W, nins = 2, data = data, type = "onestep", winitial = "identity") summary(os_I) ################################################### ### code chunk number 15: me-approx ################################################### set.seed(1) summary(impacts(os_R, type = "mc", R = 100, approximation = TRUE, pw = 6)) ################################################### ### code chunk number 16: me-approx-delta ################################################### summary(impacts(os_R, type = "delta", approximation = TRUE, pw = 6)) ################################################### ### code chunk number 17: me-approx-delta-nohet ################################################### summary(impacts(os_R, type = "delta", approximation = TRUE, pw = 6, het = FALSE)) ################################################### ### code chunk number 18: te-units-part1 ################################################### theta.hat <- coef(os_R) lambda.hat <- theta.hat["lambda"] beta.hat <- theta.hat["x"] wbeta.hat <- theta.hat["lag_x"] X <- os_R$X n <- nrow(X) W <- os_R$listw I <- Diagonal(n) A_i <- solve(I - lambda.hat * W) sigmas_u <- sqrt(diag(tcrossprod(A_i))) ################################################### ### code chunk number 19: te-units-part2 ################################################### D_i <- Diagonal(x = 1 / sigmas_u) a <- D_i %*% A_i %*% X %*% as.matrix(theta.hat[1:ncol(X)]) dfa <- Diagonal(x = dnorm(as.numeric(a))) Chat_x <- dfa %*% D_i %*% A_i %*% (I * beta.hat + wbeta.hat * W) TE_hat <- as.vector(t(rep(1, n)) %*% Chat_x) ################################################### ### code chunk number 20: Figure2 ################################################### Compare <- cbind(boston.tr$tex, TE_hat) plot(1:n, boston.tr$tex[order(boston.tr$tex)], type = "l", ylab = "Total Effects", xlab = "Census Tracts", main = "") lines(TE_hat[order(boston.tr$tex)], type = "l", col = "red") legend("topleft", c("True Total Effect", "Estimated Total Effect"), lty = c(1, 1), col = c("black", "red")) ################################################### ### code chunk number 21: ts-H ################################################### ts_H <- sbinaryGMM(y ~ x + z | x, link = "probit", listw = W, nins = 2, data = data, type = "twostep", winitial = "optimal") summary(ts_H) ################################################### ### code chunk number 22: ts-H-E ################################################### summary(ts_H, vce = "efficient") ################################################### ### code chunk number 23: me-ts ################################################### summary(impacts(ts_H, type = "delta", vce = "efficient", approximation = TRUE, pw = 6)) ################################################### ### code chunk number 24: linear ################################################### lgmm <- sbinaryLGMM(y ~ x + z | x, link = "probit", listw = W, nins = 2, data = data) summary(lgmm) ################################################### ### code chunk number 25: effect.linear ################################################### summary(impacts(lgmm, type = "delta", approximation = TRUE, pw = 6) ) ################################################### ### code chunk number 26: columbus ################################################### columbus <- st_read(system.file("shapes/columbus.shp", package = "spData")[1], quiet = TRUE) col.nb <- read.gal(system.file("weights/columbus.gal", package = "spData")[1]) W.col <- nb2listw(col.nb, style = "W") ################################################### ### code chunk number 27: binary_dependent ################################################### columbus$CRIMED <- as.numeric(columbus$CRIME > 37) ################################################### ### code chunk number 28: slm-models ################################################### slm <- CRIMED ~ INC + HOVAL probit <- glm(slm, family = binomial("probit"), data = columbus) lgmm <- sbinaryLGMM(slm, link = "probit", listw = W.col, data = columbus) osI <- sbinaryGMM(slm, link = "probit", listw = W.col, data = columbus, type = "onestep", winitial = "identity", cons.opt = TRUE, verbose = FALSE) osR <- sbinaryGMM(slm, link = "probit", listw = W.col, data = columbus, type = "onestep", winitial = "optimal", cons.opt = TRUE, verbose = FALSE) tsI <- sbinaryGMM(slm, link = "probit", listw = W.col, data = columbus, type = "twostep", winitial = "identity", cons.opt = TRUE, verbose = FALSE) tsR <- sbinaryGMM(slm, link = "probit", listw = W.col, data = columbus, type = "twostep", winitial = "optimal", cons.opt = TRUE, verbose = FALSE) tsR_E <- summary(tsR, vce = "efficient") ################################################### ### code chunk number 29: table (eval = FALSE) ################################################### ## library("memisc") ## table <- mtable("Probit" = probit, ## "LGMM" = lgmm, ## "OS-I" = osI, ## "OS-O" = osR, ## "TS-I" = tsI, ## "TS-O" = tsR, ## "TS-O-E" = tsR_E, ## summary.stats = c("N"), ## signif.symbols = c("***" = .01, "**" = 0.05, "*" = 0.1)) ## table ################################################### ### code chunk number 30: table1 ################################################### library("memisc") table <- mtable("Probit" = probit, "LGMM" = lgmm, "OS-I" = osI, "OS-O" = osR, "TS-I" = tsI, "TS-O" = tsR, "TS-O-E" = tsR_E, summary.stats = c("N"), signif.symbols = c("***" = .01, "**" = 0.05, "*" = 0.1)) toLatex(table, compact = TRUE, useBooktabs = TRUE) ################################################### ### code chunk number 31: impacts-columbus ################################################### me_delta <- impacts(tsR, type = "delta", vce = "efficient") summary(me_delta) ################################################### ### code chunk number 32: sdm ################################################### sdm <- CRIMED ~ INC + HOVAL | INC + HOVAL tsR_sdm <- sbinaryGMM(sdm, link = "probit", listw = W.col, data = columbus, type = "twostep", winitial = "optimal", cons.opt = TRUE, verbose = FALSE) summary(tsR_sdm, vce = "efficient") ################################################### ### code chunk number 33: sdm-test ################################################### library("car") coefs <- names(coef(tsR_sdm)) linearHypothesis(tsR_sdm, coefs[grep("lag_", coefs)], vcov = vcov(tsR_sdm, type = "efficient")) ################################################### ### code chunk number 34: comparison1 ################################################### library("McSpatial") data$wx <- as.matrix(W) %*% x sdm <- y ~ x + z + wx mcSlgmm <- spprobit(sdm, data = data, wmat = as.matrix(W), winst = ~ z + wx) lgmm2 <- sbinaryLGMM(sdm, link = "probit", listw = W, data = data, nins = 1) ################################################### ### code chunk number 35: comparison2 ################################################### probit <- glm(sdm, family = binomial("probit"), data = data) mcSgmm <- gmmprobit(sdm, data = data, wmat = as.matrix(W), startb = coef(probit), startrho = 0, winst = ~ z + wx) os_2 <- sbinaryGMM(sdm, link = "probit", listw = W, data = data, type = "onestep", winitial = "optimal", s.matrix = "iid", nins = 1, start = c(coef(probit), 0), reltol = 0.0001) ################################################### ### code chunk number 36: comparison3 ################################################### os_a <- sbinaryGMM(sdm, link = "probit", listw = W, data = data, type = "onestep", winitial = "optimal", nins = 1, approximation = TRUE, pw = 4, reltol = 0.0001) ################################################### ### code chunk number 37: comparison4 ################################################### library("spatialprobit") library("ProbitSpatial") bayes <- sarprobit(sdm, W = W, data = data, showProgress = FALSE) ml_cond <- ProbitSpatialFit(sdm, W = W, data = data, DGP = "SAR", method = "conditional", varcov = "varcov") ################################################### ### code chunk number 38: se-models ################################################### sbayes <- summary(bayes) summary(ml_cond, covar = TRUE) se_ml <- c(0.013307451, 0.017058855, 0.029112690, 0.036917704, 0.006378866) ################################################### ### code chunk number 39: comparison5 ################################################### library("INLA") n <- nrow(data) data$idx <- 1:n e <- eigen(W)$values re.idx <- which(abs(Im(e)) < 1e-6) rho.max <- 1 / max(Re(e[re.idx])) rho.min <- 1 / min(Re(e[re.idx])) mm <- model.matrix(sdm, data) betaprec1 <- 0.001 Q.beta1 <- Diagonal(n = ncol(mm), x = 1) Q.beta1 <- betaprec1 * Q.beta1 hyper.slm <- list( prec = list(initial = log(1), fixed = TRUE), theta2 = list(prior = "logitbeta", param = c(1, 1)) ) inlaR <- inla(y ~ -1 + f(idx, model="slm", args.slm = list(rho.min = rho.min, rho.max = rho.max, W = W, X = mm, Q.beta = Q.beta1), hyper = hyper.slm), data = data, family = "binomial", control.family = list(link = "probit"), control.compute = list(dic = TRUE, cpo = TRUE, config = TRUE) ) ################################################### ### code chunk number 40: lambda-function ################################################### ff <- function(z) { z * (rho.max - rho.min) + rho.min } lambda <- inla.tmarginal(ff, inlaR$marginals.hyperpar[[1]]) lambda.mean <- inla.zmarginal(lambda)[c(1, 2, 3, 7)]$mean lambda.se <- inla.zmarginal(lambda)[c(1, 2, 3, 7)]$sd ################################################### ### code chunk number 41: table2 ################################################### # Generate Table Table <- matrix(NA, nrow = 5, ncol = 16) Table[1:5, 1:2] <- summary(lgmm2)$Coef[, 1:2] Table[1:5, 3] <- mcSlgmm$coef Table[1:5, 4] <- mcSlgmm$se Table[1:5, 5:6] <- summary(os_2)$Coef[, 1:2] Table[1:5, 7] <- mcSgmm$coef Table[1:5, 8] <- mcSgmm$se Table[1:5, 9:10] <- summary(os_a)$Coef[, 1:2] Table[1:5, 11:12] <- sbayes[, 1:2] Table[1:5, 13] <- ml_cond@coeff Table[1:5, 14] <- se_ml Table[1:5, 15] <- c(inlaR$summary.random$idx[n + 1:(1 + 3), 2], lambda.mean) Table[1:5, 16] <- c(inlaR$summary.random$idx[n + 1:(1 + 3), 3], lambda.se) rownames(Table) <- c("Cons", "x", "z", "wx", "lambda") Table