#' If Parties NULL
#'
#' Create 3-letter identifiers for parties if 'parties' parameter in 'divisorMethods' or 'LR.Hamilton' functions is NULL.
#' @return 3-letter identifiers for parties.
#' @examples
#' (...)
#' @export
if.parties.null <- function(x) {
# SET SEED
set.seed(1)
parties <- replicate(x,
paste(sample(LETTERS, 3,
replace = TRUE), collapse = ""))
parties
}
#' Sample Election Data
#'
#' Return a pseudorandom sample (from a truncated log-normal distribution) of n election results using a variaty of paremeters such as: the number of parties (np), the total number of seats (TS) or the number of electoral districts.
#' @param seed A scalar value (a random seed) used to initialize a pseudorandom number generator.
#' @param dist A character value for a probability distribution function used to sample data. 'Uniform' (default) = 'uniform', 'Exponential' = 'exp' and 'Log-Normal' = 'lnorm' distributions are available.
#' @param np A numeric value for the number of parties.
#' @param nd A numeric value for the number of electoral districts.
#' @param ne A numeric value for the number of elections.
#' @param mean A numeric value for the location parameter (a log-normal distribution).
#' @param sd A numeric value for the scale parameter (a log-normal distribution).
#' @param rate A numeric value for the rate parameter (an exponential distribution).
#' @param max A numeric value for the upper limit of the domain for a truncated log-normal distribution used to sample votes for parties at the constituency level.
#' @param TS A numeric value for the total number of seats apportioned among electoral district (if nd (the number of district) > 1), using a selected electoral formula.
#' @param formula_dist A character value that specifies the apportionment method used to divide seats ('TS') among electoral districts.
#'
#' @return A list containing:
#' Votes_Dist_Party,
#' Seats_Dist,
#' Votes_Share_Party,
#' Votes_Total_Dist,
#' Votes_Total_Party,
#' Votes_Total.
#' @examples
#' (...)
#' @export
sampleElectionData <-
function (seed = 0,
dist = "uniform",
np,
nd,
ne,
mean = NULL,
sd = NULL,
rate = NULL,
max,
TS,
formula_dist,
...) {
# SET SEED
set.seed(seed)
x = array(dim = c(np, nd, ne))
csum = list()
rsum = list()
tsum = list()
seats_dist = list()
votes_share = list()
switch (dist,
uniform = {
for (j in seq(1, ne, by = 1)) {
for (i in seq(1, nd, by = 1)) {
x[, i, j] <- sort(floor(runif(
np, min = 0, max = max
)))
}
}
},
lnorm = {
for (j in seq(1, ne, by = 1)) {
for (i in seq(1, nd, by = 1)) {
x[, i, j] <-
sort(floor(
truncdist::rtrunc(
np,
spec = "lnorm",
a = 0,
b = max,
meanlog = mean,
sdlog = sd
)
))
}
}
},
exp = {
for (j in seq(1, ne, by = 1)) {
for (i in seq(1, nd, by = 1)) {
x[, i, j] <-
sort(floor(
truncdist::rtrunc(
np,
spec = "exp",
a = 0,
b = max,
rate = rate
)
))
}
}
})
if (nd > 1) {
for (i in seq(1, ne, by = 1)) {
csum[[i]] <- apply(x[, , i], 2, sum)
}
for (i in seq(1, ne, by = 1)) {
rsum[[i]] <- apply(x[, , i], 1, sum)
}
for (i in seq(1, ne, by = 1)) {
tsum[[i]] <- sum(x[, , i])
}
}
else {
for (i in seq(1, ne, by = 1)) {
csum[[i]] <- sum(x[, , i])
}
for (i in seq(1, ne, by = 1)) {
rsum[[i]] <- x[, , i]
}
for (i in seq(1, ne, by = 1)) {
tsum[[i]] <- sum(x[, , i])
}
}
# Use apportionment functions with order_name = FALSE
switch (
formula_dist,
hamilton = {
for (i in seq(1, ne, by = 1)) {
seats_dist[[i]] <-
LR_Hamilton(
parties = if.parties.null(nd),
votes = 1 * csum[[i]],
seats = TS,
order_name = F
)[[2]]
}
},
ad = {
for (i in seq(1, ne, by = 1)) {
seats_dist[[i]] <-
divisorMethods(
parties = if.parties.null(nd),
votes = 1 * csum[[i]],
seats = TS,
method = "ad",
order_name = F
)[[2]]
}
},
dh = {
for (i in seq(1, ne, by = 1)) {
seats_dist[[i]] <-
divisorMethods(
parties = if.parties.null(nd),
votes = 1 * csum[[i]],
seats = TS,
method = "dh",
order_name = F
)[[2]]
}
},
hh = {
for (i in seq(1, ne, by = 1)) {
seats_dist[[i]] <-
divisorMethods(
parties = if.parties.null(nd),
votes = 1 * csum[[i]],
seats = TS,
method = "hh",
order_name = F
)[[2]]
}
}
)
for (i in seq(1, ne, by = 1)) {
votes_share[[i]] <- rsum[[i]] / tsum[[i]]
}
sample <-
list(
Votes_Dist_Party = x,
Seats_Dist = seats_dist,
Votes_Share_Party = votes_share,
Votes_Total_Dist = csum,
Votes_Total_Party = rsum,
Votes_Total = tsum,
Params = c(ne, nd, np, TS)
)
out <- sample
}
# ----
#' Simulate Elections under Proportional Representation
#'
#' The function simulates election results and computes a variaty of disproportionality measures.
#'
#' @return A list containg the following data objects:
#' Seat_Excess,
#' Apportionment,
#' Disproportionality_per_elec,
#' Summary.
#' @examples
#' (...)
#' @export
simulate_E <-
function (seed,
dist = "lnorm",
np,
nd,
ne,
mean,
sd,
rate,
max,
TS,
formula,
formula_dist,
threshold,
threshold_country,
...) {
# SET SEED
set.seed(seed)
sample <-
sampleElectionData(
seed = seed,
dist = dist,
np = np,
nd = nd,
ne = ne,
mean = mean,
sd = sd,
rate = rate,
max = max,
TS = TS,
formula_dist = formula_dist
)
# Seat apportionment per district
# Return list (Party Seats SeatShare Votes VoteShare id elec dist distTS)
apportionment <-
.ProportionalRepresentation(
sample = sample,
formula = formula,
threshold = threshold,
threshold_country = threshold_country
)
apportionment$Party <- as.character(apportionment$Party)
# Determine number of votes by election
# Return df (VoteShareTotalParty VotesTotalParty elec Party)
{
vote_share <- vector("list", ne)
for (i in seq(1, ne, by = 1)) {
vote_share[[i]] = data.frame(
VoteShareTotalParty = round(sample$Votes_Share_Party[[i]], 4),
VotesTotalParty = as.integer(sample$Votes_Total_Party[[i]]),
elec = paste("e", i, sep = "")
)
}
out <- vote_share
out
for (i in seq(1, ne, by = 1)) {
out[[i]] <- dplyr::mutate(out[[i]], Party = if.parties.null(np))
out[[i]] <-
dplyr::arrange(out[[i]], desc(VoteShareTotalParty))
}
vote_share <- data.table::rbindlist(out)
}
apportionment <-
dplyr::left_join(apportionment, vote_share, by = c("Party", "elec"))
## Group districts in elections
apportionment_sum1 <-
dplyr::group_by(apportionment, elec, Party)
apportionment_sum2 <-
dplyr::summarise(
apportionment_sum1,
seats = sum(Seats),
seat_perc = sum(Seats) / TS,
distTS = distTS[1],
TS = TS[1]
)
apportionment_sum3 <-
dplyr::left_join(apportionment_sum2, vote_share, by = c("Party", "elec"))
## Compute 'ideal' shares of seats
seats_ideal <- .seatsIdeal(ne, nd, np, sample)
## Merge
merged <-
dplyr::left_join(apportionment_sum3, seats_ideal, by = c("Party", "elec"))
## Seat excesses
seat_excess <-
dplyr::mutate(
merged,
Party = as.factor(Party),
SE1_i = signif(seats - VoteShareTotalParty * TS, 3),
SE2_i = signif(seat_perc - VoteShareTotalParty, 3),
SE2_i_pp = signif(seat_perc - SeatShareIdeal, 3),
SeatShare = signif(seat_perc, 3),
VoteShare = signif(VoteShareTotalParty, 3),
StandardQuota = (SeatShareIdeal) * TS,
RSE2_i = signif((seat_perc - SeatShareIdeal) / SeatShareIdeal, 3)
)
seat_excess <-
dplyr::select(
seat_excess,
PartyID = Party,
ElectionID = elec, # ElectionID <- elec
Seats = seats,
SeatShare,
Votes = VotesTotalParty,
VoteShare,
SQ = StandardQuota,
SE1_i,
SE2_i,
SE2_i_pp,
RSE2_i
)
seat_excess <-
dplyr::mutate(seat_excess, Seats = as.integer(Seats))
### Disproportionality indexes
disp <- dplyr::group_by(seat_excess, ElectionID)
disp <-
dplyr::summarise(
disp,
meanRSE2 = signif(sum(abs(RSE2_i)) / np, digits = 3),
LHI = signif(1 / 2 * sum(abs(
SeatShare - VoteShare
)), digits = 4),
GHI = signif(sqrt(1 / 2 * sum((SeatShare - VoteShare) ^ 2
)), digits = 3),
SLI = signif(sum((
SeatShare - VoteShare
) ^ 2 / (VoteShare)), digits = 4),
ENPP = signif(1 / sum((SeatShare) ^ 2), digits = 3),
NPP = sum(SeatShare > 0)
)
summary <-
psych::describe(
disp[2:5],
quant = c(.10, .25, .75, .90),
IQR = TRUE,
skew = FALSE
)
# ----
### Results
out <-
list(
Seat_Excess = seat_excess,
Apportionment = apportionment,
Disproportionality_per_elec = disp,
Summary = summary
)
return(out)
}
# ----
#' Per-party Disproportionality Measures for Varying District Sizes and Numbers of Parties
#'
#' The function computes per-party disproportionality measures using simulated election data.
#' @export
simulate_Disp <-
function(seed = 0,
dist = "lnorm",
np,
nd = 1,
ne,
rate = 1 / 25000,
mean = 10,
sd = 1.2,
max = 100000,
formula,
formula_dist = "hh",
threshold = 0,
threshold_country = 0,
minTS = 3,
maxTS = 20,
jump = 2,
...) {
# Declare vars
sb_bw <- list()
sim <- list()
# Data
if (nd >= minTS) {
stop(
"The total number of districts ('nd') needs to be less than the minimal total number of seats ('minTS')."
)
}
else {
# set.seed(seed = seed)
for (i in seq(from = minTS, to = maxTS, by = jump)) {
sim[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
dist = dist,
np = np,
nd = nd,
ne = ne,
mean = mean,
sd = sd,
rate = rate,
max = max,
TS = i,
formula = formula,
formula_dist = formula_dist,
threshold = threshold,
threshold_country = threshold_country
)
sb_bw[[i + (1 - minTS)]] <-
dplyr::mutate(
sim[[i + (1 - minTS)]][[2]],
method = formula,
TS = i,
VoteShare = VoteShare,
SE1_i = Seats - VoteShare * distTS,
SE2_i = SeatShare - VoteShare
)
}
sb_bw <- dplyr::bind_rows(sb_bw) # return data frame
sb_bw_temp <-
dplyr::group_by(sb_bw, TS, elec, Party)
sb_bw_temp <-
dplyr::summarise(
sb_bw_temp,
SeatTotal = sum(Seats),
SeatShareTotal = sum(Seats) / TS[1],
SE2T_i = SeatShareTotal[1] - VoteShareTotalParty[1]
)
sb_bw <- dplyr::left_join(sb_bw, sb_bw_temp, by = c("Party", "elec", "TS"))
ese <- dplyr::group_by(sb_bw, TS, Party)
ese <-
dplyr::summarise(ese,
SB1_i = mean(Seats - VoteShare * distTS),
V = sum(Votes))
ese2 <- dplyr::group_by(sb_bw, TS, Party)
ese2 <-
dplyr::summarise(ese2,
SB2_i = mean(SeatShare - VoteShare),
V = sum(Votes))
ese_mean <- dplyr::group_by(sb_bw, TS, Party)
ese_mean <-
dplyr::summarise(ese_mean,
SB1_i = mean(Seats - VoteShare * distTS),
V = sum(Votes))
ese_mean <- dplyr::group_by(ese_mean, Party)
ese_mean <-
dplyr::summarise(ese_mean, ESB1 = mean(SB1_i), TV = sum(V))
# Return list
bias_data <-
list(
sb_bw = sb_bw,
ese = ese,
ese2 = ese2,
ese_mean = ese_mean,
sim = sim
)
return(bias_data)
}
}
# ----
#' Plot Disproportionality Measures Related to District Sizes and Numbers of Parties
#'
#' The function plots the relationships between, on the one hand, values of a variaty of disproportionality measures and, on the other hand, district sizes or numbers of parties.
#' @export
plot_Disp <-
function(bias_data, tse = c(0, 5 / 12, -1 / 12, -4 / 12), ...)
{
# Plots
sb_bw_plot1 <-
ggplot2::ggplot(data = bias_data$sb_bw) + ggplot2::geom_boxplot(
ggplot2::aes(
x = Party,
y = (SeatShare * distTS - VoteShare * distTS),
fill = factor(TS)
),
lwd = 0.25,
fatten = 0.4,
outlier.size = 0.3
) + ggplot2::ylab("SE1_i(DM)") + viridis::scale_fill_viridis(
discrete = TRUE,
name = "DM",
option = "D",
begin = 0.5
) + ggplot2::geom_hline(yintercept = tse) + ggplot2::theme_classic() + ggplot2::theme(
panel.grid.major = ggplot2::element_line(size = .3, color = "red"),
#increase size of axis lines
axis.line = ggplot2::element_line(size =
.3, color = "black"),
axis.ticks = ggplot2::element_line(size =
0.3, color = "black"),
#increase the font size
text = ggplot2::element_text(size =
12)
)
sb_bw_plot2 <-
ggplot2::ggplot(data = bias_data$sb_bw) + ggplot2::geom_boxplot(
ggplot2::aes(
x = Party,
y = (Seats - VoteShare * distTS),
fill = factor(TS)
),
lwd = 0.25,
fatten = 0.4,
outlier.size = 0.3
) + ggplot2::ylab("SE1_i(DM)") + viridis::scale_fill_viridis(
discrete = TRUE,
name = "DM",
option = "D",
begin = 0.5
) + ggplot2::geom_hline(yintercept = tse) + ggplot2::theme_classic() + ggplot2::theme(
panel.grid.major = ggplot2::element_line(size = .3, color = "red"),
#increase size of axis lines
axis.line = ggplot2::element_line(size =
.3, color = "black"),
axis.ticks = ggplot2::element_line(size =
.3, color = "black"),
#increase the font size
text = ggplot2::element_text(size =
12)
)
sb_bw_plot3 <-
ggplot2::ggplot(data = bias_data$sb_bw) + ggplot2::geom_boxplot(
ggplot2::aes(
x = "",
y = (SeatShare - VoteShare),
fill = factor(TS)
),
lwd = 0.25,
fatten = 0.4,
outlier.size = 0.3
) + ggplot2::facet_wrap(~ Party) + ggplot2::ylab("SE2_i(DM)") + viridis::scale_fill_viridis(
discrete = TRUE,
name = "DM",
option = "D",
begin = 0.6
) + ggplot2::geom_hline(yintercept = c(0, 0.1), colour = "blue") + ggplot2::theme_classic() + ggplot2::theme(
panel.grid.major = ggplot2::element_line(size = .3, color = "red"),
#increase size of axis lines
axis.line = ggplot2::element_line(size =
.3, color = "black"),
axis.ticks = ggplot2::element_line(size =
.3, color = "black"),
#increase the font size
text = ggplot2::element_text(size =
12)
) + ggplot2::stat_summary(
ggplot2::aes(
x = "",
y = (SeatShare - VoteShare),
fill = factor(TS)
),
fun.y = "mean",
shape = 3,
geom = "point",
size = 2,
position = ggplot2::position_dodge(width = 0.75),
color = "black"
)
ese_plot <-
ggplot2::ggplot(data = bias_data$ese) + ggplot2::geom_point(
ggplot2::aes(
x = Party,
y = SB1_i,
colour = factor(TS)
),
size = 4,
alpha = 1 / 2
) + ggplot2::ylab("B_i1(DM)") + ggplot2::facet_grid(~ V) + viridis::scale_color_viridis(name =
"DM", discrete = TRUE) + ggplot2::theme_classic() + ggplot2::geom_hline(yintercept = tse)
ese_plot2 <-
ggplot2::ggplot(data = bias_data$ese2) + ggplot2::geom_point(
ggplot2::aes(
x = Party,
y = SB2_i,
colour = factor(TS)
),
size = 4,
alpha = 1 / 2
) + ggplot2::ylab("B_i2(DM)") + ggplot2::facet_grid(~ V) + viridis::scale_color_viridis(name =
"DM", discrete = TRUE) + ggplot2::theme_classic() + ggplot2::geom_hline(yintercept = c(0))
ese_mean_plot <-
ggplot2::ggplot(data = bias_data$ese_mean) + ggplot2::geom_point(
ggplot2::aes(
x = Party,
y = ESB1,
colour = factor(TV)
),
size = 4,
alpha = 1 / 2
) + ggplot2::ylab("B_i1(DM)") + ggplot2::theme_classic() + ggplot2::geom_hline(yintercept = tse) + viridis::scale_color_viridis(name =
"TV", discrete = TRUE)
# Return list
sb <-
list(sb_bw_plot2,
sb_bw_plot3,
ese_plot,
ese_plot2,
ese_mean_plot)
return(sb)
}
# ----
#' Aggregate-Level Disproportionality Measures
#'
#' The function computes aggregate-level measures of disproportionality. It also models relationships between values of aggregate-level disproportionality measures and district sizes (if model == TRUE).
#' @export
Disp2 <- function(seed = 0,
np = 3,
nd = 1,
ne = 100,
dist = "lnorm",
rate = 1 / 25000,
mean = 10,
sd = 1.2,
max = 100000,
minTS = 3,
maxTS = 20,
jump = 1,
threshold = 0,
threshold_country = 0,
start_C = 0.2,
start_alpha = -0.2,
model = TRUE,
...)
{
dsl <- list()
dmsl <- list()
ddh <- list()
dhh <- list()
dad <- list()
dhamilton <- list()
dimp <- list()
dda <- list()
dh_seatbias <- list()
sl_seatbias <- list()
msl_seatbias <- list()
h_seatbias <- list()
ad_seatbias <- list()
hh_seatbias <- list()
imp_seatbias <- list()
da_seatbias <- list()
out <- list()
for (i in seq(from = minTS, to = maxTS, by = jump)) {
dh_seatbias[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
np = np,
nd = nd,
ne = ne,
dist = dist,
rate = rate,
mean = mean,
sd = sd,
max = max,
TS = i,
formula = "dh",
formula_dist = "hh",
threshold = threshold,
threshold_country = threshold_country
)
ddh[[i]] <-
dplyr::mutate(dh_seatbias[[i + (1 - minTS)]][[3]],
method = "DH",
DM = i,
NP = np)
}
for (i in seq(from = minTS, to = maxTS, by = jump)) {
sl_seatbias[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
np = np,
nd = nd,
ne = ne,
dist = dist,
rate = rate,
mean = mean,
sd = sd,
max = max,
TS = i,
formula = "sl",
formula_dist = "hh",
threshold = threshold,
threshold_country = threshold_country
)
dsl[[i]] <-
dplyr::mutate(sl_seatbias[[i + (1 - minTS)]][[3]],
method = "SL",
DM = i,
NP = np)
}
for (i in seq(from = minTS, to = maxTS, by = jump)) {
msl_seatbias[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
np = np,
nd = nd,
ne = ne,
dist = dist,
rate = rate,
mean = mean,
sd = sd,
max = max,
TS = i,
formula = "msl",
formula_dist = "hh",
threshold = threshold,
threshold_country = threshold_country
)
dmsl[[i]] <-
dplyr::mutate(msl_seatbias[[i + (1 - minTS)]][[3]],
method = "MSL",
DM = i,
NP = np)
}
for (i in seq(from = minTS, to = maxTS, by = jump)) {
h_seatbias[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
np = np,
nd = nd,
ne = ne,
dist = dist,
rate = rate,
mean = mean,
sd = sd,
max = max,
TS = i,
formula = "hamilton",
formula_dist = "hh",
threshold = threshold,
threshold_country = threshold_country
)
dhamilton[[i]] <-
dplyr::mutate(h_seatbias[[i + (1 - minTS)]][[3]],
method = "H",
DM = i,
NP = np)
}
for (i in seq(from = minTS, to = maxTS, by = jump)) {
hh_seatbias[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
np = np,
nd = nd,
ne = ne,
dist = dist,
rate = rate,
mean = mean,
sd = sd,
max = max,
TS = i,
formula = "hh",
formula_dist = "hh",
threshold = threshold,
threshold_country = threshold_country
)
dhh[[i]] <-
dplyr::mutate(hh_seatbias[[i + (1 - minTS)]][[3]],
method = "HH",
DM = i,
NP = np)
}
for (i in seq(from = minTS, to = maxTS, by = jump)) {
ad_seatbias[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
np = np,
nd = nd,
ne = ne,
dist = dist,
rate = rate,
mean = mean,
sd = sd,
max = max,
TS = i,
formula = "ad",
formula_dist = "hh",
threshold = threshold,
threshold_country = threshold_country
)
dad[[i]] <-
dplyr::mutate(ad_seatbias[[i + (1 - minTS)]][[3]],
method = "A",
DM = i,
NP = np)
}
# ---- Danish and HA Imperiali
for (i in seq(from = minTS, to = maxTS, by = jump)) {
imp_seatbias[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
np = np,
nd = nd,
ne = ne,
dist = dist,
rate = rate,
mean = mean,
sd = sd,
max = max,
TS = i,
formula = "imperiali",
formula_dist = "hh",
threshold = threshold,
threshold_country = threshold_country
)
dimp[[i]] <-
dplyr::mutate(imp_seatbias[[i + (1 - minTS)]][[3]],
method = "Imperiali",
DM = i,
NP = np)
}
for (i in seq(from = minTS, to = maxTS, by = jump)) {
da_seatbias[[i + (1 - minTS)]] <-
simulate_E(
seed = seed,
np = np,
nd = nd,
ne = ne,
dist = dist,
rate = rate,
mean = mean,
sd = sd,
max = max,
TS = i,
formula = "danish",
formula_dist = "hh",
threshold = threshold,
threshold_country = threshold_country
)
dda[[i]] <-
dplyr::mutate(da_seatbias[[i + (1 - minTS)]][[3]],
method = "Danish",
DM = i,
NP = np)
}
# ----
lghi_dh <- dplyr::bind_rows(ddh)
lghi_sl <- dplyr::bind_rows(dsl)
lghi_msl <- dplyr::bind_rows(dmsl)
lghi_hh <- dplyr::bind_rows(dhh)
lghi_ad <- dplyr::bind_rows(dad)
lghi_hamilton <- dplyr::bind_rows(dhamilton)
lghi_imp <- dplyr::bind_rows(dimp)
lghi_da <- dplyr::bind_rows(dda)
# ----
# Models
## DH
if (model == TRUE) {
# model_dh <- glm( GHI ~ DM, family = gaussian(link = "log"), data = dplyr::filter(lghi_dh, method == "DH") )
model_dh <-
nls(
GHI ~ C * exp(alpha * DM),
start = list(C = start_C, alpha = start_alpha),
data = dplyr::filter(lghi_dh, method == "DH")
)
lghi_dh = dplyr::mutate(lghi_dh, GHI_predicted = predict(model_dh))
# ----
## SL
# model_sl <- lm( I(log(GHI)) ~ DM, data = dplyr::filter(lghi_sl, method == "SL") )
model_sl <-
nls(
GHI ~ C * exp(alpha * DM),
start = list(C = start_C, alpha = start_alpha),
data = dplyr::filter(lghi_sl, method == "SL")
)
lghi_sl = dplyr::mutate(lghi_sl, GHI_predicted = predict(model_sl))
# ----
## MSL
model_msl <-
nls(
GHI ~ C * exp(alpha * DM),
start = list(C = start_C, alpha = start_alpha),
data = dplyr::filter(lghi_msl, method == "MSL")
)
lghi_msl = dplyr::mutate(lghi_msl, GHI_predicted = predict(model_msl))
# ----
## H
model_h <-
nls(
GHI ~ C * exp(alpha * DM),
start = list(C = start_C, alpha = start_alpha),
data = dplyr::filter(lghi_hamilton, method == "H")
)
lghi_hamilton = dplyr::mutate(lghi_hamilton, GHI_predicted = predict(model_h))
# ----
## Danish
model_danish <-
nls(
GHI ~ C * exp(alpha * DM),
start = list(C = start_C, alpha = start_alpha),
data = dplyr::filter(lghi_da, method == "Danish")
)
lghi_da = dplyr::mutate(lghi_da, GHI_predicted = predict(model_danish))
# ----
## Imperiali
model_imperiali <-
nls(
GHI ~ C * exp(alpha * DM),
start = list(C = start_C, alpha = start_alpha),
data = dplyr::filter(lghi_imp, method == "Imperiali")
)
lghi_imp = dplyr::mutate(lghi_imp, GHI_predicted = predict(model_imperiali))
# ----
## Adams
model_ad <-
nls(
GHI ~ C * exp(alpha * DM),
start = list(C = start_C, alpha = start_alpha),
data = dplyr::filter(lghi_ad, method == "A")
)
lghi_ad = dplyr::mutate(lghi_ad, GHI_predicted = predict(model_ad))
# ----
## HH
model_hh <-
nls(
GHI ~ C * exp(alpha * DM),
start = list(C = start_C, alpha = start_alpha),
data = dplyr::filter(lghi_hh, method == "HH")
)
lghi_hh = dplyr::mutate(lghi_hh, GHI_predicted = predict(model_hh))
# ----
}
lghi_all <-
dplyr::bind_rows(lghi_dh,
lghi_sl,
lghi_msl,
lghi_hamilton,
lghi_da,
lghi_imp,
lghi_ad,
lghi_hh)
lghi_all <- lghi_all[gtools::mixedorder(lghi_all$ElectionID), ]
if (model == TRUE) {
out <-
list(
summary = lghi_all,
Model_DH = model_dh,
Model_SL = model_sl,
Model_MSL = model_msl,
Model_Hamilton = model_h,
Model_Danish = model_danish,
Model_Imperiali = model_imperiali,
DH = dh_seatbias,
SL = sl_seatbias,
MSL = msl_seatbias,
H = h_seatbias,
AD = ad_seatbias,
HH = hh_seatbias,
I = imp_seatbias,
D = da_seatbias
)
} else {
out <-
list(
summary = lghi_all,
DH = dh_seatbias,
SL = sl_seatbias,
MSL = msl_seatbias,
H = h_seatbias,
AD = ad_seatbias,
HH = hh_seatbias,
I = imp_seatbias,
D = da_seatbias
)
}
return(out)
}
# ----
#' Plot Aggregate-Level Disproportionality Measures
#'
#' The function plots relationships between values of aggregate-level disproportionality measures and district sizes.
#' @export
plot_Disp2 <-
function(data = NULL,
methods = c("DH", "SL", "H", "Imperiali"),
...) {
lghi_all <- data[['summary']]
lghi_all <- filter(lghi_all, lghi_all$method %in% methods)
# Plots
# plot_GHI <-
# ggplot2::ggplot(data = lghi_all) + ggplot2::geom_boxplot(
# ggplot2::aes(
# x = factor(DM),
# y = GHI,
# fill = factor(method)
# ),
# lwd = 0.25,
# fatten = 0.4,
# outlier.size = 0.6
# ) + viridis::scale_fill_viridis(option = "C",
# discrete = TRUE,
# begin = 0.4) + ggplot2::xlab("DM") + ggplot2::ylab("GHI") + ggplot2::labs(fill = "Method") + ggplot2::geom_hline(
# yintercept = c(0.1),
# size = 0.35,
# linetype = "longdash",
# colour = "blue"
# ) + ggplot2::theme_classic() + ggplot2::theme(
# panel.grid.major = ggplot2::element_line(size = 0.35, color = "red"),
# #increase size of axis lines
# axis.line = ggplot2::element_line(size =
# 0.35, color = "black"),
# axis.ticks = ggplot2::element_line(size =
# 0.35, color = "black"),
# #Adjust legend position to maximize space, use a vector of proportion
# #across the plot and up the plot where you want the legend.
# #You can also use "left", "right", "top", "bottom", for legends on t
# #he side of the plot
# legend.position = c(.85, .7),
# #increase the font size
# text = ggplot2::element_text(size =
# 12)
# )
if ("GHI_predicted" %in% names(lghi_all) == TRUE) {
plot_GHI <-
ggplot2::ggplot(data = lghi_all) + ggplot2::geom_boxplot(
ggplot2::aes(x = as.factor(DM),
y = GHI),
lwd = 0.25,
fatten = 0.4,
outlier.size = 0.6
) + ggplot2::geom_line(
ggplot2::aes(
x = as.factor(DM),
y = GHI_predicted,
group = 1
),
size = 0.35,
colour = "blue"
) + ggplot2::facet_wrap(~ method) + viridis::scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.4) + ggplot2::xlab("DM") + ggplot2::ylab("GHI") + ggplot2::labs(fill = "Method") + ggplot2::geom_hline(
yintercept = c(0.1, 0.05),
size = 0.35,
linetype = "longdash",
colour = "blue"
) + ggplot2::geom_vline(
xintercept = c(3, 7),
size = 0.45,
linetype = "longdash",
colour = "green"
) + ggplot2::theme_classic() + ggplot2::theme(
panel.grid.major.y = ggplot2::element_line(size = 0.1, color = "red"),
axis.line = ggplot2::element_line(size = 0.35, color = "black"),
axis.ticks = ggplot2::element_line(size = 0.35, color = "black"),
text = ggplot2::element_text(size = 12)
)
} else {
plot_GHI <-
ggplot2::ggplot(data = lghi_all) + ggplot2::geom_boxplot(
ggplot2::aes(x = as.factor(DM),
y = GHI),
lwd = 0.25,
fatten = 0.4,
outlier.size = 0.6
) + ggplot2::facet_wrap(~ method) + viridis::scale_fill_viridis(option = "C",
discrete = TRUE,
begin = 0.4) + ggplot2::xlab("DM") + ggplot2::ylab("GHI") + ggplot2::labs(fill = "Method") + ggplot2::geom_hline(
yintercept = c(0.1, 0.05),
size = 0.35,
linetype = "longdash",
colour = "blue"
) + ggplot2::geom_vline(
xintercept = c(3, 7),
size = 0.45,
linetype = "longdash",
colour = "green"
) + ggplot2::theme_classic() + ggplot2::theme(
panel.grid.major.y = ggplot2::element_line(size = 0.1, color = "red"),
axis.line = ggplot2::element_line(size = 0.35, color = "black"),
axis.ticks = ggplot2::element_line(size = 0.35, color = "black"),
text = ggplot2::element_text(size = 12)
)
}
plot_NPP <-
ggplot2::ggplot(data = lghi_all) + ggplot2::geom_count(ggplot2::aes(x = DM,
y = NPP),
colour = "red",
alpha = 0.7) + ggplot2::facet_wrap( ~ method) + ggplot2::xlab("DM") + ggplot2::ylab("NPP") + ggplot2::theme_classic() + ggplot2::theme(
panel.grid.major = ggplot2::element_line(size = 0.1, color = "red"),
axis.line = ggplot2::element_line(size = 0.35, color = "black"),
axis.ticks = ggplot2::element_line(size = 0.35, color = "black"),
text = ggplot2::element_text(size = 12)
)
plot_ENPP <-
ggplot2::ggplot(data = lghi_all) + ggplot2::geom_count(
ggplot2::aes(x = DM,
y = ENPP),
colour = "red",
alpha = 0.8,
shape = 1
) + ggplot2::facet_wrap( ~ method) + ggplot2::xlab("DM") + ggplot2::ylab("ENPP") + ggplot2::theme_classic() + ggplot2::theme(
panel.grid.major = ggplot2::element_line(size = 0.1, color = "red"),
axis.line = ggplot2::element_line(size = 0.35, color = "black"),
axis.ticks = ggplot2::element_line(size = 0.35, color = "black"),
text = ggplot2::element_text(size = 12)
)
# scatter_GHI <-
# ggplot2::ggplot(data = dplyr::filter(lghi_all, method %in% c("DH", "SL", "MSL", "H") )) + ggplot2::geom_jitter(ggplot2::aes(
# x = DM,
# y = GHI,
# color = factor(method)
# ), size = 1.5, alpha = 0.5, width = 0.2) + viridis::scale_color_viridis(option = "D", discrete = TRUE) + ggplot2::xlab("DM") + ggplot2::ylab("GHI") + ggplot2::labs(color = "Method") + ggplot2::geom_hline(yintercept = c(0.1), size = 0.35, linetype = "longdash", colour = "blue") + ggplot2::theme_classic() + ggplot2::theme(
# #increase size of axis lines
# axis.line = ggplot2::element_line(size =
# 0.35, color = "black"),
# axis.ticks = ggplot2::element_line(size =
# 0.35, color = "black"),
# #Adjust legend position to maximize space, use a vector of proportion
# #across the plot and up the plot where you want the legend.
# #You can also use "left", "right", "top", "bottom", for legends on t
# #he side of the plot
# legend.position = c(.85, .7),
# #increase the font size
# text = ggplot2::element_text(size =
# 12)
# ) + ggplot2::geom_line( ggplot2::aes(x = DM, y = GHI_predicted, color = factor(method)), size=1, alpha = 0.8)
# ----
out <-
list(plot_GHI = plot_GHI,
plot_NPP = plot_NPP,
plot_ENPP = plot_ENPP)
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.