A Replication Code for “Buying Votes across Borders?: A List Experiment on Mexican Immigrants in the US”

Authors
Affiliations
Jaehyun Song

Kansai University

Takeshi Iida

Doshisha University

Yuriko Takahashi

Waseda University

Jesús Tovar Mendoza

Autonomous University of Mexico State

Updated

July 28, 2022

1 Setup

pacman::p_load(tidyverse, list, summarytools, gt)
pacman::p_load_gh("JaehyunSong/BalanceR")
df <- read.csv("ReplicationData_CJPS.csv")

List_df <- df |> 
    filter(Treat != "Direct") |>
    mutate(Treat = ifelse(Treat == "Treat", 1, 0)) |>
    select(-Y_D)

DQ_df <- df |> 
  filter(Treat == "Direct") |>
  select(-Treat)

2 Descriptive Statistics (Tab. 1)

2.1 Direct group

DQ_df |>
    select(-ID) |>
    dfSummary() |>
    print(digits = 3, method = "render", headings = FALSE)
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 Y [integer]
Mean (sd) : 1.7 (1)
min ≤ med ≤ max:
0 ≤ 1 ≤ 5
IQR (CV) : 1 (0.6)
0:14(4.4%)
1:160(50.3%)
2:82(25.8%)
3:50(15.7%)
4:8(2.5%)
5:4(1.3%)
 318 (100.0%) 0 (0.0%)
2 Y_D [integer]
Min : 0
Mean : 0.2
Max : 1
0:255(80.2%)
1:63(19.8%)
 318 (100.0%) 0 (0.0%)
3 Female [integer]
Min : 0
Mean : 0.6
Max : 1
0:120(37.7%)
1:198(62.3%)
 318 (100.0%) 0 (0.0%)
4 Age [numeric]
Mean (sd) : 0.1 (1)
min ≤ med ≤ max:
-1.2 ≤ 0 ≤ 3.8
IQR (CV) : 1.3 (6.7)
 43 distinct values 318 (100.0%) 0 (0.0%)
5 Dual [integer]
Min : 0
Mean : 0.2
Max : 1
0:241(75.8%)
1:77(24.2%)
 318 (100.0%) 0 (0.0%)
6 Educ [integer]
Mean (sd) : 0.4 (1.4)
min ≤ med ≤ max:
-3 ≤ 1 ≤ 3
IQR (CV) : 2 (3.4)
-3:7(2.2%)
-2:8(2.5%)
-1:94(29.6%)
0:31(9.7%)
1:112(35.2%)
2:50(15.7%)
3:16(5.0%)
 318 (100.0%) 0 (0.0%)
7 Married [integer]
Min : 0
Mean : 0.5
Max : 1
0:144(45.3%)
1:174(54.7%)
 318 (100.0%) 0 (0.0%)
8 Income [integer]
Mean (sd) : 0.2 (1.9)
min ≤ med ≤ max:
-3 ≤ 0 ≤ 3
IQR (CV) : 3 (9.8)
-3:30(9.4%)
-2:38(11.9%)
-1:53(16.7%)
0:49(15.4%)
1:60(18.9%)
2:45(14.2%)
3:43(13.5%)
 318 (100.0%) 0 (0.0%)
9 Full [integer]
Min : 0
Mean : 0.5
Max : 1
0:161(50.6%)
1:157(49.4%)
 318 (100.0%) 0 (0.0%)
10 Live [integer]
Mean (sd) : 0.3 (1.5)
min ≤ med ≤ max:
-3 ≤ 0 ≤ 2
IQR (CV) : 3 (4.4)
-3:7(2.2%)
-2:37(11.6%)
-1:52(16.4%)
0:71(22.3%)
1:49(15.4%)
2:102(32.1%)
 318 (100.0%) 0 (0.0%)
11 HTA_10mi [integer]
Mean (sd) : 24.6 (46)
min ≤ med ≤ max:
0 ≤ 5.5 ≤ 217
IQR (CV) : 20.8 (1.9)
 69 distinct values 318 (100.0%) 0 (0.0%)

2.2 Control group

List_df |>
    filter(Treat == 0) |>
    select(-c(ID, Treat)) |>
    dfSummary() |>
    print(digits = 3, method = "render", headings = FALSE)
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 Y [integer]
Mean (sd) : 1.6 (1.1)
min ≤ med ≤ max:
0 ≤ 2 ≤ 4
IQR (CV) : 1 (0.7)
0:51(16.6%)
1:94(30.6%)
2:105(34.2%)
3:39(12.7%)
4:18(5.9%)
 307 (100.0%) 0 (0.0%)
2 Female [integer]
Min : 0
Mean : 0.6
Max : 1
0:121(39.4%)
1:186(60.6%)
 307 (100.0%) 0 (0.0%)
3 Age [numeric]
Mean (sd) : 0.3 (1)
min ≤ med ≤ max:
-1.2 ≤ 0.2 ≤ 3.5
IQR (CV) : 1.3 (3.8)
 41 distinct values 307 (100.0%) 0 (0.0%)
4 Dual [integer]
Min : 0
Mean : 0.3
Max : 1
0:215(70.0%)
1:92(30.0%)
 307 (100.0%) 0 (0.0%)
5 Educ [integer]
Mean (sd) : 0.3 (1.4)
min ≤ med ≤ max:
-3 ≤ 0 ≤ 3
IQR (CV) : 2 (4.9)
-3:6(2.0%)
-2:16(5.2%)
-1:90(29.3%)
0:43(14.0%)
1:92(30.0%)
2:46(15.0%)
3:14(4.6%)
 307 (100.0%) 0 (0.0%)
6 Married [integer]
Min : 0
Mean : 0.6
Max : 1
0:127(41.4%)
1:180(58.6%)
 307 (100.0%) 0 (0.0%)
7 Income [integer]
Mean (sd) : 0.1 (1.9)
min ≤ med ≤ max:
-3 ≤ 0 ≤ 3
IQR (CV) : 3 (18)
-3:34(11.1%)
-2:35(11.4%)
-1:51(16.6%)
0:56(18.2%)
1:48(15.6%)
2:42(13.7%)
3:41(13.4%)
 307 (100.0%) 0 (0.0%)
8 Full [integer]
Min : 0
Mean : 0.4
Max : 1
0:181(59.0%)
1:126(41.0%)
 307 (100.0%) 0 (0.0%)
9 Live [integer]
Mean (sd) : 0.1 (1.6)
min ≤ med ≤ max:
-3 ≤ 0 ≤ 2
IQR (CV) : 3 (13.6)
-3:17(5.5%)
-2:39(12.7%)
-1:63(20.5%)
0:53(17.3%)
1:42(13.7%)
2:93(30.3%)
 307 (100.0%) 0 (0.0%)
