Estimate Voting Probabilities
We use data from the First Round of the Chilean Presidential Election
of 2021, where at each ballot box we have the voters of eight age ranges
(groups), and candidate votes obtained. The function
get_eim_chile() loads this data at either the country,
region, or electoral district level.
library(fastei)
eim_apo <- get_eim_chile(elect_district = "APOQUINDO")
eim_apo
## eim ecological inference model
## Candidates' vote matrix (X) [b x c]:
## C1 C2 C3 C4 C5 C6 C7 C8
## 37M 36 91 19 48 2 7 2 5
## 38M 29 103 17 54 1 6 3 3
## 39M 32 81 13 73 1 5 5 1
## 40M 25 95 21 38 3 4 2 2
## 41M 31 86 19 48 0 4 5 2
## .
## .
## .
## Group-level voter matrix (W) [b x g]:
## X18.19 X20.29 X30.39 X40.49 X50.59 X60.69 X70.79 X80.
## 37M 5 27 11 50 23 39 36 19
## 38M 5 33 11 48 22 42 36 19
## 39M 4 19 11 71 29 30 35 11
## 40M 6 21 6 49 36 20 29 23
## 41M 3 27 6 55 23 29 33 19
## .
## .
## .
As shown in the output, it returns an eim object that contains two
matrices: the number of votes per candidate and the number of votes per
group. The rows of both matrices correspond to the specific ballot-box,
and the columns are candidates and groups, respectively. This eim object
is used as input to run the EM algorithm that estimates the voting
probabilities. In this example, it uses the default method
mult, which is the most efficient in terms of runtime.
Running the algorithm is done by calling run_em().
eim_apo <- run_em(eim_apo)
round(eim_apo$prob, 4)
## C1 C2 C3 C4 C5 C6 C7 C8
## X18.19 0.1802 0.3812 0.0859 0.3191 0.0105 0.0187 0.0040 0.0003
## X20.29 0.2022 0.3427 0.0492 0.3472 0.0086 0.0163 0.0188 0.0149
## X30.39 0.1621 0.4130 0.0418 0.3493 0.0030 0.0081 0.0188 0.0039
## X40.49 0.0919 0.4898 0.0786 0.3271 0.0021 0.0041 0.0029 0.0034
## X50.59 0.1486 0.4141 0.0939 0.2893 0.0011 0.0345 0.0150 0.0035
## X60.69 0.1007 0.4989 0.0907 0.2744 0.0018 0.0218 0.0071 0.0045
## X70.79 0.1051 0.5208 0.0844 0.2355 0.0060 0.0219 0.0172 0.0093
## X80. 0.0759 0.5152 0.1341 0.2169 0.0061 0.0314 0.0118 0.0086
Note that each row corresponds to the probability that a demographic
group (g) voted for a candidate (c). It is
worth noting how the estimated probabilities differ substantially across
groups.
Standard deviation estimates
We can compute the standard deviation of the estimated probabilities
using bootstrapping. This can be done with the function
bootstrap().
eim_apo <- bootstrap(eim_apo, seed = 42, nboot = 30)
round(eim_apo$sd, 4)
## C1 C2 C3 C4 C5 C6 C7 C8
## X18.19 0.0292 0.0688 0.0087 0.0404 0.0043 0.0062 0.0059 0.0018
## X20.29 0.0091 0.0113 0.0042 0.0073 0.0013 0.0037 0.0029 0.0013
## X30.39 0.0097 0.0113 0.0051 0.0126 0.0010 0.0020 0.0033 0.0014
## X40.49 0.0077 0.0134 0.0078 0.0094 0.0013 0.0018 0.0017 0.0017
## X50.59 0.0115 0.0215 0.0091 0.0173 0.0008 0.0066 0.0049 0.0012
## X60.69 0.0105 0.0199 0.0103 0.0146 0.0010 0.0046 0.0031 0.0020
## X70.79 0.0125 0.0220 0.0099 0.0171 0.0029 0.0058 0.0060 0.0037
## X80. 0.0178 0.0310 0.0187 0.0248 0.0041 0.0070 0.0047 0.0036
The standard deviations obtained in the district “Apoquindo” are low
in general. One reason for this is the high number of ballot boxes in
this district. In contrast, standard deviations of estimated
probabilities in districts with fewer ballot boxes, such as “Navidad”,
are larger.
eim_nav <- get_eim_chile(elect_district = "NAVIDAD")
eim_nav <- bootstrap(eim_nav, seed = 42, nboot = 30)
round(eim_nav$sd, 4)
## C1 C2 C3 C4 C5 C6 C7 C8
## X18.19 0.1092 0.2340 0.0002 0.3233 0.1278 0.0955 0.2173 0.0859
## X20.29 0.1556 0.0754 0.0612 0.1184 0.0344 0.0009 0.1017 0.0056
## X30.39 0.0896 0.0708 0.0597 0.0816 0.0209 0.0333 0.0409 0.0040
## X40.49 0.0904 0.1111 0.0462 0.0629 0.0098 0.0521 0.0272 0.0072
## X50.59 0.0815 0.0816 0.0812 0.0795 0.0013 0.0328 0.0158 0.0255
## X60.69 0.1047 0.1017 0.1196 0.0768 0.0152 0.0794 0.0199 0.0010
## X70.79 0.1534 0.1398 0.1374 0.0815 0.0242 0.0447 0.0364 0.0001
## X80. 0.1526 0.0842 0.2060 0.0635 0.1175 0.2250 0.2422 0.1073
It is possible to see the difference in a visual way by plotting the
standard deviations of the estimated probabilities across groups.
library(ggplot2)
library(reshape2)
library(viridis)
plot_district <- function(matrix1, district1, matrix2, district2, sd = FALSE) {
value <- ifelse(sd == FALSE, "prob", "sd")
df1 <- melt(matrix1)
df2 <- melt(matrix2)
df1$Matrix <- district1
df2$Matrix <- district2
combined_df <- rbind(df1, df2)
color <- ifelse(value == "prob", "plasma", "viridis")
# Add text to each cell of the matrix
combined_df$label <- sprintf("%.3f", combined_df$value)
combined_df$text_color <- ifelse(combined_df$value > round(max(combined_df$value) * 0.75 + min(combined_df$value) * 0.25, 2), "black", "white")
districts <- sort(c(district1, district2))
# Call the plot
ggplot(combined_df, aes(x = Var2, y = Var1, fill = value)) +
geom_tile() +
geom_text(aes(label = label, color = text_color), size = 2.5) +
scale_fill_viridis(
name = value,
option = color,
) +
scale_color_identity() +
facet_wrap(~Matrix) +
coord_fixed() +
theme_bw() +
labs(
title = ifelse(value == "prob",
paste("Estimated probabilities in districts:", districts[1], "and", districts[2]),
paste("Standard deviation of estimated probabilities in districts:", districts[1], "and", districts[2])
),
x = "Candidates' votes", y = "Voters' age range", fill = value
)
}
plot_district(
matrix1 = eim_nav$sd, district1 = "Navidad",
matrix2 = eim_apo$sd, district2 = "Apoquindo", sd = TRUE
)
Reduce Estimation Error using Group Aggregation
Demographic groups can be merged to have probability estimates with
lower error. A greedy strategy that maximizes the variability of the
distribution of groups across ballot boxes is used, which ensures
standard deviations are below a specific threshold. The package provides
the following function for the latter:
eim_nav_proxy <- get_agg_proxy(eim_nav, seed = 6, sd_threshold = 0.03, sd_statistic = "maximum")
eim_nav_proxy$group_agg
## [1] 4 8
As shown in the output, the heuristic finds a group aggregation that
merges the first and last four age ranges, namely, macro-groups with
voters aged 18-49, and older than 50. We can evaluate the effectiveness
of this grouping by comparing the mean standard deviation to the
original formulation:
mean(eim_nav$sd) - mean(eim_nav_proxy$sd)
## [1] 0.07432595
It is possible to see the aggregated standard deviations
visually:
plot_matrix <- function(mat, sd = FALSE, y_labels = NULL) {
# Initial configurations
if (!sd) mat <- t(mat)
df <- reshape2::melt(mat)
colnames(df) <- c("Row", "Column", "Value")
df$Row <- factor(df$Row, levels = rev(sort(unique(df$Row))))
df$Column <- factor(df$Column, levels = sort(unique(df$Column)))
if (!sd) {
df$Label <- sprintf("%d", df$Value)
title_text <- "Voters distribution"
x_lab <- "Ballot Box"
y_lab <- "Dem. Group"
fill_lab <- "Voters"
df$text_color <- ifelse(df$Value > 30, "black", "white")
option <- "inferno"
start <- 0.5
limits <- NULL
} else {
df$Label <- sprintf("%.3f", df$Value)
title_text <- "Standard deviation of estimated probabilities on district: Navidad"
x_lab <- "Candidates' votes"
y_lab <- "Voters' age range"
fill_lab <- "sd"
df$text_color <- ifelse(df$Value > 0.13, "black", "white")
option <- "viridis"
start <- 0
limits <- c(0, 0.1)
}
# Plot
p <- ggplot(df, aes(x = Column, y = Row, fill = Value)) +
geom_tile() +
geom_text(aes(label = Label, color = text_color), size = 3) +
scale_color_identity() +
scale_fill_viridis_c(option = option, begin = start, limits = limits) +
coord_fixed() +
theme_bw() +
theme(axis.text.y = element_text(size = 7), axis.text.x = element_text(size = 7)) +
labs(
title = title_text,
x = x_lab,
y = y_lab,
fill = fill_lab
)
# Add custom y-axis labels if provided
if (!is.null(y_labels)) {
p <- p + scale_y_discrete(labels = y_labels)
}
p
}
We can visualize the standard deviations of the estimated
probabilities.
plot_matrix(eim_nav_proxy$sd, sd = TRUE, y_labels = c("X18.49", "X50."))
An exhaustive algorithm is also included in the package. It explores
all combinations of adjacent groups in order to maximize the
log-likelihood subject to having standard deviations below a given
threshold. It might require substantial computation time; therefore it
is recommended to use the default method, “mult”.
eim_nav_opt <- get_agg_opt(eim_nav, seed = 0, sd_threshold = 0.03, sd_statistic = "maximum")
eim_nav_opt$group_agg
## [1] 2 4 8
The optimal group aggregation differs slightly from one obtained
before. In this case, the group aggregation consists in the three age
range: 18-29, 30-49, and older than 70.
We can visualize the standard deviations of the estimated
probabilities.
plot_matrix(eim_nav_opt$sd, sd = TRUE, y_labels = c("X18.29", "X30.49", "X50."))
Test difference between estimates
A relevant question is how significantly different are the
probability estimates of two sets of data, such as two different
districts. For instance, we can compare the estimated probabilities of
the district “LO BARNECHEA” and “LA GRANJA”, whose voting tendencies are
expected to differ.
eim_gra <- get_eim_chile("LA GRANJA")
eim_gra <- run_em(eim_gra)
eim_bar <- get_eim_chile("LO BARNECHEA")
eim_bar <- run_em(eim_bar)
plot_district(eim_gra$prob, "La Granja", eim_bar$prob, "Lo Barnechea")
A Wald test can be applied to each component of the two probability
matrix estimates:
comparison2 <- waldtest(
object1 = eim_gra,
object2 = eim_bar,
method = "mult",
nboot = 30,
seed = 42,
)
round(comparison2$pvals, 3)
## C1 C2 C3 C4 C5 C6 C7 C8
## X18.19 0.000 0.000 0.208 0.000 0.533 0.033 0.443 0.878
## X20.29 0.000 0.000 0.005 0.000 0.000 0.000 0.000 0.823
## X30.39 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.009
## X40.49 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## X50.59 0.000 0.000 0.000 0.000 0.182 0.000 0.000 0.180
## X60.69 0.000 0.000 0.001 0.000 0.069 0.000 0.633 0.426
## X70.79 0.748 0.003 0.000 0.714 0.456 0.923 0.015 0.542
## X80. 0.992 0.154 0.071 0.747 0.978 0.336 0.326 0.370
Conversely, one may examine the alternative scenario in which ballot
boxes within a given region are partitioned into equally sized groups.
In this configuration, both groups are expected to exhibit identical
probability distributions.
eimRM <- get_eim_chile(region = "METROPOLITANA DE SANTIAGO")
n <- nrow(eimRM$X)
n_train <- floor(n * 0.5)
train_indices <- sample(seq_len(n), n_train)
eimRMhalf1 <- eim(X = eimRM$X[train_indices, ], W = eimRM$W[train_indices, ])
eimRMhalf2 <- eim(X = eimRM$X[-train_indices, ], W = eimRM$W[-train_indices, ])
eimRMhalf1 <- run_em(eimRMhalf1)
eimRMhalf2 <- run_em(eimRMhalf2)
plot_district(eimRMhalf1$prob, "Half 1", eimRMhalf2$prob, "Half 2")
comparison2 <- waldtest(
object1 = eimRMhalf1,
object2 = eimRMhalf2,
nboot = 10,
seed = 42
)
round(comparison2$pvals, 4)
## C1 C2 C3 C4 C5 C6 C7 C8
## X18.19 0.1052 0.1768 0.9525 0.0459 0.0615 0.5609 0.3218 0.3814
## X20.29 0.0240 0.0584 0.1939 0.1279 0.3934 0.2638 0.8407 0.0418
## X30.39 0.0173 0.0695 0.8241 0.6153 0.9132 0.9874 0.9704 0.5675
## X40.49 0.9661 0.7262 0.7254 0.0750 0.0199 0.0098 0.3868 0.2005
## X50.59 0.7181 0.4912 0.4623 0.5796 0.6915 0.9020 0.4608 0.6597
## X60.69 0.0684 0.0034 0.6953 0.7221 0.8145 0.5284 0.2461 0.0607
## X70.79 0.9179 0.7549 0.8788 0.4923 0.2300 0.3565 0.7537 0.0595
## X80. 0.9998 0.6619 0.4891 0.6027 0.8372 0.9438 0.8549 0.1317
Simulating Election Results
It is possible to simulate an artificial election and get its values
as an eim object with the simulated probability. This is useful for
testing the package’s performance and comparing the different methods
available.
eim_sim <- simulate_election(num_ballots = 15, num_groups = 2, num_candidates = 3, seed = 42)
eim_sim
## eim ecological inference model
## Candidates' vote matrix (X) [b x c]:
## [,1] [,2] [,3]
## [1,] 23 31 46
## [2,] 22 35 43
## [3,] 30 20 50
## [4,] 24 35 41
## [5,] 24 29 47
## .
## .
## .
## Group-level voter matrix (W) [b x g]:
## [,1] [,2]
## [1,] 75 25
## [2,] 77 23
## [3,] 73 27
## [4,] 80 20
## [5,] 73 27
## .
## .
## .
The real voting probabilities of each group for each candidate can be
obtained as follows:
## [,1] [,2] [,3]
## [1,] 0.2588692 0.3615164 0.3796145
## [2,] 0.2913730 0.1829977 0.5256293
In this case, we can access these values since this is a simulation;
however, this is not possible when working with real data.
The estimated voting probabilities are:
eim_sim <- run_em(eim_sim)
eim_sim$prob
## [,1] [,2] [,3]
## [1,] 0.2398303 0.3292528 0.4309169
## [2,] 0.3014838 0.1827996 0.5157166
It is possible to provide a specific voting probability for each
candidate for each group.
input_probability <- matrix(c(0.9, 0.05, 0.05, 0.2, 0.3, 0.5), nrow = 2, byrow = TRUE)
input_probability
## [,1] [,2] [,3]
## [1,] 0.9 0.05 0.05
## [2,] 0.2 0.30 0.50
eim_sim2 <- simulate_election(
num_ballots = 30, num_groups = 2, num_candidates = 3, seed = 42,
prob = input_probability
)
eim_sim2
## eim ecological inference model
## Candidates' vote matrix (X) [b x c]:
## [,1] [,2] [,3]
## [1,] 69 14 17
## [2,] 75 12 13
## [3,] 68 10 22
## [4,] 64 16 20
## [5,] 72 12 16
## .
## .
## .
## Group-level voter matrix (W) [b x g]:
## [,1] [,2]
## [1,] 79 21
## [2,] 78 22
## [3,] 73 27
## [4,] 70 30
## [5,] 73 27
## .
## .
## .
## [,1] [,2] [,3]
## [1,] 0.9 0.05 0.05
## [2,] 0.2 0.30 0.50
There is a parameter, lambda, that controls for the heterogeneity of
voters’ groups across ballot boxes. A value of lambda = 0 indicates that
voters are all randomly assigned in ballot boxes; in contrast to lambda
= 1, in which case voters are assigned in ballot boxes according to
their demographic group. This is explained in detail in the
documentation in simulate_election().
eim_sim3 <- simulate_election(
num_ballots = 20, num_groups = 4, num_candidates = 2, seed = 42,
lambda = 0.1
)
Having the following estimated probabilities:
## [,1] [,2]
## [1,] 0.3686854 0.6313146
## [2,] 0.6627212 0.3372788
## [3,] 0.2112870 0.7887130
## [4,] 0.1805033 0.8194967
On the other side, it is possible to simulate an heterogeneous
sample
eim_sim4 <- simulate_election(
num_ballots = 20, num_groups = 4, num_candidates = 2, seed = 42,
lambda = 0.9
)
With the underlying probabilities
## [,1] [,2]
## [1,] 0.3731878 0.6268122
## [2,] 0.5431853 0.4568147
## [3,] 0.9365629 0.0634371
## [4,] 0.8711210 0.1288790
