# Code to run examples #
# Mou, Liu, and Xu 2023 #


#################
# installation #
#################

### From CRAN
# remove.packages('panelView') # remove old version
# install.packages("panelView") 

### From Github
# install.packages("remotes", repos = "http://cran.r-project.org") # if not already installed
# remotes::install_github("xuyiqing/panelView") 

library("ggplot2")
library("panelView")
data("panelView", package = "panelView")

# Set the path where the figures will be saved
fig_path <- "R_Figures"


####################################################
# Plotting treatment conditions and missing values #
####################################################


# Figure 1(a)
p1a <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                 data = turnout, index = c("abb", "year"), 
                 xlab = "Year", ylab = "State")

ggsave(file.path(fig_path, "figure1a.png"), p1a)

# Figure 1(b)
p1b <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                 data = turnout, index = c("abb", "year"), 
                 xlab = "Year", ylab = "State", 
                 by.timing = TRUE, legend.labs = c("No EDR", "EDR"), 
                 background = "white", 
                 cex.main = 20, cex.axis= 8, 
                 cex.lab = 12, cex.legend = 12)

ggsave(file.path(fig_path, "figure1b.png"), p1b)

# Figure 1(c)
p1c <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                 data = turnout, index = c("abb", "year"), 
                 xlab = "Year", ylab = "State", 
                 pre.post = TRUE)

ggsave(file.path(fig_path, "figure1c.png"), p1c)

panelview(D = "policy_edr",
          data = turnout, index = c("abb", "year"),
          xlab = "Year", ylab = "State",
          pre.post = TRUE,
          main = "EDR Reform", axis.lab = "time")

# Figure 2 (with time gaps)
turnout2 <- turnout[!(turnout$year=="1924" | turnout$year=="1928" 
  | turnout$year == "1940"),]

p2 <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                data = turnout2, index = c("abb", "year"), 
                type = "treat", leave.gap = TRUE)

ggsave(file.path(fig_path, "figure2.png"), p2)

## with Treatment reversal #
# Figure 3 
p3 <- panelview(Capacity ~ demo + lnpop + lngdp, data = capacity, 
                index = c("ccode", "year"), axis.lab.gap = c(2,0), 
                main = "Democracy and State Capacity", show.id = 1:100, 
                axis.lab = "time", report.missing = TRUE)

ggsave(file.path(fig_path, "figure3.png"), p3)

## high-dimensional data; plot unique treatment histories
# Figure 4(a)
p4a <- panelview(Capacity ~ demo + lnpop + lngdp, data = capacity, 
                 index = c("ccode", "year"), axis.lab.gap = c(2,0), axis.lab.angle = 90,
                 collapse.history = TRUE) # sorted by cohort size

ggsave(file.path(fig_path, "figure4a.png"), p4a)

# Figure 4(b)
p4b <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                 data = turnout, index = c("abb", "year"), axis.lab.angle = 90,
                 by.timing = TRUE, collapse.history = TRUE) # sorted by treatment timing

ggsave(file.path(fig_path, "figure4b.png"), p4b)

## sort by cohort size
panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
          data = turnout, index = c("abb", "year"),
          by.timing = FALSE, collapse.history = TRUE) 

## More than two treatment conditions #
# Figure 5(a)
demo2 <- rep(0, dim(capacity)[1])
demo2[which(capacity$polity2 < -0.5)] <- -1
demo2[which(capacity$polity2 > 0.5)] <- 1
capacity$demo2 <- demo2

p5a <- panelview(Capacity ~ demo2 + lngdp, 
                 data = capacity, index = c("ccode", "year"), 
                 axis.lab.gap = c(2, 10), main = "Regime Type")

ggsave(file.path(fig_path, "figure5a.png"), p5a)

# Figure 5(b) 
p5b <- panelview(Capacity ~ polity2 + lngdp, 
                 data = capacity, index = c("ccode", "year"), 
                 axis.lab.gap = c(2, 10), main = "Regime Type", 
                 gridOff = TRUE)

ggsave(file.path(fig_path, "figure5b.png"), p5b)