10 HTA_10mi [integer]
Mean (sd) : 21.9 (40.5)
min ≤ med ≤ max:
0 ≤ 6 ≤ 217
IQR (CV) : 19 (1.9)
 63 distinct values 307 (100.0%) 0 (0.0%)

2.3 Treatment group

List_df |>
    filter(Treat == 1) |>
    select(-c(ID, Treat)) |>
    dfSummary() |>
    print(digits = 3, method = "render", headings = FALSE)
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 Y [integer]
Mean (sd) : 1.9 (1.3)
min ≤ med ≤ max:
0 ≤ 2 ≤ 5
IQR (CV) : 2 (0.7)
0:44(14.0%)
1:72(22.9%)
2:109(34.7%)
3:57(18.2%)
4:15(4.8%)
5:17(5.4%)
 314 (100.0%) 0 (0.0%)
2 Female [integer]
Min : 0
Mean : 0.6
Max : 1
0:124(39.5%)
1:190(60.5%)
 314 (100.0%) 0 (0.0%)
3 Age [numeric]
Mean (sd) : 0.2 (1.1)
min ≤ med ≤ max:
-1.2 ≤ -0.1 ≤ 4.7
IQR (CV) : 1.4 (5.9)
 45 distinct values 314 (100.0%) 0 (0.0%)
4 Dual [integer]
Min : 0
Mean : 0.3
Max : 1
0:226(72.0%)
1:88(28.0%)
 314 (100.0%) 0 (0.0%)
5 Educ [integer]
Mean (sd) : 0.3 (1.4)
min ≤ med ≤ max:
-3 ≤ 0 ≤ 3
IQR (CV) : 2 (5.3)
-3:7(2.2%)
-2:16(5.1%)
-1:100(31.8%)
0:37(11.8%)
1:86(27.4%)
2:52(16.6%)
3:16(5.1%)
 314 (100.0%) 0 (0.0%)
6 Married [integer]
Min : 0
Mean : 0.6
Max : 1
0:131(41.7%)
1:183(58.3%)
 314 (100.0%) 0 (0.0%)
7 Income [integer]
Mean (sd) : 0.2 (1.9)
min ≤ med ≤ max:
-3 ≤ 0 ≤ 3
IQR (CV) : 3 (8.5)
-3:29(9.2%)
-2:31(9.9%)
-1:69(22.0%)
0:35(11.1%)
1:54(17.2%)
2:55(17.5%)
3:41(13.1%)
 314 (100.0%) 0 (0.0%)
8 Full [integer]
Min : 0
Mean : 0.4
Max : 1
0:174(55.4%)
1:140(44.6%)
 314 (100.0%) 0 (0.0%)
9 Live [integer]
Mean (sd) : 0.2 (1.5)
min ≤ med ≤ max:
-3 ≤ 0 ≤ 2
IQR (CV) : 3 (8.1)
-3:9(2.9%)
-2:40(12.7%)
-1:69(22.0%)
0:57(18.2%)
1:43(13.7%)
2:96(30.6%)
 314 (100.0%) 0 (0.0%)
10 HTA_10mi [integer]
Mean (sd) : 27.7 (48.1)
min ≤ med ≤ max:
0 ≤ 5 ≤ 221
IQR (CV) : 28.8 (1.7)
 69 distinct values 314 (100.0%) 0 (0.0%)

2.4 Balance Check

BlcChk_AOV <- function(Cov) {
  cmd    <- paste0("summary(aov(", Cov,
                   " ~ Treat, data = df))[[1]]$`Pr(>F)`[1]")
  pvalue <- eval(parse(text = cmd))
  
  return(paste0(sprintf("%8s", Cov), ": ", sprintf("%.4f", pvalue), "\n"))
}

for (cov in names(df)[c(3, 5:12)]) {
  cat(BlcChk_AOV(cov))
}
       Y: 0.0005
  Female: 0.8779
     Age: 0.3929
    Dual: 0.2594
    Educ: 0.3983
 Married: 0.5476
  Income: 0.7272
    Full: 0.1100
    Live: 0.1973

Balance check using standardized biases.

BlcChk <- df |>
  BalanceR(group = "Treat", 
           cov   = c("Female", "Age", "Dual", "Educ", "Married", 
                     "Income", "Full", "Live", "HTA_10mi"))

summary(BlcChk)
  Covariate Abs_Maximum_SB
1    Female          3.604
2       Age         11.067
3      Dual         12.973
4      Educ          9.639
5   Married          7.907
6    Income          6.183
7      Full         16.793
8      Live         14.073
9  HTA_10mi         13.148
plot(BlcChk, vline = 25, simplify = TRUE, abs = TRUE) +
    scale_x_continuous(breaks = seq(0, 25, 5), labels = seq(0, 25, 5))

3 Correlation HTA (in 10 miles) and other covariates (Fig. A1)

Cor_df <- df |> 
    select(Female:Live, HTA = HTA_10mi) |>
    cor(method = "spearman") |>
    as_tibble()

names(Cor_df) <- c("Female", "Age", "Citizenship",
                   "Education", "Marital\nStatus", "Income",
                   "Full-time", "Length of\nUS Residence", "# of HTAs\nin 10mi.")

Cor_Plot <- Cor_df |>
    abs() |>
    mutate(Cov = names(Cor_df), .before = Female) |>
    pivot_longer(Female:"# of HTAs\nin 10mi.",
                 names_to = "Cov2",
                 values_to = "Cor") |>
    mutate(Cov = case_when(Cov == "Info_Pol" ~ "Info\n(Pol)",
                           Cov == "Info_HTA" ~ "Info\n(HTA)",
                           TRUE              ~ Cov),
           Cov  = fct_inorder(Cov),
           Cov2 = fct_inorder(Cov2),
           Cov2 = fct_rev(Cov2)) |>
    ggplot(aes(x = Cov, y = Cov2)) +
    geom_tile(aes(fill = Cor), color = "black") +
    geom_text(aes(label = sprintf("%.3f", Cor))) +
    scale_fill_gradient2(low = "white", high = "black") +
    labs(fill = "Spearman's correlation coef.\n(Absolute value)") +
    theme_bw() +
    theme(legend.position = "bottom",
          axis.title      = element_blank(),
          panel.border    = element_blank(),
          panel.grid      = element_blank(),
          axis.ticks      = element_blank())

Cor_Plot

4 Testing of No-Design Effect (Tab. A4)

