#=================================== # # Main Regressions for: Do monetary subjective well-being evaluations vary across space? # Comparing continuous and discrete spatial heterogeneity # # By: Mauricio Sarrias # # Date: May, 22, 2014 # # Last Modification: June, 27, 2018 # #=================================== # ## Load Packages ---- library("Rchoice") library("foreign") library("memisc") library("stargazer") library("msm") library("lmtest") library("plotrix") ## Load and arrange data --- load("Res.RData") ### Some functions ===== ### function for CV cv.Rchoice <- function(object, wrt = NULL, digits = max(3, getOption("digits") - 2)){ if (is.null(wrt)) stop("CV needs the variable in the denominator: wrt") beta.hat <- coef(object) posi <- match(wrt, names(beta.hat)) form <- c() b <- c() namesb <- names(beta.hat)[-c(posi)] for (i in 1:length(beta.hat)) { if (i != posi) { b <- c(b, beta.hat[i] / beta.hat[posi]) form <- c(form, paste("~", "x", i, "/", "x", posi, sep = "")) } } names(b) <- namesb std.err <- c() for (i in 1:length(form)) { std.err <- c(std.err, msm::deltamethod(as.formula(form[i]), beta.hat, vcov(object), ses = TRUE)) } z <- b / std.err p <- 2 * (1 - pnorm(abs(z))) tablewtp <- cbind(b, std.err, z, p) colnames(tablewtp) <- c("Estimate", "Std. Error", "t-value", "Pr(>|t|)") cat(paste("\nCompesating variation respect to: ", wrt, "\n\n")) printCoefmat(tablewtp, digits = digits) } ## function for test ben_akiva <- function(mod1, mod2, forL){ # H0: mod1 is the correct specification L1 <- as.numeric(logLik(mod1)) L2 <- as.numeric(logLik(mod2)) K1 <- length(coef(mod1)) K2 <- length(coef(mod2)) theta <- c(coef(forL)["mean.constant"], rep(0, ncol(model.matrix(forL)) - 1)) L0 <- sum(Rchoice:::lnbinary(theta = theta, y = forL$mf$slife_binA, X = model.matrix(forL), link = "probit")) rho_1 <- 1 - ((L1 - K1) / L0) rho_2 <- 1 - ((L2 - K2) / L0) z <- rho_2 - rho_1 p <- pnorm(-sqrt( -2 * z * L0 + (K2 - K1) ), mean = 0, sd = 1) out <- list(z = z, p = p) return(out) } #=========================================# # I. Probit models (Fixed Parameters) ==== #=========================================# bin_A <- Rchoice(slife_binA ~ linch + edu + age + sex + dlstatus2 + dlstatus3 + dcivil1 + hsize + disability + medtreat + accident + noise + airpol + watpol + waspol + robprob + vigprob + alcdrogprob + trafdrogprob + ldens + lmdinc + urban, data = health, family = binomial('probit')) summary(bin_A) AIC(bin_A) BIC(bin_A) ### Compute CV for fixed parameter model cv.Rchoice(bin_A, wrt = "linch") #======================================================# # II. Probit models with Random Parameters ===== #======================================================# # All Normal and commune variables (Warning: it takes around 15 min) ==== ran_allnAcv <- Rchoice(slife_binA ~ linch + edu + age + sex + dlstatus2 + dlstatus3 + dcivil1 + hsize + disability + medtreat + accident + noise + airpol + watpol + waspol + robprob + vigprob + alcdrogprob + trafdrogprob + ldens + lmdinc + urban, data = health, family = binomial('probit'), ranp = c(constant = "n", edu = "n", age = "n", sex = "n", dlstatus2 = "n", dlstatus3 = "n", dcivil1 = "n", hsize = "n", disability = "n", medtreat = "n", accident = "n", noise = "n", airpol = "n", watpol = "n", waspol = "n", robprob = "n", vigprob = "n", alcdrogprob = "n", trafdrogprob = "n", ldens = "n", lmdinc = "n", urban = "n"), R = 100, panel = TRUE, index = "idc", print.level = 3, method = "bfgs") summary(ran_allnAcv) AIC(ran_allnAcv) BIC(ran_allnAcv) length(coef(ran_allnAcv)) ## Compensated variations cv.Rchoice(ran_allnAcv, wrt = "linch") ### Summary Statistics (or Table 1) stargazer(as.data.frame(cbind(health$slife_binA, model.matrix(ran_allnAcv)))) ## LR test lrtest(bin_A, ran_allnAcv) # Correlated parameters Warning: it takes time) ==== ran_CrvA <- Rchoice(slife_binA ~ linch + edu + age + sex + dlstatus2 + dlstatus3 + dcivil1 + hsize + disability + medtreat + accident + noise + airpol + watpol + waspol + robprob + vigprob + alcdrogprob + trafdrogprob + ldens + lmdinc + urban, data = health, family = binomial('probit'), ranp = c(constant = "n", edu = "n", age = "n", sex = "n", dlstatus2 = "n", dlstatus3 = "n", dcivil1 = "n", hsize = "n", disability = "n", medtreat = "n", accident = "n", noise = "n", airpol = "n", watpol = "n", waspol = "n", robprob = "n", vigprob = "n", alcdrogprob = "n", trafdrogprob = "n", ldens = "n", lmdinc = "n", urban = "n"), R = 100, panel = TRUE, index = "idc", print.level = 3, init.ran = 0.1, method = "bhhh", correlation = TRUE) summary(ran_CrvA) ## Compute standard deviation and ses vcov(ran_CrvA, what = "ranp", type = "sd", se = TRUE) AIC(ran_CrvA) BIC(ran_CrvA) # Truncated and Log-normal distribution ==== health$iage <- -1 * health$age health$idlstatus2 <- -1 * health$dlstatus2 health$imedtreat <- -1 * health$medtreat health$inoise <- -1 * health$noise health$idisability <- -1 * health$disability health$iaccident <- -1 * health$accident health$iairpol <- -1 * health$airpol health$iwatpol <- -1 * health$watpol health$ivispol <- -1 * health$vispol health$iwaspol <- -1 * health$waspol health$irobprob <- -1 * health$robprob health$ivigprob <- -1 * health$vigprob health$ialcdrogprob <- -1 * health$alcdrogprob health$itrafdrogprob <- -1 * health$trafdrogprob ran_tnCorAcv3 <- Rchoice(slife_binA ~ linch + edu + age + sex + idlstatus2 + dlstatus3 + dcivil1 + hsize + idisability + imedtreat + iaccident + inoise + iairpol + iwatpol + iwaspol + irobprob + ivigprob + ialcdrogprob + itrafdrogprob + ldens + lmdinc + urban, data = health, family = binomial('probit'), ranp = c(constant = "n", edu = "ln", age = "n", sex = "n", idlstatus2 = "ln", dlstatus3 = "n", dcivil1 = "n", hsize = "n", idisability = "ln", imedtreat = "ln", iaccident = "ln", inoise = "cn", iairpol = "cn", iwatpol = "cn", iwaspol = "cn", irobprob = "ln", ivigprob = "ln", ialcdrogprob = "ln", itrafdrogprob = "ln", ldens = "n", lmdinc = "n", urban = "n"), R = 100, panel = TRUE, index = "idc", print.level = 3, correlation = TRUE, method = "bhhh", init.ran = 0) summary(ran_tnCorAcv3) vcov(ran_tnCorAcv3, what = "ranp", type = "sd", se = TRUE) AIC(ran_tnCorAcv3) BIC(ran_tnCorAcv3) length(coef(ran_tnCorAcv3)) # Function to compute SE of complex standard errors complex_se <- function(x, var, deno){ all_coef <- coef(x) n_bbar <- paste("mean", var, sep = ".") b_bar <- all_coef[n_bbar] pos_bb <- which(names(all_coef) == n_bbar) Ka <- length(x$ranp) nr <- names(x$ranp) names.stds <- c() for (i in 1:Ka) names.stds <- c(names.stds, paste('sd', nr[i], nr[i:Ka], sep = '.')) stds.hat <- all_coef[names.stds] i <- which(names(x$ranp) == var) j <- 1 wf <- c() while (j <= i) { wf <- c(wf, paste("sd", nr[j], nr[i], sep = ".")) j <- j + 1 } pos_wff <- c() for (i in 1:length(wf)) pos_wff <- c(pos_wff,which(names(all_coef) == wf[i])) x.var <- paste(paste("x", pos_wff, "^2", sep = "")) x.var.s <- paste(x.var, collapse = '+') sd_hat <- sqrt(sum(all_coef[wf]^2)) # Now we compute the mean of log and se # mean of log-normal mean <- exp(b_bar + (abs(sd_hat)^2)/2) form_mean <- paste("~exp(x", pos_bb, "+", "(", x.var.s, ")", "/2", ")", sep = "") mean.se <- msm::deltamethod(as.formula(form_mean), mean = coef(x), cov = vcov(x), ses = TRUE) se <- mean * sqrt(exp(sd_hat^2) - 1) form_se <- paste("~exp(x", pos_bb, "+", "(", x.var.s, ")", "/2", ")", "* sqrt(exp(", x.var.s,") - 1)", sep = "") se.se <- msm::deltamethod(as.formula(form_se), mean = coef(x), cov = vcov(x), ses = TRUE) pos_deno <- which(names(coef(x)) == deno) b_deno <- coef(x)[pos_deno] cv <- exp(b_bar + (abs(sd_hat)^2)/2) / b_deno form_cv <- paste("~exp(x", pos_bb, "+", "(", x.var.s, ")", "/2", ")", "* sqrt(exp(", x.var.s,") - 1) / ", paste("x", pos_deno, sep = ""), sep = "") cv.se <- msm::deltamethod(as.formula(form_cv), mean = coef(x), cov = vcov(x), ses = TRUE) cat("mean of ", var, "mean = ", mean, "se = ", paste("se = ", paste0('"=', round(mean.se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(mean / mean.se))), "\n") cat("sd of ", var, "mean = ", se, "se = ", paste("se = ", paste0('"=', round(se.se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(se / se.se))), "\n") cat("cv ", var, "mean = ", cv, "se = ", paste("se = ", paste0('"=', round(cv.se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(cv / cv.se))), "\n") } complex_se(ran_tnCorAcv3, "edu", "linch") complex_se(ran_tnCorAcv3, "idlstatus2", "linch") complex_se(ran_tnCorAcv3, "idisability", "linch") complex_se(ran_tnCorAcv3, "imedtreat", "linch") complex_se(ran_tnCorAcv3, "iaccident", "linch") complex_se(ran_tnCorAcv3, "irobprob", "linch") complex_se(ran_tnCorAcv3, "ivigprob", "linch") complex_se(ran_tnCorAcv3, "ialcdrogprob", "linch") complex_se(ran_tnCorAcv3, "itrafdrogprob", "linch") #SE for the rest of CVs names <- names(coef(ran_tnCorAcv3))[-c(1, 2)] names for (i in 1:length(names)) { mean <- -coef(ran_tnCorAcv3)[names[i]] / coef(ran_tnCorAcv3)["linch"] form <- paste(paste("~ - x", i + 2, sep = ""), paste("/ x1")) se <- deltamethod(as.formula(form), mean = coef(ran_tnCorAcv3), cov = vcov(ran_tnCorAcv3), ses = TRUE) cat("for variable", names[i], "cv = ", mean, "z = ", mean/se, paste("se = ", paste0('"', round(se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(mean/se))), "\n") } #======================================================# ### III. Latent Class Models: Q = 2 ==== #======================================================# Xm <- model.matrix(bin_A) ym <- health$slife_binA source("max_bin_lc.R") id <- health$idc H <- cbind(health$x_stub/1000000, health$y_stub/1000000) H <- H[!duplicated(id), ] colnames(H) <- c("u", "v") # Initial values Q <- 2 mean <- coef(bin_A) init.shift <- seq(-0.02, 0.02, length.out = Q) lc.mean <- c() for (i in 1:Q) lc.mean <- c(lc.mean, mean + init.shift[i]) lc.names <- c() for (i in 1:Q) lc.names <- c(lc.names, paste('class', i, names(mean), sep = '.')) names(lc.mean) <- lc.names # H HE <- H[rep(seq_len(nrow(H)), each = Q), ] class <- factor(rep(1:Q, length(unique(id)))) cldata <- as.data.frame(cbind(HE, factor(class), row.names = NULL)) hv <- model.matrix(formula(~ u:class+ v:class, lhs = 0), cldata)[, -1, drop = FALSE] colnames(hv) lev1 <- levels(class)[1] lev1 <- paste("class", lev1, sep = "") toremove <- unlist(lapply(as.list(c("u", "v")), function(x) paste(lev1, x, sep = ":"))) revtoremove <- unlist(lapply(as.list(c("u", "v")), function(x) paste(x, lev1, sep = ":"))) toremove <- colnames(hv) %in% c(toremove, revtoremove) hv <- hv[, !toremove, drop = FALSE] names.hv <- colnames(hv) Hl <- vector(length = Q, mode = "list") names(Hl) <- levels(class) for (i in levels(class)) { Hl[[i]] <- hv[class == i, , drop = FALSE] } # Initial valeus s_class <- rep(0, ncol(hv)) names(s_class) <- colnames(hv) start <- c(lc.mean, s_class) library("maxLik") source("max_bin_lc.R") rep.row <- function(x, n){ matrix(rep(x,each = n),nrow = n) } out_2 <- maxLik(ml_bin_lc, start = start, method = "bhhh", y = ym, X = Xm, H = Hl, Q = Q, id = id, weights = 1, gradient = TRUE, get.bi = FALSE, iterlim = 5000, print.level = 3) summary(out_2) b.hat_A2 <- coef(out_2) out.hat_A2 <- ml_bin_lc(b.hat_A2, y = ym, X = Xm, H = Hl, id = id, weights = 1, gradient = FALSE, get.bi = TRUE, Q = 2) bi_A2 <- attr(out.hat_A2, "bi") Qnq_A2 <- attr(out.hat_A2, "Qir") colMeans(attr(out.hat_A2, "Wnq")) K <- ncol(bi_A2) N <- length(unique(id)) mean_A2 <- mean.sq_A2 <- array(NA, dim = c(N, Q = 2, K)) for (j in 1:K) { mean_A2[, , j] <- repRows(bi_A2[, j] / bi_A2[, "linch"], N) * Qnq_A2 mean.sq_A2[, , j] <- (repRows(bi_A2[, j] / bi_A2[, "linch"], N) ^ 2) * Qnq_A2 } mean_A2 <- apply(mean_A2, c(1,3), sum) mean.sq_A2 <- apply(mean.sq_A2, c(1,3), sum) sd.est_A2 <- suppressWarnings(sqrt(mean.sq_A2 - mean_A2 ^ 2)) colnames(mean_A2) <- colnames(mean.sq_A2) <- colnames(sd.est_A2) <- colnames(bi_A2) #se class 1 names <- names(coef(out_2)[1:23]) names names_s <- names[-c(1, 2)] for (i in 1:length(names_s)) { mean <- -coef(out_2)[names_s[i]] / coef(out_2)["class.1.linch"] form <- paste(paste("~ - x", i + 2, sep = ""), paste("/ x2")) se <- deltamethod(as.formula(form), mean = coef(out_2), cov = vcov(out_2), ses = TRUE) cat("for variable", names_s[i], "cv = ", mean, "z = ", mean / se, paste("se = ", paste0('"', round(se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(mean/se))), "\n") } #se class 2 names <- names(coef(out_2)[24:46]) names names_s <- names[-c(1, 2)] for (i in 1:length(names_s)) { mean <- -coef(out_2)[names_s[i]] / coef(out_2)["class.2.linch"] form <- paste(paste("~ - x", 23 + i + 2, sep = ""), paste("/ x25")) se <- deltamethod(as.formula(form), mean = coef(out_2), cov = vcov(out_2), ses = TRUE) cat("for variable", names_s[i], "cv = ", mean, "z = ", mean / se, paste("se = ", paste0('"', round(se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(mean/se))), "\n") } #======================================================# ### IV. Latent Class Models: Q = 3 ==== #======================================================# # Initial values Q <- 3 mean <- coef(bin_A) init.shift <- seq(-0.002, 0.002, length.out = Q) lc.mean <- c() for (i in 1:Q) lc.mean <- c(lc.mean, mean + init.shift[i]) lc.names <- c() for (i in 1:Q) lc.names <- c(lc.names, paste('class', i, names(mean), sep = '.')) names(lc.mean) <- lc.names # H HE <- H[rep(seq_len(nrow(H)), each = Q), ] class <- factor(rep(1:Q, length(unique(id)))) cldata <- as.data.frame(cbind(HE, factor(class), row.names = NULL)) hv <- model.matrix(formula(~ u:class+ v:class, lhs = 0), cldata)[, -1, drop = FALSE] colnames(hv) lev1 <- levels(class)[1] lev1 <- paste("class", lev1, sep = "") toremove <- unlist(lapply(as.list(c("u", "v")), function(x) paste(lev1, x, sep = ":"))) revtoremove <- unlist(lapply(as.list(c("u", "v")), function(x) paste(x, lev1, sep = ":"))) toremove <- colnames(hv) %in% c(toremove, revtoremove) hv <- hv[, !toremove, drop = FALSE] names.hv <- colnames(hv) Hl <- vector(length = Q, mode = "list") names(Hl) <- levels(class) for (i in levels(class)) { Hl[[i]] <- hv[class == i, , drop = FALSE] } # Initial valeus s_class <- rep(0, ncol(hv)) names(s_class) <- colnames(hv) start <- c(lc.mean, s_class) out_3 <- maxLik(ml_bin_lc, start = start, method = "bhhh", y = ym, X = Xm, H = Hl, Q = Q, id = id, weights = 1, gradient = TRUE, get.bi = FALSE, iterlim = 5000, print.level = 3) summary(out_3) b.hat_A3 <- coef(out_3) out.hat_A3 <- ml_bin_lc(b.hat_A3, y = ym, X = Xm, H = Hl, id = id, weights = 1, gradient = FALSE, get.bi = TRUE, Q = 3) bi_A3 <- attr(out.hat_A3, "bi") Qnq_A3 <- attr(out.hat_A3, "Qir") colMeans(attr(out.hat_A3, "Wnq")) K <- ncol(bi_A3) N <- length(unique(id)) mean_A3 <- mean.sq_A3 <- array(NA, dim = c(N, Q, K)) for (j in 1:K) { mean_A3[, , j] <- repRows(bi_A3[, j] / bi_A3[, "linch"], N) * Qnq_A3 mean.sq_A3[, , j] <- (repRows(bi_A3[, j] / bi_A3[, "linch"], N) ^2) * Qnq_A3 } mean_A3 <- apply(mean_A3, c(1,3), sum) mean.sq_A3 <- apply(mean.sq_A3, c(1,3), sum) sd.est_A3 <- suppressWarnings(sqrt(mean.sq_A3 - mean_A3 ^ 2)) colnames(mean_A3) <- colnames(mean.sq_A3) <- colnames(sd.est_A3) <- colnames(bi_A3) ### SE names <- names(coef(out_3)[1:23]) names_s <- names[-c(1, 2)] for (i in 1:length(names_s)) { mean <- -coef(out_3)[names_s[i]] / coef(out_3)["class.1.linch"] form <- paste(paste("~ - x", i + 2, sep = ""), paste("/ x2")) se <- deltamethod(as.formula(form), mean = coef(out_3), cov = vcov(out_3), ses = TRUE) cat("for variable", names_s[i], "cv = ", mean, "z = ", mean / se, paste("se = ", paste0('"', round(se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(mean/se))), "\n") } #se class 2 names <- names(coef(out_3)[24:46]) names_s <- names[-c(1, 2)] for (i in 1:length(names_s)) { mean <- -coef(out_3)[names_s[i]] / coef(out_3)["class.2.linch"] form <- paste(paste("~ - x", 23 + i + 2, sep = ""), paste("/ x25")) se <- deltamethod(as.formula(form), mean = coef(out_3), cov = vcov(out_3), ses = TRUE) cat("for variable", names_s[i], "cv = ", mean, "z = ", mean / se, paste("se = ", paste0('"', round(se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(mean/se))), "\n") } #se class 3 names <- names(coef(out_3)[47:69]) names_s <- names[-c(1, 2)] for (i in 1:length(names_s)) { mean <- -coef(out_3)[names_s[i]] / coef(out_3)["class.3.linch"] form <- paste(paste("~ - x", 46 + i + 2, sep = ""), paste("/ x48")) se <- deltamethod(as.formula(form), mean = coef(out_3), cov = vcov(out_3), ses = TRUE) cat("for variable", names_s[i], "cv = ", mean, "z = ", mean/se, paste("se = ", paste0('"', round(se, 3)),'"') , "p-value", 2 * (1 - pnorm(abs(mean/se))), "\n") } ### Test ==== ben_akiva(ran_allnAcv, out_2, forL = ran_allnAcv) ben_akiva(out_2, ran_CrvA, forL = ran_CrvA) ben_akiva(out_2,ran_tnCorAcv3, forL = ran_tnCorAcv3) #======================================================# ### V. Plots ==== #======================================================# m_edu_c1 <- -effect.Rchoice(ran_allnAcv, par = "edu", effect = "cv", wrt = "linch", ind = TRUE)$mean m_edu_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "edu", effect = "cv", wrt = "linch", ind = TRUE)$mean m_edu_q2 <- -mean_A2[, "edu"] m_edu_q3 <- -mean_A3[, "edu"] m_age_c1 <- -effect.Rchoice(ran_allnAcv, par = "age", effect = "cv", wrt = "linch", ind = TRUE)$mean m_age_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "age", effect = "cv", wrt = "linch", ind = TRUE)$mean m_age_q2 <- -mean_A2[, "age"] m_age_q3 <- -mean_A3[, "age"] m_sex_c1 <- -effect.Rchoice(ran_allnAcv, par = "sex", effect = "cv", wrt = "linch", ind = TRUE)$mean m_sex_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "sex", effect = "cv", wrt = "linch", ind = TRUE)$mean m_sex_q2 <- -mean_A2[, "sex"] m_sex_q3 <- -mean_A3[, "sex"] m_unem_c1 <- -effect.Rchoice(ran_allnAcv, par = "dlstatus2", effect = "cv", wrt = "linch", ind = TRUE)$mean s_unem_c1 <- effect.Rchoice(ran_allnAcv, par = "dlstatus2", effect = "cv", wrt = "linch", ind = TRUE)$sd.est lb_unem_c1 <- m_unem_c1 - 1.96 * s_unem_c1 ub_unem_c1 <- m_unem_c1 + 1.96 * s_unem_c1 m_unem_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "idlstatus2", effect = "cv", wrt = "linch", ind = TRUE)$mean s_unem_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "idlstatus2", effect = "cv", wrt = "linch", ind = TRUE)$sd.est lb_unem_c2 <- m_unem_c2 - 1.96 * s_unem_c2 ub_unem_c2 <- m_unem_c2 + 1.96 * s_unem_c2 m_unem_q2 <- -mean_A2[, "dlstatus2"] s_unem_q2 <- sd.est_A2[, "dlstatus2"] lb_unem_q2 <- m_unem_q2 - 1.96 * s_unem_q2 ub_unem_q2 <- m_unem_q2 + 1.96 * s_unem_q2 m_unem_q3 <- -mean_A3[, "dlstatus2"] s_unem_q3 <- sd.est_A3[, "dlstatus2"] lb_unem_q3 <- m_unem_q3 - 1.96 * s_unem_q3 ub_unem_q3 <- m_unem_q3 + 1.96 * s_unem_q3 sum((lb_unem_c1 < 0) & (0 < ub_unem_c1)) / length(m_unem_c1) sum((lb_unem_c2 < 0) & (0 < ub_unem_c2)) / length(m_unem_c2) sum((lb_unem_q2 < 0) & (0 < ub_unem_q2)) / length(m_unem_q2) 1 - sum((lb_unem_q3 < 0) & (0 < ub_unem_q3)) / length(m_unem_q3) library("msm") mean_un <- -coef(bin_A)["dlstatus2"] / coef(bin_A)["linch"] se_un <- deltamethod( ~ x16 / x2, mean = coef(bin_A), cov = vcov(bin_A), ses = TRUE) up_un <- mean_un + 1.96 * se_un lw_un <- mean_un - 1.96 * se_un par(mfrow = c(2, 2), mar = c(1, 1, 3, 1), oma = c(4, 4, 0, 0)) plotCI(1:length(m_unem_c1), m_unem_c1[order(m_unem_c1)], ui = ub_unem_c1[order(m_unem_c1)], li = lb_unem_c1[order(m_unem_c1)], lty = 2, pch = 19, xlab = "", ylab = expression(hat(E)(beta[c])), scol = "gray", col = "black", gap = TRUE, cex = 0.8, main = "A: Normal", ylim = c(-3, 4)) abline(h = 0, col = "red") abline(h = mean_un, col = "blue", lty = 2, lwd = 2) abline(h = up_un, col = "blue", lty = 3, lwd = 2) abline(h = lw_un, col = "blue", lty = 3, lwd = 2) plotCI(1:length(m_unem_c2), m_unem_c2[order(m_unem_c2)], ui = ub_unem_c2[order(m_unem_c2)], li = lb_unem_c2[order(m_unem_c2)], lty = 2, pch = 19, xlab = "", ylab = expression(hat(E)(beta[c])), scol = "gray", col = "black", gap = TRUE, cex = 0.8, main = "B: Log-Normal", ylim = c(-3, 4)) abline(h = 0, col = "red") abline(h = mean_un, col = "blue", lty = 2, lwd = 2) abline(h = up_un, col = "blue", lty = 3, lwd = 2) abline(h = lw_un, col = "blue", lty = 3, lwd = 2) plotCI(1:length(m_unem_q2), m_unem_q2[order(m_unem_q2)], ui = ub_unem_q2[order(m_unem_q2)], li = lb_unem_q2[order(m_unem_q2)], lty = 2, pch = 19, xlab = "", ylab = expression(hat(E)(beta[c])), scol = "gray", col = "black", gap = TRUE, cex = 0.8, main = "C:Dicrete Q = 2", ylim = c(-0.5, 1.5)) abline(h = 0, col = "red") abline(h = mean_un, col = "blue", lty = 2, lwd = 2) abline(h = up_un, col = "blue", lty = 3, lwd = 2) abline(h = lw_un, col = "blue", lty = 3, lwd = 2) plotCI(1:length(m_unem_q3), m_unem_q3[order(m_unem_q3)], ui = ub_unem_q3[order(m_unem_q3)], li = lb_unem_q3[order(m_unem_q3)], lty = 2, pch = 19, xlab = "", ylab = expression(hat(E)(beta[c])), scol = "gray", col = "black", gap = TRUE, cex = 0.8, main = "Discrete Q = 3", ylim = c(-0.5, 1.5)) abline(h = 0, col = "red") abline(h = mean_un, col = "blue", lty = 2, lwd = 2) abline(h = up_un, col = "blue", lty = 3, lwd = 2) abline(h = lw_un, col = "blue", lty = 3, lwd = 2) dev.off() m_olf_c1 <- -effect.Rchoice(ran_allnAcv, par = "dlstatus3", effect = "cv", wrt = "linch", ind = TRUE)$mean m_olf_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "dlstatus3", effect = "cv", wrt = "linch", ind = TRUE)$mean m_olf_q2 <- -mean_A2[, "dlstatus3"] m_olf_q3 <- -mean_A3[, "dlstatus3"] m_marr_c1 <- -effect.Rchoice(ran_allnAcv, par = "dcivil1", effect = "cv", wrt = "linch", ind = TRUE)$mean m_marr_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "dcivil1", effect = "cv", wrt = "linch", ind = TRUE)$mean m_marr_q2 <- -mean_A2[, "dcivil1"] m_marr_q3 <- -mean_A3[, "dcivil1"] m_hsize_c1 <- -effect.Rchoice(ran_allnAcv, par = "hsize", effect = "cv", wrt = "linch", ind = TRUE)$mean m_hsize_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "hsize", effect = "cv", wrt = "linch", ind = TRUE)$mean m_hsize_q2 <- -mean_A2[, "hsize"] m_hsize_q3 <- -mean_A3[, "hsize"] m_dis_c1 <- -effect.Rchoice(ran_allnAcv, par = "disability", effect = "cv", wrt = "linch", ind = TRUE)$mean s_dis_c1 <- effect.Rchoice(ran_allnAcv, par = "disability", effect = "cv", wrt = "linch", ind = TRUE)$sd.est lb_dis_c1 <- m_dis_c1 - 1.96 * s_dis_c1 ub_dis_c1 <- m_dis_c1 + 1.96 * s_dis_c1 m_dis_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "idisability", effect = "cv", wrt = "linch", ind = TRUE)$mean s_dis_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "idisability", effect = "cv", wrt = "linch", ind = TRUE)$sd.est lb_dis_c2 <- m_dis_c2 - 1.96 * s_dis_c2 ub_dis_c2 <- m_dis_c2 + 1.96 * s_dis_c2 m_dis_q2 <- -mean_A2[, "disability"] s_dis_q2 <- sd.est_A2[, "disability"] lb_dis_q2 <- m_dis_q2 - 1.96 * s_dis_q2 ub_dis_q2 <- m_dis_q2 + 1.96 * s_dis_q2 m_dis_q3 <- -mean_A3[, "disability"] s_dis_q3 <- sd.est_A3[, "disability"] lb_dis_q3 <- m_dis_q3 - 1.96 * s_dis_q3 ub_dis_q3 <- m_dis_q3 + 1.96 * s_dis_q3 # Plot for paper library("msm") mean_dis <- -coef(bin_A)["disability"] / coef(bin_A)["linch"] se_dis <- deltamethod( ~ x10 / x2, mean = coef(bin_A), cov = vcov(bin_A), ses = TRUE) up_dis <- mean_dis + 1.96 * se_dis lw_dis <- mean_dis - 1.96 * se_dis par(mfrow = c(2, 2), mar = c(1, 1, 3, 1), oma = c(4, 4, 0, 0)) plotCI(1:length(m_dis_c1), m_dis_c1[order(m_dis_c1)], ui = ub_dis_c1[order(m_dis_c1)], li = lb_dis_c1[order(m_dis_c1)], lty = 2, pch = 19, xlab = "", ylab = expression(hat(E)(beta[c])), scol = "gray", col = "black", gap = TRUE, cex = 0.8, main = "A: Normal", ylim = c(-3, 4)) abline(h = 0, col = "red") abline(h = mean_dis, col = "blue", lty = 2, lwd = 2) abline(h = up_dis, col = "blue", lty = 3, lwd = 2) abline(h = lw_dis, col = "blue", lty = 3, lwd = 2) plotCI(1:length(m_dis_c2), m_dis_c2[order(m_dis_c2)], ui = ub_dis_c2[order(m_dis_c2)], li = lb_dis_c2[order(m_dis_c2)], lty = 2, pch = 19, xlab = "", ylab = expression(hat(E)(beta[c])), scol = "gray", col = "black", gap = TRUE, cex = 0.8, main = "B: Log-Normal", ylim = c(-3, 4)) abline(h = 0, col = "red") abline(h = mean_dis, col = "blue", lty = 2, lwd = 2) abline(h = up_dis, col = "blue", lty = 3, lwd = 2) abline(h = lw_dis, col = "blue", lty = 3, lwd = 2) plotCI(1:length(m_dis_q2), m_dis_q2[order(m_dis_q2)], ui = ub_dis_q2[order(m_dis_q2)], li = lb_dis_q2[order(m_dis_q2)], lty = 2, pch = 19, xlab = "", ylab = expression(hat(E)(beta[c])), scol = "gray", col = "black", gap = TRUE, cex = 0.8, main = "C: Discrete Q = 2", ylim = c(-3, 4)) abline(h = 0, col = "red") abline(h = mean_dis, col = "blue", lty = 2, lwd = 2) abline(h = up_dis, col = "blue", lty = 3, lwd = 2) abline(h = lw_dis, col = "blue", lty = 3, lwd = 2) plotCI(1:length(m_dis_q3), m_dis_q3[order(m_dis_q3)], ui = ub_dis_q3[order(m_dis_q3)], li = lb_dis_q3[order(m_dis_q3)], lty = 2, pch = 19, xlab = "", ylab = expression(hat(E)(beta[c])), scol = "gray", col = "black", gap = TRUE, cex = 0.8, main = "D: Discrete Q = 3", ylim = c(-3, 4)) abline(h = 0, col = "red") abline(h = mean_dis, col = "blue", lty = 2, lwd = 2) abline(h = up_dis, col = "blue", lty = 3, lwd = 2) abline(h = lw_dis, col = "blue", lty = 3, lwd = 2) dev.off() m_med_c1 <- -effect.Rchoice(ran_allnAcv, par = "medtreat", effect = "cv", wrt = "linch", ind = TRUE)$mean m_med_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "imedtreat", effect = "cv", wrt = "linch", ind = TRUE)$mean m_med_q2 <- -mean_A2[, "medtreat"] m_med_q3 <- -mean_A3[, "medtreat"] m_acc_c1 <- -effect.Rchoice(ran_allnAcv, par = "accident", effect = "cv", wrt = "linch", ind = TRUE)$mean m_acc_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "iaccident", effect = "cv", wrt = "linch", ind = TRUE)$mean m_acc_q2 <- -mean_A2[, "accident"] m_acc_q3 <- -mean_A3[, "accident"] m_noise_c1 <- -effect.Rchoice(ran_allnAcv, par = "noise", effect = "cv", wrt = "linch", ind = TRUE)$mean m_noise_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "inoise", effect = "cv", wrt = "linch", ind = TRUE)$mean m_noise_q2 <- -mean_A2[, "noise"] m_noise_q3 <- -mean_A3[, "noise"] m_watpol_c1 <- -effect.Rchoice(ran_allnAcv, par = "watpol", effect = "cv", wrt = "linch", ind = TRUE)$mean m_watpol_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "iwatpol", effect = "cv", wrt = "linch", ind = TRUE)$mean m_watpol_q2 <- -mean_A2[, "watpol"] m_watpol_q3 <- -mean_A3[, "watpol"] m_waspol_c1 <- -effect.Rchoice(ran_allnAcv, par = "waspol", effect = "cv", wrt = "linch", ind = TRUE)$mean m_waspol_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "iwaspol", effect = "cv", wrt = "linch", ind = TRUE)$mean m_waspol_q2 <- -mean_A2[, "waspol"] m_waspol_q3 <- -mean_A3[, "waspol"] m_robprob_c1 <- -effect.Rchoice(ran_allnAcv, par = "robprob", effect = "cv", wrt = "linch", ind = TRUE)$mean m_robprob_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "irobprob", effect = "cv", wrt = "linch", ind = TRUE)$mean m_robprob_q2 <- -mean_A2[, "robprob"] m_robprob_q3 <- -mean_A3[, "robprob"] m_vigprob_c1 <- -effect.Rchoice(ran_allnAcv, par = "vigprob", effect = "cv", wrt = "linch", ind = TRUE)$mean m_vigprob_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "ivigprob", effect = "cv", wrt = "linch", ind = TRUE)$mean m_vigprob_q2 <- -mean_A2[, "vigprob"] m_vigprob_q3 <- -mean_A3[, "vigprob"] m_alcdrogprob_c1 <- -effect.Rchoice(ran_allnAcv, par = "alcdrogprob", effect = "cv", wrt = "linch", ind = TRUE)$mean m_alcdrogprob_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "ialcdrogprob", effect = "cv", wrt = "linch", ind = TRUE)$mean m_alcdrogprob_q2 <- -mean_A2[, "alcdrogprob"] m_alcdrogprob_q3 <- -mean_A3[, "alcdrogprob"] m_trafdrogprob_c1 <- -effect.Rchoice(ran_allnAcv, par = "trafdrogprob", effect = "cv", wrt = "linch", ind = TRUE)$mean m_trafdrogprob_c2 <- effect.Rchoice(ran_tnCorAcv3, par = "itrafdrogprob", effect = "cv", wrt = "linch", ind = TRUE)$mean m_trafdrogprob_q2 <- -mean_A2[, "trafdrogprob"] m_trafdrogprob_q3 <- -mean_A3[, "trafdrogprob"] m_ldens_c1 <- -effect.Rchoice(ran_allnAcv, par = "ldens", effect = "cv", wrt = "linch", ind = TRUE)$mean m_ldens_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "ldens", effect = "cv", wrt = "linch", ind = TRUE)$mean m_ldens_q2 <- -mean_A2[, "ldens"] m_ldens_q3 <- -mean_A3[, "ldens"] m_lmdinc_c1 <- -effect.Rchoice(ran_allnAcv, par = "lmdinc", effect = "cv", wrt = "linch", ind = TRUE)$mean m_lmdinc_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "lmdinc", effect = "cv", wrt = "linch", ind = TRUE)$mean m_lmdinc_q2 <- -mean_A2[, "lmdinc"] m_lmdinc_q3 <- -mean_A3[, "lmdinc"] m_urban_c1 <- -effect.Rchoice(ran_allnAcv, par = "urban", effect = "cv", wrt = "linch", ind = TRUE)$mean m_urban_c2 <- -effect.Rchoice(ran_tnCorAcv3, par = "urban", effect = "cv", wrt = "linch", ind = TRUE)$mean m_urban_q2 <- -mean_A2[, "urban"] m_urban_q3 <- -mean_A3[, "urban"] #Histograms ==== par(mfrow = c(3, 2), mar = c(1, 1, 3, 1), oma = c(4, 4, 0, 0)) boxplot(as.data.frame(cbind(m_edu_c1, m_edu_c2, m_edu_q2, m_edu_q3)), main = "A: Schooling", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_unem_c1, m_unem_c2, m_unem_q2, m_unem_q3)), main = "B: Unemployed", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_dis_c1, m_dis_c2, m_dis_q2, m_dis_q3)), main = "C: Disability", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_acc_c1, m_acc_c2, m_acc_q2, m_acc_q3)), main = "D: Accident", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_waspol_c1, m_waspol_c2, m_waspol_q2, m_waspol_q3)), main = "E: Garbage Pollution", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_trafdrogprob_c1, m_trafdrogprob_c2, m_trafdrogprob_q2, m_trafdrogprob_q3)), main = "F: Drug trafficking", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") dev.off() par(mfrow = c(4, 4), mar = c(1, 1, 3, 1), oma = c(4, 4, 0, 0)) boxplot(as.data.frame(cbind(m_age_c1, m_age_c2, m_age_q2, m_age_q3)), main = "A: Age", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_sex_c1, m_sex_c2, m_sex_q2, m_sex_q3)), main = "B: Male", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_olf_c1, m_olf_c2, m_olf_q2, m_olf_q3)), main = "C: OLF", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_marr_c1, m_marr_c2, m_marr_q2, m_marr_q3)), main = "D: Married", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_hsize_c1, m_hsize_c2, m_hsize_q2, m_hsize_q3)), main = "E: Household size", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_med_c1, m_med_c2, m_med_q2, m_med_q3)), main = "F: Medical Treatment", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_noise_c1, m_noise_c2, m_noise_q2, m_noise_q3)), main = "A: Noise Pollution", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_watpol_c1, m_watpol_c2, m_watpol_q2, m_watpol_q3)), main = "G: Water Pollution", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_robprob_c1, m_robprob_c2, m_robprob_q2, m_robprob_q3)), main = "H: Robbery", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_vigprob_c1, m_vigprob_c2, m_vigprob_q2, m_vigprob_q3)), main = "I: Poor Surveillance", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_alcdrogprob_c1, m_alcdrogprob_c2, m_alcdrogprob_q2, m_alcdrogprob_q3)), main = "J: Street drugs", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_ldens_c1, m_ldens_c2, m_ldens_q2, m_ldens_q3)), main = "K: Log Density", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_lmdinc_c1, m_lmdinc_c2, m_lmdinc_q2, m_lmdinc_q3)), main = "L: Log Median Income", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") boxplot(as.data.frame(cbind(m_urban_c1, m_urban_c2, m_urban_q2, m_urban_q3)), main = "M: Urban", names = c("C1", "C2", "D1", "D2")) abline(h = 0, col = "red") dev.off() #### End =====