## Showing missing values only 
# Figure 6
library(dplyr)
p6 <- capacity %>% panelview(Capacity ~ 1, index = c("ccode", "year"), 
                             axis.lab="off", type = "missing")

ggsave(file.path(fig_path, "figure6.png"), p6)

# Or
capacity %>% panelview(Y = "Capacity", index = c("ccode", "year"), 
                       axis.lab.gap = c(2, 10), type = "missing")

# additionally, plot aggregate treatment adoption
tr.avg <- aggregate(turnout$policy_edr, by=list(turnout$year), FUN=mean)
plot(tr.avg, type = "l", xlab = "Year", ylab = "EDR", ylim = c(0, 0.2))


##########################################################
# Plotting the temporal dynamics of an (outcome) variable
##########################################################

## Continuous outcome
# Figure 7(a)
p7a <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                 data = turnout, index = c("abb", "year"), 
                 type = "outcome", main = "EDR Reform and Turnout", 
                 legend.labs = c("Control States", "Treated States (before EDR)", 
                                 "Treated States (after EDR)"))

ggsave(file.path(fig_path, "figure7a.png"), p7a)

# Figure 7(b)
p7b <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor,
                 data = turnout, index = c("abb", "year"),
                 type = "outcome", main = "", by.group = TRUE, 
                 cex.main = 20, cex.main.sub = 15)

ggsave(file.path(fig_path, "figure7b.png"), p7b)

# Figure 8(a)
p8a <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                 data = turnout, index = c("abb", "year"), 
                 type = "outcome", main = "EDR Reform and Turnout (AL, AR, CT)", 
                 id = c("AL", "AR", "CT"))

ggsave(file.path(fig_path, "figure8a.png"), p8a)

# Figure 8(b)
p8b <- panelview(turnout ~ 1, data = turnout, 
                 index = c("abb", "year"), type = "outcome", 
                 main = "Turnout", ylim = c(0, 100), 
                 xlab = "Year", ylab = "Turnout")

ggsave(file.path(fig_path, "figure8b.png"), p8b)

## Discrete outcome

# Set random seed for exactly reproducible jittering
set.seed(0)

# Figure 9(a)
p9a <- panelview(Y ~ D, data = simdata, 
                 index = c("id", "time"), by.group = FALSE, 
                 outcome.type = "discrete", 
                 type = "outcome", xlim = c(8, 15))

ggsave(file.path(fig_path, "figure9a.png"), p9a)

# Figure 9(b)
p9b <- panelview(Y ~ D, data = simdata, 
                 index = c("id", "time"), type = "outcome",  
                 outcome.type = "discrete", by.group = TRUE,
                 xlim = c(8, 15))

ggsave(file.path(fig_path, "figure9b.png"), p9b)

########################################
# Visualizing bivariate relationships #
########################################

## Plotting bivariate relationships in aggregate
# Figure 10(a)
p10a <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor,
                  data = turnout, index = c("abb", "year"), 
                  type = "bivariate", style = c("c", "b"), ylab = "Turnout")

ggsave(file.path(fig_path, "figure10a.png"), p10a)

# Figure 10(b)
p10b <- panelview(Capacity ~ polity2, data = capacity, index = c("country", "year"), 
                  ylim = list(c(-1,0.5), c(-0.2, 0.4)), type = "bivar")

ggsave(file.path(fig_path, "figure10b.png"), p10b)

## Plotting bivariate relationship by each unit
# Figure 11(a)
p11a <- panelview(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                  data = turnout, index = c("abb", "year"),
                  type = "bivar", by.unit = TRUE,
                  show.id = 10:21, ylab = "Turnout")

ggsave(file.path(fig_path, "figure11a.png"), p11a)

# Figure 11(b)
p11b <- panelview(Capacity ~ polity2, data = capacity,
                  index = c("country", "year"),
                  type = "bivar", by.unit = TRUE,
                  show.id = c(20:25, 40:45))

ggsave(file.path(fig_path, "figure11b.png"), p11b)

sessionInfo()