ict.test(List_df$Y, List_df$Treat, J = 4)

Test for List Experiment Design Effects

Estimated population proportions 
                          est.   s.e.
pi(Y_i(0) = 0, Z_i = 1) 0.0260 0.0289
pi(Y_i(0) = 1, Z_i = 1) 0.1029 0.0394
pi(Y_i(0) = 2, Z_i = 1) 0.0978 0.0338
pi(Y_i(0) = 3, Z_i = 1) 0.0433 0.0217
pi(Y_i(0) = 4, Z_i = 1) 0.0541 0.0128
pi(Y_i(0) = 0, Z_i = 0) 0.1401 0.0196
pi(Y_i(0) = 1, Z_i = 0) 0.2033 0.0345
pi(Y_i(0) = 2, Z_i = 0) 0.2442 0.0382
pi(Y_i(0) = 3, Z_i = 0) 0.0838 0.0280
pi(Y_i(0) = 4, Z_i = 0) 0.0045 0.0185

 Y_i(0) is the (latent) count of 'yes' responses to the control items. Z_i is the (latent) binary response to the sensitive item.

Bonferroni-corrected p-value
If this value is below alpha, you reject the null hypothesis of no design effect. If it is above alpha, you fail to reject the null.

Sensitive Item 1 
               1
Est. SE
J(0, 1) 0.026 0.029
J(1, 1) 0.103 0.039
J(2, 1) 0.098 0.034
J(3, 1) 0.043 0.022
J(4, 1) 0.054 0.013
J(0, 0) 0.140 0.020
J(1, 0) 0.203 0.035
J(2, 0) 0.244 0.038
J(3, 0) 0.084 0.028
J(4, 0) 0.004 0.018
**Diff.-in-means** 0.324 0.095

5 Difference in Means Estimator (Tab. 2)

df |>
    group_by(Treat) |>
    summarise(N     = n(),
              Count = mean(Y),
              SD    = sd(Y),
              SE    = SD / sqrt(N)) |>
    mutate(Treat = factor(Treat, levels = c("Direct", "Control", "Treat"))) |>
    arrange(Treat) |>
    gt() |>
    fmt_number(columns = 3:5, decimals = 3) |>
    cols_align(columns = 1, align = "left")
Treat N Count SD SE
Direct 318 1.654 0.966 0.054
Control 307 1.606 1.087 0.062
Treat 314 1.930 1.285 0.072
 
dim_est <- aov(Y ~ Treat, data = df) |>
    TukeyHSD()

dim_est$Treat
                     diff         lwr       upr        p adj
Direct-Control 0.04822486 -0.16210261 0.2585523 0.8524916337
Treat-Control  0.32407311  0.11308863 0.5350576 0.0009551948
Treat-Direct   0.27584826  0.06671703 0.4849795 0.0057187215
DQ_mean <- mean(DQ_df$Y_D)            # Estimated %
DQ_se   <- t.test(DQ_df$Y_D)$stderr   # Standard Error

6 Multivariate Regression (Tab. 3 and A5)

#####################
# Direct item group #
#####################

## logistic regression
D_Fit <- glm(Y_D5 ~ Female + Age + Dual + Educ + 
                 Married + Income + Full + Live +
                 HTA_10mi,
             family = binomial("logit"), data = DQ_df)

############################
# List experimental groups #
############################

## difference in means (linear model)
ML_Fit1 <- ictreg(Y ~ 1, treat = "Treat", J = 4, 
                  data = List_df, method = "lm",
                  maxIter = 1000000)
summary(ML_Fit1)

## difference in means (maximum likelihood)
ML_Fit2 <- ictreg(Y ~ 1, treat = "Treat", J = 4, 
                  data = List_df, method = "ml",
                  maxIter = 1000000)
summary(ML_Fit2)

## with covariates (maximum likelihood)

### basic model
ML_Fit3 <- ictreg(Y ~ Female + Age + Dual + Educ + 
                     Married + Income + Full + Live +
                     HTA_10mi, treat = "Treat", J = 4, 
                  fit.start = "lm",
                  fit.sensitive = "bayesglm",
                  fit.nonsensitive = "nls",
                  data = List_df, method = "ml",
                  maxIter = 1000000)
summary(ML_Fit3)

### with robust standard error (Blair and Imai 2011)
ML_Fit4 <- ictreg(Y ~ Female + Age + Dual + Educ + 
                     Married + Income + Full + Live +
                     HTA_10mi, treat = "Treat", J = 4, 
                  robust = TRUE,
                  fit.start = "lm",
                  fit.sensitive = "bayesglm",
                  fit.nonsensitive = "nls",
                  data = List_df, method = "ml",
                  maxIter = 1000000)
summary(ML_Fit4)

### with top-biased measurement error
ML_Fit5 <- ictreg(Y ~ Female + Age + Dual + Educ + 
                     Married + Income + Full + Live +
                     HTA_10mi, treat = "Treat", J = 4, 
                  error = "topcode",
                  fit.start = "lm",
                  fit.sensitive = "bayesglm",
                  fit.nonsensitive = "nls",
                  data = List_df, method = "ml",
                  maxIter = 1000000)
summary(ML_Fit5)

### with unifrom measurement error
ML_Fit6 <- ictreg(Y ~ Female + Age + Dual + Educ + 
                     Married + Income + Full + Live +
                     HTA_10mi, treat = "Treat", J = 4, 
                  error = "uniform",
                  fit.start = "lm",
                  fit.sensitive = "bayesglm",
                  fit.nonsensitive = "nls",
                  data = List_df, method = "ml",
                  maxIter = 1000000)
summary(ML_Fit6)
Coef_C <- bind_cols(tibble(Name   = names(ML_Fit3$par.control)),
                    tibble(Coef_1 = c(ML_Fit1$par.control, rep(NA, 9))),
                    tibble(SE_1   = c(ML_Fit1$se.control, rep(NA, 9))),
                    tibble(Coef_2 = c(ML_Fit2$par.control, rep(NA, 9))),
                    tibble(SE_2   = c(ML_Fit2$se.control, rep(NA, 9))),
                    tibble(Coef_3 = ML_Fit3$par.control,
                           SE_3   = ML_Fit3$se.control),
                    tibble(Coef_4 = ML_Fit4$par.control,
                           SE_4   = ML_Fit4$se.control),
                    tibble(Coef_5 = ML_Fit5$par.control,
                           SE_5   = ML_Fit5$se.control),
                    tibble(Coef_6 = ML_Fit6$par.control,
                           SE_6   = ML_Fit6$se.control)) |>
    mutate(Name = if_else(Name == "(Intercept)", "Intercept", Name))

Coef_T <- bind_cols(tibble(Name   = names(ML_Fit3$par.treat)),
                    tibble(Coef_1 = c(ML_Fit1$par.treat, rep(NA, 9))),
                    tibble(SE_1   = c(ML_Fit1$se.treat, rep(NA, 9))),
                    tibble(Coef_2 = c(ML_Fit2$par.treat, rep(NA, 9))),
                    tibble(SE_2   = c(ML_Fit2$se.treat, rep(NA, 9))),
                    tibble(Coef_3 = ML_Fit3$par.treat,
                           SE_3   = ML_Fit3$se.treat),
                    tibble(Coef_4 = ML_Fit4$par.treat,
                           SE_4   = ML_Fit4$se.treat),
                    tibble(Coef_5 = ML_Fit5$par.treat,
                           SE_5   = ML_Fit5$se.treat),
                    tibble(Coef_6 = ML_Fit6$par.treat,
                           SE_6   = ML_Fit6$se.treat)) |>
    mutate(Name = if_else(Name == "(Intercept)", "Intercept", Name))

bind_rows(Coef_T, Coef_C) |>
    gt() |>
    tab_spanner(label = "Model 1", columns = c(Coef_1, SE_1)) |>
    tab_spanner(label = "Model 2", columns = c(Coef_2, SE_2)) |>
    tab_spanner(label = "Model 3", columns = c(Coef_3, SE_3)) |>
    tab_spanner(label = "Model 4", columns = c(Coef_4, SE_4)) |>
    tab_spanner(label = "Model 5", columns = c(Coef_5, SE_5)) |>
    tab_spanner(label = "Model 6", columns = c(Coef_6, SE_6)) |>
    cols_label("Name" = "",
               "Coef_1" = "Coef", "SE_1" = "SE",
               "Coef_2" = "Coef", "SE_2" = "SE",
               "Coef_3" = "Coef", "SE_3" = "SE",
               "Coef_4" = "Coef", "SE_4" = "SE",
               "Coef_5" = "Coef", "SE_5" = "SE",
               "Coef_6" = "Coef", "SE_6" = "SE") |>
    fmt_number(columns = 2:13, decimals = 3) |>
    fmt_missing(columns = 2:13, missing_text = "") |>
    tab_row_group(label = "(Control Item)", rows = 11:20) |>
    tab_row_group(label = "(Sensitive Item)", rows = 1:10)
Warning: The `fmt_missing()` function is deprecated and will soon be removed
* Use the `sub_missing()` function instead
Model 1 Model 2 Model 3 Model 4 Model 5 Model 6
Coef SE Coef SE Coef SE Coef SE Coef SE Coef SE
(Sensitive Item)
Intercept 0.324 0.096 −0.588 0.259 −1.749 0.800 −2.161 1.160 −2.003 1.488 −2.392 2.076
Female 0.411 0.704 0.617 1.461 0.184 1.377 0.481 1.343
Age 0.230 0.292 0.273 0.249 0.247 0.434 0.107 0.589
Dual −0.378 0.754 −0.306 1.502 −1.969 2.276 −2.360 2.436
Educ −0.181 0.286 −0.165 0.526 −0.231 0.522 −0.306 0.702
Married −0.146 0.710 −0.068 1.682 −0.590 1.433 −0.764 1.641
Income 0.022 0.196 0.009 0.377 0.062 0.393 0.179 0.520
Full 1.237 0.682 1.595 0.786 0.700 1.145 0.612 1.291
Live −0.717 0.238 −0.792 0.255 −1.029 0.633 −1.317 0.654
HTA_10mi 0.012 0.006 0.014 0.008 0.027 0.021 0.031 0.017
(Control Item)
Intercept 1.606 0.062 −0.417 0.050 −0.586 0.118 −0.561 0.130 −0.745 0.136 −0.779 0.165
Female 0.184 0.108 0.170 0.145 0.221 0.122 0.252 0.125
Age −0.219 0.053 −0.222 0.053 −0.218 0.054 −0.235 0.062
Dual −0.160 0.116 −0.165 0.153 −0.061 0.129 −0.056 0.147
Educ 0.096 0.044 0.091 0.066 0.138 0.049 0.155 0.057
Married 0.101 0.111 0.095 0.171 0.186 0.136 0.213 0.147
Income 0.092 0.033 0.094 0.052 0.070 0.039 0.088 0.044
Full 0.120 0.109 0.096 0.121 0.150 0.125 0.147 0.131
Live 0.014 0.033 0.017 0.043 0.039 0.036 0.037 0.041
HTA_10mi 0.000 0.001 0.000 0.001 −0.001 0.001 −0.001 0.002

Add row below of the table above

  • Method: Linear, ML, ML, ML, ML, ML
  • Covariates: No, No, Yes, Yes, Yes, Yes
  • Adjustment: No, No, No, Robust SE, Top-code error, Uniform error

7 Comparison of Estimated Proportion of Vote-buying by Models (Fig. A4)

VB_df <- tibble(Model = list(ML_Fit1, ML_Fit2, ML_Fit3, 
                             ML_Fit4, ML_Fit5, ML_Fit6)) |>
    mutate(Pred = map(Model, ~predict(.x, avg = TRUE, se.fit = TRUE)),
           Est  = map(Pred, ~.x$fit),
           SE   = map(Pred, ~.x$se)) |>
    unnest(Est:SE) |>
    select(!Model:Pred) |>
    #add_row(data.frame(Est = DQ_mean, SE  = DQ_se), .before = 1) |>
    mutate(Model = c(#"Direct Question",
        "Difference in Means",
        "Maximum Likelihood\n(without covariates)",
        "Maximum Likelihood\n(with covariates)",
        "Maximum Likelihood\n(with covariates & robust SEs)",
        "Top-biased Error",
        "Unifrom Error"),
        Model = fct_rev(fct_inorder(Model)),
        .before = Est) |>
    mutate(Lower = Est + qnorm(0.05) * SE,
           Upper = Est + qnorm(0.95) * SE)

VB_df |>
    gt() |>
    fmt_number(columns = 2:5, decimals = 3) |>
    cols_align(align = c("left"), columns = 1)
Model Est SE Lower Upper
Difference in Means 0.324 0.096 0.167 0.481
Maximum Likelihood (without covariates) 0.357 0.059 0.259 0.455
Maximum Likelihood (with covariates) 0.320 0.056 0.228 0.413
Maximum Likelihood (with covariates & robust SEs) 0.324 0.168 0.048 0.600
Top-biased Error 0.216 0.079 0.086 0.346
Unifrom Error 0.207 0.074 0.084 0.329
 
VB_plot <- VB_df |>
    ggplot() +
    geom_vline(xintercept = 0, color = "red") +
    geom_pointrange(aes(x = Est, xmin = Lower, xmax = Upper, y = Model)) +
    labs(x = "Estimated Prevalnce w/ 90% CI") +
    coord_cartesian(xlim = c(0, 1)) +
    theme_bw() +
    theme(axis.title.y = element_blank())
VB_plot

8 Effect Sizes (Fig. 1)

Pred_Fem <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                    newdata.diff = mutate(List_df, Female = 0),
                    newdata      = mutate(List_df, Female = 1))
Pred_Fem <- tibble(Value = c("Female", "Male", "Diff."), 
                   Est   = Pred_Fem$fit$V1, 
                   SE    = Pred_Fem$se.fit)

Pred_Age <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                    newdata.diff = mutate(List_df, Age = -1.2),
                    newdata      = mutate(List_df, Age = 4.7))
Pred_Age <- tibble(Value = c("68 yrs", "18 yrs", "Diff."), 
                   Est   = Pred_Age$fit$V1, 
                   SE    = Pred_Age$se.fit)

Pred_Dual <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                     newdata.diff = mutate(List_df, Dual = 0),
                     newdata      = mutate(List_df, Dual = 1))
Pred_Dual <- tibble(Value = c("Dual", "Mexican", "Diff."), 
                    Est   = Pred_Dual$fit$V1, 
                    SE    = Pred_Dual$se.fit)

Pred_Edu <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                    newdata.diff = mutate(List_df, Educ = -3),
                    newdata      = mutate(List_df, Educ = 3))
Pred_Edu <- tibble(Value = c("University", "Elementary", "Diff."), 
                   Est   = Pred_Edu$fit$V1, 
                   SE    = Pred_Edu$se.fit)

Pred_Marr <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                     newdata.diff = mutate(List_df, Married = 0),
                     newdata      = mutate(List_df, Married = 1))
Pred_Marr <- tibble(Value = c("Married", "Unmarried", "Diff."), 
                    Est   = Pred_Marr$fit$V1, 
                    SE    = Pred_Marr$se.fit)

Pred_Inc <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                    newdata.diff = mutate(List_df, Income = -3),
                    newdata      = mutate(List_df, Income = 3))
Pred_Inc <- tibble(Value = c("$100,000~", "~$10,000", "Diff."), 
                   Est   = Pred_Inc$fit$V1, 
                   SE    = Pred_Inc$se.fit)

Pred_Full <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                     newdata.diff = mutate(List_df, Full = 0),
                     newdata      = mutate(List_df, Full = 1))
Pred_Full <- tibble(Value = c("Full-time", "Other", "Diff."), 
                    Est   = Pred_Full$fit$V1, 
                    SE    = Pred_Full$se.fit)

Pred_Live <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                     newdata.diff = mutate(List_df, Live = -3),
                     newdata      = mutate(List_df, Live = 2))
Pred_Live <- tibble(Value = c("20 yrs~", "~1 yrs", "Diff."), 
                    Est   = Pred_Live$fit$V1, 
                    SE    = Pred_Live$se.fit)

Pred_HTA <- predict(ML_Fit4, avg = TRUE, se.fit = TRUE,
                    newdata.diff = mutate(List_df, HTA_10mi = 0),
                    newdata      = mutate(List_df, HTA_10mi = 221))
Pred_HTA <- tibble(Value = c("221", "0", "Diff."), 
                   Est   = Pred_HTA$fit$V1, 
                   SE    = Pred_HTA$se.fit)

Pred_df <- bind_rows(list("Gender" = Pred_Fem,
                          "Age" = Pred_Age,
                          "Citizenship" = Pred_Dual,
                          "Education" = Pred_Edu,
                          "Marital Status" = Pred_Marr,
                          "Income" = Pred_Inc,
                          "Employment Status" = Pred_Full,
                          "Length of US Residency" = Pred_Live,
                          "# of HTA within 10 miles" = Pred_HTA),
                     .id = "Cov") |>
    mutate(Lower = Est + qnorm(0.05) * SE,
           Upper = Est + qnorm(0.95) * SE,
           Cov   = fct_inorder(Cov),
           Value = fct_inorder(Value),
           Value = fct_relevel(Value, "Diff."),
           Value = fct_shift(Value, 1))

gt(Pred_df) |>
    cols_align(align = "left", columns = 1:2) |>
    fmt_number(columns = 3:6, decimals = 3)
Cov Value Est SE Lower Upper
Gender Female 0.365 0.104 0.193 0.536
Gender Male 0.268 0.269 −0.174 0.710
Gender Diff. 0.097 0.203 −0.237 0.430
Age 68 yrs 0.537 0.252 0.122 0.952
Age 18 yrs 0.268 0.160 0.006 0.531
Age Diff. 0.269 0.233 −0.114 0.652
Citizenship Dual 0.290 0.309 −0.217 0.798
Citizenship Mexican 0.339 0.115 0.149 0.529
Citizenship Diff. −0.048 0.224 −0.416 0.319
Education University 0.258 0.332 −0.287 0.804
Education Elementary 0.417 0.184 0.113 0.720
Education Diff. −0.159 0.478 −0.944 0.627
Marital Status Married 0.319 0.269 −0.122 0.761
Marital Status Unmarried 0.330 0.080 0.199 0.462
Marital Status Diff. −0.011 0.268 −0.452 0.430
Income $100,000~ 0.328 0.079 0.199 0.457
Income ~$10,000 0.319 0.349 −0.255 0.893
Income Diff. 0.009 0.360 −0.584 0.601
Employment Status Full-time 0.476 0.167 0.202 0.750
Employment Status Other 0.215 0.162 −0.052 0.482
Employment Status Diff. 0.261 0.102 0.093 0.429
Length of US Residency 20 yrs~ 0.111 0.061 0.011 0.211
Length of US Residency ~1 yrs 0.776 0.247 0.370 1.183
Length of US Residency Diff. −0.665 0.223 −1.033 −0.298
# of HTA within 10 miles 221 0.777 0.136 0.554 1.000
# of HTA within 10 miles 0 0.270 0.181 −0.028 0.568
# of HTA within 10 miles Diff. 0.507 0.245 0.104 0.910
 
Cov_Plot <- Pred_df |>
    mutate(Sig = if_else(Lower * Upper > 0, "Sig", "Insig")) |>
    ggplot() +
    geom_hline(yintercept = 0, color = "red") +
    geom_pointrange(aes(x = Value, y = Est, ymin = Lower, ymax = Upper,
                        alpha = Sig)) +
    scale_alpha_manual(values = c("Sig" = 1, "Insig" = 0.35)) +
    labs(y = "Estimated Prevalnce w/ 90% CI") +
    guides(alpha = "none") +
    facet_wrap(~Cov, scales = "free_x") +
    theme_bw() +
    theme(axis.title.x = element_blank())

Cov_Plot

9 Comparison of Estimated Effect Sizes by Models (Fig. A5)

Pred_base <- tibble(Model = list(ML_Fit3, ML_Fit4, ML_Fit5, ML_Fit6))

Pred_Fem <-  Pred_base |>
    mutate(Pred = pmap(list(Model, 0, 1), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, Female = ..2),
                                newdata      = mutate(List_df, Female = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "Female\n(Min: Male, Max: Female)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_Age <-  Pred_base |>
    mutate(Pred = pmap(list(Model, -1.2, 4.7), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, Age = ..2),
                                newdata      = mutate(List_df, Age = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "Age\n(Min: 18, Max: 68)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_Dual <- Pred_base |>
    mutate(Pred = pmap(list(Model, 0, 1), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, Dual = ..2),
                                newdata      = mutate(List_df, Dual = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "Citizenship\n(Min: Mexican, Max: Dual)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_Edu <- Pred_base |>
    mutate(Pred = pmap(list(Model, -3, 3), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, Educ = ..2),
                                newdata      = mutate(List_df, Educ = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "Education\n(Min: Elementary, Max: University)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_Marr <- Pred_base |>
    mutate(Pred = pmap(list(Model, 0, 1), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, Married = ..2),
                                newdata      = mutate(List_df, Married = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "Marital Status\n(Min: Unmarried, Max: Married)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_Inc <- Pred_base |>
    mutate(Pred = pmap(list(Model, -3, 3), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, Income = ..2),
                                newdata      = mutate(List_df, Income = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "Income\n(Min: ~$10,000, Max: $100,000~)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_Full <- Pred_base |>
    mutate(Pred = pmap(list(Model, 0, 1), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, Full = ..2),
                                newdata      = mutate(List_df, Full = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "Employment Status\n(Min: Other, Max: Full-time)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_Live <- Pred_base |>
    mutate(Pred = pmap(list(Model, -3, 2), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, Live = ..2),
                                newdata      = mutate(List_df, Live = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "Length of US Residency\n(Min: ~1 yrs, Max: 20 yrs~)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_HTA <- Pred_base |>
    mutate(Pred = pmap(list(Model, 0, 221), 
                       ~predict(..1, avg = TRUE, se.fit = TRUE,
                                newdata.diff = mutate(List_df, HTA_10mi = ..2),
                                newdata      = mutate(List_df, HTA_10mi = ..3))),
           Pred = map(Pred, 
                      ~tibble(Cov  = "# of HTAs within 10 miles\n(Min: 0, Max: 221)",
                              Name = c("Min", "Max", "Diff (Max - Min)"), 
                              Coef = .x$fit$V1, 
                              SE   = .x$se.fit)),
           Model = paste0("Fit", 3:6)) |>
    unnest(Pred)

Pred_df_A <- bind_rows(Pred_Fem, Pred_Age, Pred_Dual, Pred_Edu,
                       Pred_Marr, Pred_Inc, Pred_Full, Pred_Live, Pred_HTA) |>
    mutate(Model = case_when(Model == "Fit3" ~ "ML",
                             Model == "Fit4" ~ "Robust ML",
                             Model == "Fit5" ~ "Top-biased",
                             TRUE            ~ "Uniform"),
           Model = fct_rev(fct_inorder(Model)),
           Cov   = fct_inorder(Cov),
           Lower = Coef + qnorm(0.05) * SE,
           Upper = Coef + qnorm(0.95) * SE) 

gt(Pred_df_A) |>
    cols_align(align = "left", columns = 1:3) |>
    fmt_number(columns = 4:7, decimals = 3)
Model Cov Name Coef SE Lower Upper
ML Female (Min: Male, Max: Female) Min 0.349 0.066 0.240 0.458
ML Female (Min: Male, Max: Female) Max 0.281 0.093 0.127 0.435
ML Female (Min: Male, Max: Female) Diff (Max - Min) 0.068 0.113 −0.118 0.253
Robust ML Female (Min: Male, Max: Female) Min 0.365 0.104 0.193 0.536
Robust ML Female (Min: Male, Max: Female) Max 0.268 0.269 −0.174 0.710
Robust ML Female (Min: Male, Max: Female) Diff (Max - Min) 0.097 0.203 −0.237 0.430
Top-biased Female (Min: Male, Max: Female) Min 0.224 0.126 0.017 0.432
Top-biased Female (Min: Male, Max: Female) Max 0.205 0.076 0.080 0.331
Top-biased Female (Min: Male, Max: Female) Diff (Max - Min) 0.019 0.147 −0.222 0.261
Uniform Female (Min: Male, Max: Female) Min 0.225 0.110 0.044 0.407
Uniform Female (Min: Male, Max: Female) Max 0.183 0.071 0.066 0.300
Uniform Female (Min: Male, Max: Female) Diff (Max - Min) 0.042 0.120 −0.156 0.240
ML Age (Min: 18, Max: 68) Min 0.510 0.240 0.115 0.904
ML Age (Min: 18, Max: 68) Max 0.271 0.083 0.135 0.407
ML Age (Min: 18, Max: 68) Diff (Max - Min) 0.238 0.298 −0.251 0.728
Robust ML Age (Min: 18, Max: 68) Min 0.537 0.252 0.122 0.952
Robust ML Age (Min: 18, Max: 68) Max 0.268 0.160 0.006 0.531
Robust ML Age (Min: 18, Max: 68) Diff (Max - Min) 0.269 0.233 −0.114 0.652
Top-biased Age (Min: 18, Max: 68) Min 0.356 0.297 −0.134 0.845
Top-biased Age (Min: 18, Max: 68) Max 0.186 0.089 0.040 0.332
Top-biased Age (Min: 18, Max: 68) Diff (Max - Min) 0.170 0.336 −0.384 0.723
Uniform Age (Min: 18, Max: 68) Min 0.255 0.261 −0.175 0.684
Uniform Age (Min: 18, Max: 68) Max 0.196 0.112 0.012 0.379
Uniform Age (Min: 18, Max: 68) Diff (Max - Min) 0.059 0.341 −0.501 0.619
ML Citizenship (Min: Mexican, Max: Dual) Min 0.277 0.110 0.097 0.458
ML Citizenship (Min: Mexican, Max: Dual) Max 0.339 0.059 0.242 0.437
ML Citizenship (Min: Mexican, Max: Dual) Diff (Max - Min) −0.062 0.119 −0.258 0.134
Robust ML Citizenship (Min: Mexican, Max: Dual) Min 0.290 0.309 −0.217 0.798
Robust ML Citizenship (Min: Mexican, Max: Dual) Max 0.339 0.115 0.149 0.529
Robust ML Citizenship (Min: Mexican, Max: Dual) Diff (Max - Min) −0.048 0.224 −0.416 0.319
Top-biased Citizenship (Min: Mexican, Max: Dual) Min 0.098 0.098 −0.063 0.259
Top-biased Citizenship (Min: Mexican, Max: Dual) Max 0.276 0.098 0.114 0.438
Top-biased Citizenship (Min: Mexican, Max: Dual) Diff (Max - Min) −0.178 0.137 −0.403 0.047
Uniform Citizenship (Min: Mexican, Max: Dual) Min 0.087 0.072 −0.032 0.207
Uniform Citizenship (Min: Mexican, Max: Dual) Max 0.269 0.114 0.082 0.456
Uniform Citizenship (Min: Mexican, Max: Dual) Diff (Max - Min) −0.181 0.153 −0.432 0.070
ML Education (Min: Elementary, Max: University) Min 0.245 0.133 0.026 0.465
ML Education (Min: Elementary, Max: University) Max 0.428 0.169 0.149 0.706
ML Education (Min: Elementary, Max: University) Diff (Max - Min) −0.182 0.284 −0.649 0.284
Robust ML Education (Min: Elementary, Max: University) Min 0.258 0.332 −0.287 0.804
Robust ML Education (Min: Elementary, Max: University) Max 0.417 0.184 0.113 0.720
Robust ML Education (Min: Elementary, Max: University) Diff (Max - Min) −0.159 0.478 −0.944 0.627
Top-biased Education (Min: Elementary, Max: University) Min 0.158 0.119 −0.037 0.354
Top-biased Education (Min: Elementary, Max: University) Max 0.310 0.281 −0.152 0.771
Top-biased Education (Min: Elementary, Max: University) Diff (Max - Min) −0.151 0.369 −0.758 0.455
Uniform Education (Min: Elementary, Max: University) Min 0.142 0.115 −0.047 0.332
Uniform Education (Min: Elementary, Max: University) Max 0.315 0.322 −0.215 0.845
Uniform Education (Min: Elementary, Max: University) Diff (Max - Min) −0.173 0.415 −0.856 0.510
ML Marital Status (Min: Unmarried, Max: Married) Min 0.310 0.084 0.171 0.449
ML Marital Status (Min: Unmarried, Max: Married) Max 0.335 0.077 0.208 0.461
ML Marital Status (Min: Unmarried, Max: Married) Diff (Max - Min) −0.025 0.120 −0.222 0.172
Robust ML Marital Status (Min: Unmarried, Max: Married) Min 0.319 0.269 −0.122 0.761
Robust ML Marital Status (Min: Unmarried, Max: Married) Max 0.330 0.080 0.199 0.462
Robust ML Marital Status (Min: Unmarried, Max: Married) Diff (Max - Min) −0.011 0.268 −0.452 0.430
Top-biased Marital Status (Min: Unmarried, Max: Married) Min 0.189 0.076 0.064 0.314
Top-biased Marital Status (Min: Unmarried, Max: Married) Max 0.252 0.153 0.001 0.503
Top-biased Marital Status (Min: Unmarried, Max: Married) Diff (Max - Min) −0.063 0.166 −0.337 0.210
Uniform Marital Status (Min: Unmarried, Max: Married) Min 0.176 0.066 0.067 0.284
Uniform Marital Status (Min: Unmarried, Max: Married) Max 0.246 0.147 0.004 0.488
Uniform Marital Status (Min: Unmarried, Max: Married) Diff (Max - Min) −0.070 0.159 −0.333 0.192
ML Income (Min: ~$10,000, Max: $100,000~) Min 0.331 0.096 0.173 0.489
ML Income (Min: ~$10,000, Max: $100,000~) Max 0.309 0.129 0.097 0.521
ML Income (Min: ~$10,000, Max: $100,000~) Diff (Max - Min) 0.022 0.196 −0.301 0.345
Robust ML Income (Min: ~$10,000, Max: $100,000~) Min 0.328 0.079 0.199 0.457
Robust ML Income (Min: ~$10,000, Max: $100,000~) Max 0.319 0.349 −0.255 0.893
Robust ML Income (Min: ~$10,000, Max: $100,000~) Diff (Max - Min) 0.009 0.360 −0.584 0.601
Top-biased Income (Min: ~$10,000, Max: $100,000~) Min 0.236 0.169 −0.042 0.513
Top-biased Income (Min: ~$10,000, Max: $100,000~) Max 0.197 0.119 0.001 0.393
Top-biased Income (Min: ~$10,000, Max: $100,000~) Diff (Max - Min) 0.039 0.245 −0.364 0.442
Uniform Income (Min: ~$10,000, Max: $100,000~) Min 0.256 0.199 −0.071 0.583
Uniform Income (Min: ~$10,000, Max: $100,000~) Max 0.163 0.098 0.002 0.324
Uniform Income (Min: ~$10,000, Max: $100,000~) Diff (Max - Min) 0.093 0.266 −0.345 0.532
ML Employment Status (Min: Other, Max: Full-time) Min 0.443 0.088 0.298 0.589
ML Employment Status (Min: Other, Max: Full-time) Max 0.232 0.066 0.123 0.340
ML Employment Status (Min: Other, Max: Full-time) Diff (Max - Min) 0.211 0.108 0.033 0.389
Robust ML Employment Status (Min: Other, Max: Full-time) Min 0.476 0.167 0.202 0.750
Robust ML Employment Status (Min: Other, Max: Full-time) Max 0.215 0.162 −0.052 0.482
Robust ML Employment Status (Min: Other, Max: Full-time) Diff (Max - Min) 0.261 0.102 0.093 0.429
Top-biased Employment Status (Min: Other, Max: Full-time) Min 0.259 0.145 0.020 0.497
Top-biased Employment Status (Min: Other, Max: Full-time) Max 0.184 0.063 0.080 0.289
Top-biased Employment Status (Min: Other, Max: Full-time) Diff (Max - Min) 0.074 0.141 −0.158 0.307
Uniform Employment Status (Min: Other, Max: Full-time) Min 0.238 0.113 0.052 0.423
Uniform Employment Status (Min: Other, Max: Full-time) Max 0.183 0.080 0.051 0.314
Uniform Employment Status (Min: Other, Max: Full-time) Diff (Max - Min) 0.055 0.123 −0.148 0.258
ML Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Min 0.117 0.050 0.034 0.200
ML Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Max 0.755 0.120 0.558 0.952
ML Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Diff (Max - Min) −0.638 0.150 −0.884 −0.391
Robust ML Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Min 0.111 0.061 0.011 0.211
Robust ML Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Max 0.776 0.247 0.370 1.183
Robust ML Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Diff (Max - Min) −0.665 0.223 −1.033 −0.298
Top-biased Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Min 0.058 0.055 −0.032 0.148
Top-biased Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Max 0.704 0.179 0.410 0.998
Top-biased Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Diff (Max - Min) −0.647 0.217 −1.003 −0.290
Uniform Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Min 0.041 0.040 −0.024 0.107
Uniform Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Max 0.758 0.174 0.471 1.044
Uniform Length of US Residency (Min: ~1 yrs, Max: 20 yrs~) Diff (Max - Min) −0.716 0.203 −1.051 −0.382
ML # of HTAs within 10 miles (Min: 0, Max: 221) Min 0.758 0.171 0.477 1.038
ML # of HTAs within 10 miles (Min: 0, Max: 221) Max 0.268 0.064 0.164 0.373
ML # of HTAs within 10 miles (Min: 0, Max: 221) Diff (Max - Min) 0.489 0.200 0.160 0.819
Robust ML # of HTAs within 10 miles (Min: 0, Max: 221) Min 0.777 0.136 0.554 1.000
Robust ML # of HTAs within 10 miles (Min: 0, Max: 221) Max 0.270 0.181 −0.028 0.568
Robust ML # of HTAs within 10 miles (Min: 0, Max: 221) Diff (Max - Min) 0.507 0.245 0.104 0.910
Top-biased # of HTAs within 10 miles (Min: 0, Max: 221) Min 0.910 0.160 0.647 1.172
Top-biased # of HTAs within 10 miles (Min: 0, Max: 221) Max 0.139 0.088 −0.006 0.285
Top-biased # of HTAs within 10 miles (Min: 0, Max: 221) Diff (Max - Min) 0.770 0.226 0.398 1.143
Uniform # of HTAs within 10 miles (Min: 0, Max: 221) Min 0.914 0.124 0.710 1.117
Uniform # of HTAs within 10 miles (Min: 0, Max: 221) Max 0.130 0.072 0.012 0.248
Uniform # of HTAs within 10 miles (Min: 0, Max: 221) Diff (Max - Min) 0.784 0.161 0.518 1.049
 
Cov_plot_A <- Pred_df_A |>
    ggplot(aes(x = Coef, xmin = Lower, xmax = Upper, y = Model, color = Name)) +
    geom_vline(xintercept = 0, color = "red") +
    geom_pointrange(position = position_dodge2(0.75)) +
    labs(x = "Estimated Prevalnce w/ 90% CI", y = "", color = "") +
    guides(color = guide_legend(reverse = TRUE)) +
    facet_wrap(~ Cov) +
    #coord_cartesian(xlim = c(-0.2, 1)) +
    theme_bw() +
    theme(legend.position = "bottom",
          axis.title.y = element_blank())

Cov_plot_A

10 Session Information

sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] BalanceR_0.7.7     gt_0.6.0           summarytools_1.0.1 list_9.2.4        
 [5] sandwich_3.0-2     forcats_0.5.1      stringr_1.4.0      dplyr_1.0.9       
 [9] purrr_0.3.4        readr_2.1.2        tidyr_1.2.0        tibble_3.1.8      
[13] ggplot2_3.3.6      tidyverse_1.3.2   

loaded via a namespace (and not attached):
 [1] nlme_3.1-158        matrixStats_0.62.0  fs_1.5.2           
 [4] lubridate_1.8.0     httr_1.4.3          tools_4.2.1        
 [7] backports_1.4.1     utf8_1.2.2          R6_2.5.1           
[10] DBI_1.1.3           colorspace_2.0-3    withr_2.5.0        
[13] tidyselect_1.1.2    compiler_4.2.1      textshaping_0.3.6  
[16] cli_3.3.0           rvest_1.0.2         pacman_0.5.1       
[19] xml2_1.3.3          sass_0.4.2          labeling_0.4.2     
[22] checkmate_2.1.0     scales_1.2.0        mvtnorm_1.1-3      
[25] quadprog_1.5-8      systemfonts_1.0.4   digest_0.6.29      
[28] minqa_1.2.4         rmarkdown_2.14      base64enc_0.1-3    
[31] pkgconfig_2.0.3     htmltools_0.5.3     lme4_1.1-30        
[34] dbplyr_2.2.1        fastmap_1.1.0       htmlwidgets_1.5.4  
[37] rlang_1.0.4         readxl_1.4.0        pryr_0.1.5         
[40] rstudioapi_0.13     VGAM_1.1-7          farver_2.1.1       
[43] generics_0.1.3      zoo_1.8-10          jsonlite_1.8.0     
[46] googlesheets4_1.0.0 magrittr_2.0.3      rapportools_1.1    
[49] Matrix_1.4-1        Rcpp_1.0.9          munsell_0.5.0      
[52] fansi_1.0.3         abind_1.4-5         lifecycle_1.0.1    
[55] stringi_1.7.8       yaml_2.3.5          MASS_7.3-58        
[58] plyr_1.8.7          gamlss.dist_6.0-3   grid_4.2.1         
[61] crayon_1.5.1        lattice_0.20-45     haven_2.5.0        
[64] splines_4.2.1       pander_0.6.5        hms_1.1.1          
[67] magick_2.7.3        knitr_1.39          pillar_1.8.0       
[70] tcltk_4.2.1         boot_1.3-28         corpcor_1.6.10     
[73] reshape2_1.4.4      codetools_0.2-18    stats4_4.2.1       
[76] magic_1.6-0         reprex_2.0.1        glue_1.6.2         
[79] evaluate_0.15       modelr_0.1.8        vctrs_0.4.1        
[82] nloptr_2.0.3        tzdb_0.3.0          cellranger_1.1.0   
[85] gtable_0.3.0        assertthat_0.2.1    xfun_0.31          
[88] broom_1.0.0         coda_0.19-4         ragg_1.2.2         
[91] googledrive_2.0.0   gargle_1.2.0        arm_1.12-2         
[94] ellipsis_0.3.2