---
title: "SIKAP JEAS"
output: 
  html_document:
    toc: true
    toc_float: true
    toc_collapsed: false
    toc_depth: 3
    number_sections: true
    theme: lumen
---

**version: `r format(Sys.time(), '%d %B, %Y -- %H:%M:%S')`**

```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

source("./extras/packages.R")
library(summarytools)
library(ggridges)
library(nnet)

rm(list = ls())

load("sikap_w1_w58.Rdata")
dat <- haven::zap_labels(dat)

source("./extras/prepare_vars.R")


```


# Data Description


N per week:

```{r, echo=FALSE, message=FALSE, warning=FALSE}
table(dat$week)
dat$weekfac = factor(dat$week,
                     levels = c(min(dat$week):max(dat$week)))
```

Total N:

```{r, echo=FALSE, message=FALSE, warning=FALSE}
nrow(dat)
```

Unique Respondents:

```{r, echo=FALSE, message=FALSE, warning=FALSE}
length(unique(dat$pid))
```

Respondents by number of repeat participation:

```{r, echo=FALSE, message=FALSE, warning=FALSE}
t = table(dat$counter)
t - c(t[-1],0)
```

Heatmap of weeks

```{r, echo=FALSE}

vec_pid = dat$pid[dat$counter == 2]

df = dat
df = subset(df, pid %in% vec_pid & counter %in% c(1,2,3,4,5,6))
df = df[, c("pid", "week", "counter")]

d = df %>% group_by(pid) %>%
  summarise(week1 = mean(week[counter == 1]),
         week2 = mean(week[counter == 2], na.rm=TRUE),
         week3 = mean(week[counter == 3], na.rm=TRUE),
         week4 = mean(week[counter == 4], na.rm=TRUE),
         week5 = mean(week[counter == 5], na.rm=TRUE),
         week6 = mean(week[counter == 6], na.rm=TRUE),
         
         diff12 = week2 - week1,
         diff23 = week3 - week2,
         diff34 = week4 - week3,
         diff45 = week5 - week4,
         diff56 = week6 - week5)

table(d$diff12)
table(d$diff23)
table(d$diff34)
table(d$diff45)
table(d$diff56)

# View(d[which(d$diff12 <9 | d$diff23 <9),])

d_counts <- list(
  n12 = d %>% count(week1, week2),
  n23 = d %>% count(week2, week3),
  n34 = d %>% count(week3, week4),
  n45 = d %>% count(week4, week5),
  n56 = d %>% count(week5, week6)
)

d_counts$n12 = subset(d_counts$n12, !is.na(week1) & !is.na(week2))
d_counts$n23 = subset(d_counts$n23, !is.na(week2) & !is.na(week3))
d_counts$n34 = subset(d_counts$n34, !is.na(week3) & !is.na(week4))
d_counts$n45 = subset(d_counts$n45, !is.na(week4) & !is.na(week5))
d_counts$n56 = subset(d_counts$n56, !is.na(week5) & !is.na(week6))

# Heatmap
p <- ggplot(d_counts$n12, aes(x = week1, y = week2, fill = n)) +
  geom_tile() +
  scale_x_continuous(name = "First Participation",
                     limits = c(1,46),
                     breaks = seq(1,46, 3),
                     labels = seq(1,46, 3)) +
  scale_y_continuous(name = "Second Participation",
                     limits = c(1,58),
                     breaks = seq(1,58, 3),
                     labels = seq(1,58, 3)) +
  scale_fill_gradient(low = "white", high = "black") +
  labs(x = "Week 1", y = "Week 2", fill = "Count") +
  theme_minimal()

p
```

```{r, echo=TRUE}

# d_counts$n12[which(d_counts$n12$n == max(d_counts$n12$n)), c("week1", "week2", "n")]
# d_counts$n12[which(d_counts$n12$n == min(d_counts$n12$n)), c("week1", "week2", "n")]
# 
# d_counts$n23[which(d_counts$n23$n == max(d_counts$n23$n)), c("week2", "week3", "n")]
# d_counts$n23[which(d_counts$n23$n == min(d_counts$n23$n)), c("week2", "week3", "n")]
# 
# d_counts$n34[which(d_counts$n34$n == max(d_counts$n34$n)), c("week3", "week4", "n")]
# d_counts$n34[which(d_counts$n34$n == min(d_counts$n34$n)), c("week4", "week4", "n")]
# 
# d_counts$n45[which(d_counts$n45$n == max(d_counts$n45$n)), c("week4", "week5", "n")]
# d_counts$n45[which(d_counts$n45$n == min(d_counts$n45$n)), c("week4", "week5", "n")]
# 
# d_counts$n56[which(d_counts$n56$n == max(d_counts$n56$n)), c("week5", "week6", "n")]
# d_counts$n56[which(d_counts$n56$n == min(d_counts$n56$n)), c("week6", "week6", "n")]



### combining week combinations

colnames(d_counts$n12) = colnames(d_counts$n23) = colnames(d_counts$n34) = colnames(d_counts$n45) = colnames(d_counts$n56) = c("w1", "w2", "n")

d_comb = bind_rows(d_counts$n12, d_counts$n23, d_counts$n34, d_counts$n45, d_counts$n56)

d_comb = d_comb %>% group_by(w1, w2) %>%
            dplyr::summarise(n = sum(n))


d_comb[which(d_comb$n == max(d_comb$n)), c("w1", "w2", "n")]
d_comb[which(d_comb$n == min(d_comb$n)), c("w1", "w2", "n")]

```


# Descriptive

## Demographic

```{r, echo=FALSE}

df = dat %>%
  group_by(weeklab) %>%
  dplyr::summarise(n = n(),
                   pct_female = round(mean(female)*100,2),
                   age_median = median(age),
                   age_18_24 = round(mean(age %in% 18:24)*100,2),
                   age_25_34 = round(mean(age %in% 25:34)*100,2),
                   age_35_44 = round(mean(age %in% 35:44)*100,2),
                   age_45_54 = round(mean(age %in% 45:54)*100,2),
                   age_55 = round(mean(age >= 55)*100,2),
                   pct_muslim = round(mean(religionlab3pt == "Islam")*100,2),
                   pct_christian = round(mean(religionlab3pt == "Christianity")*100,2),
                   pct_others = round(mean(religionlab3pt == "Others")*100,2),
                   pct_edu_less_hs = round(mean(edulab3pt == "Lower than HS")*100,2),
                   pct_edu_hs = round(mean(edulab3pt == "High School")*100,2),
                   pct_edu_more_hs = round(mean(edulab3pt == "College or Higher")*100,2),
                   pct_region_java = round(mean(region == "Java & Bali")*100,2),
                   pct_region_sumatera = round(mean(region == "Sumatera")*100,2),
                   pct_region_eastern = round(mean(region == "Central & Eastern")*100,2),
                   pct_boomers = round(mean(age_group == "Boomers (born 1944-1964)", na.rm=TRUE)*100,2),
                   pct_genx = round(mean(age_group == "Gen-X (born 1965-1980)", na.rm=TRUE)*100,2),
                   pct_millenials = round(mean(age_group == "Millennials (born 1981-1996)", na.rm=TRUE)*100,2),
                   pct_genz = round(mean(age_group == "Gen-Z (born after 1996)", na.rm=TRUE)*100,2),
                   pct_anies = round(mean(polpref_votepres == 1, na.rm=TRUE)*100, 2),
                   pct_prabowo = round(mean(polpref_votepres == 2, na.rm=TRUE)*100, 2),
                   pct_ganjar = round(mean(polpref_votepres == 3, na.rm=TRUE)*100, 2)
                   )

last = dat %>%
  dplyr::summarise(weeklab = "Total",
                   n = n(),
                   pct_female = round(mean(female)*100,2),
                   age_median = median(age),
                   age_18_24 = round(mean(age %in% 18:24)*100,2),
                   age_25_34 = round(mean(age %in% 25:34)*100,2),
                   age_35_44 = round(mean(age %in% 35:44)*100,2),
                   age_45_54 = round(mean(age %in% 45:54)*100,2),
                   age_55 = round(mean(age >= 55)*100,2),
                   pct_muslim = round(mean(religionlab3pt == "Islam")*100,2),
                   pct_christian = round(mean(religionlab3pt == "Christianity")*100,2),
                   pct_others = round(mean(religionlab3pt == "Others")*100,2),
                   pct_edu_less_hs = round(mean(edulab3pt == "Lower than HS")*100,2),
                   pct_edu_hs = round(mean(edulab3pt == "High School")*100,2),
                   pct_edu_more_hs = round(mean(edulab3pt == "College or Higher")*100,2),
                   pct_region_java = round(mean(region == "Java & Bali")*100,2),
                   pct_region_sumatera = round(mean(region == "Sumatera")*100,2),
                   pct_region_eastern = round(mean(region == "Central & Eastern")*100,2),
                   pct_boomers = round(mean(age_group == "Boomers (born 1944-1964)", na.rm=TRUE)*100,2),
                   pct_genx = round(mean(age_group == "Gen-X (born 1965-1980)", na.rm=TRUE)*100,2),
                   pct_millenials = round(mean(age_group == "Millennials (born 1981-1996)", na.rm=TRUE)*100,2),
                   pct_genz = round(mean(age_group == "Gen-Z (born after 1996)", na.rm=TRUE)*100,2),
                   pct_anies = round(mean(polpref_votepres == 1, na.rm=TRUE)*100, 2),
                   pct_prabowo = round(mean(polpref_votepres == 2, na.rm=TRUE)*100, 2),
                   pct_ganjar = round(mean(polpref_votepres == 3, na.rm=TRUE)*100, 2)
                   )
df = data.frame(rbind(df, last))

rownames(df) = NULL

df[,c("weeklab", "n", "pct_female", "pct_muslim", "pct_christian", "pct_others")]

df[,c("weeklab", "age_median", "age_18_24", "age_25_34", "age_35_44", "age_45_54", "age_55")]

df[, c("weeklab", "pct_edu_less_hs", "pct_edu_hs", "pct_edu_more_hs")]

df[,c("weeklab", "pct_region_java", "pct_region_sumatera", "pct_region_eastern")]

df[,c("weeklab", "pct_boomers", "pct_genx", "pct_millenials", "pct_genz")]

```

## Distribution of Respondents

```{r, echo=FALSE}

df_prov = dat %>%
  group_by(provincelab) %>%
  dplyr::summarise(n = n(),
                   prop = n/nrow(dat)*100)

map = gadm(country = "IDN", path = getwd(), level = 1)

# create nmprov
map$provincelab = map$NAME_1

### change province name in gdam to follow SIKAP
map$provincelab = car::recode(map$provincelab,
                              "'Aceh' = 'Nanggroe Aceh';
                              'Bangka Belitung' = 'Bangka-Belitung';
                              'Jakarta Raya' = 'DKI Jakarta';
                              'Nusa Tenggara Barat' = 'NTB';
                              'Nusa Tenggara Timur' = 'NTT';
                              'Yogyakarta' = 'D.I. Yogyakarta'")


### merge gadm and ifls
map = merge(map, df_prov,
             by = c("provincelab"),
                all.x = TRUE, all.y = FALSE)


p_dist = ggplot(map) + 
  geom_spatvector(aes(fill = prop)) +
  scale_fill_distiller(name = "Percent", 
                       palette = "YlGn",
                       breaks = pretty_breaks(n = 10),
                       trans = "reverse") +
  labs(subtitle = "Geographic Distribution of SIKAP Respondents") +
  theme_bw()

```

## Distribution of Respondents vs BPS

https://sulut.bps.go.id/indicator/12/958/1/jumlah-penduduk-menurut-provinsi-di-indonesia.html

```{r, echo=FALSE}

bps = read.csv("./extras/bps_penduduk.csv", header = TRUE)
bps = bps[, c("provincelab", "X2022")]
colnames(bps) = c("provincelab", "n_bps")

bps$prop_bps = bps$n_bps/sum(bps$n_bps)*100

### change province name in bps to follow SIKAP (and matched with gdam)
bps$provincelab = car::recode(bps$provincelab,
                              "'Aceh' = 'Nanggroe Aceh';
                              'Kep. Riau' = 'Kepulauan Riau';
                              'Kep. Bangka Belitung' = 'Bangka-Belitung';
                              'Nusa Tenggara Barat' = 'NTB';
                              'Nusa Tenggara Timur' = 'NTT';
                              'DI Yogyakarta' = 'D.I. Yogyakarta'")


### merge gadm and ifls
map = merge(map, bps,
             by = c("provincelab"),
                all.x = TRUE, all.y = FALSE)

### difference between SIKAP percentage and BPS percentage
map$diff = map$prop - map$prop_bps

p_diff = ggplot(map) + 
  geom_spatvector(aes(fill = diff)) +
  scale_fill_distiller(name = "Percentage\nPoint", 
                       palette = "YlGn",
                       breaks = pretty_breaks(n = 10),
                       trans = "reverse") +
  labs(subtitle = "Discrepancy (SIKAP Data Minus Official Statistic)") +
  theme_bw()


```

```{r, echo=FALSE}
p <- p_diff
p

ggsave(filename = "./output/figure_1.pdf", plot = p, width = 12, height = 10)
ggsave(filename = "./output/figure_1.tiff", plot = p, width = 12, height = 10, dpi = 300)

```

# Vote Share

```{r, echo=FALSE}


# prop.table(table(dat$polpref_votepreslab3pt[dat$week == 11]))

d = read.xlsx("./extras/sikap_jeas.xlsx", sheet = 1)
d$source = factor(d$source,
                  levels = c("SIKAP (Week 11)",
                             "Indikator's January Survey",
                             "Indikator's Exit Poll",
                             "Official Result"))
d$candidate = factor(d$candidate,
                     levels = c("Anies Baswedan",
                                "Prabowo Subianto",
                                "Ganjar Pranowo"))


p <- ggplot(d, aes(x = candidate, y = vote, group = source)) +
  geom_bar(aes(fill = source), stat = "identity", position = position_dodge(.9)) +
  geom_text(aes(label = round(vote,2)), vjust = -0.2, size = 3, position = position_dodge(.9)) +
  scale_x_discrete(name = "Candidate") +
  scale_y_continuous(name = "Vote Share",
                     limits = c(0,62),
                     breaks = seq(0,100,10),
                     labels = seq(0,100,10)) +
  theme_bw() +
  scale_fill_manual(values = c("#d6d6d6", "#ababab", "#757575", "#454545")) +
  labs(fill = NULL) +
  theme(legend.position = "bottom")

p

ggsave(p, filename = "./output/figure_2.pdf", width = 7, height = 4)
ggsave(p, filename = "./output/figure_2.tiff", width = 7, height = 4, dpi = 300)
```


# Approval and Favorability

## Approval

```{r, echo=FALSE}

m = lm(polpref_jokowisatis2pt  ~ 
         femalefac + age_group + edulab3pt + religionlab3pt + region + javanese + weeklab,
       data = dat)

df = ggeffect(m, terms = "weeklab")
df$predicted = df$predicted*100
df$conf.low = df$conf.low*100
df$conf.high = df$conf.high*100

p_approval <- ggplot(df, aes(x = x, y = predicted)) +
  geom_line(group = 1) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), group = 1, alpha = .1) +
  scale_x_discrete(name = "Survey Week") +
  scale_y_continuous(name = "Percent Respondents",
                     limits = c(50, 100),
                     breaks = seq(-15,100,5),
                     labels = seq(-15,100,5)) +
  geom_hline(yintercept = 50) +
  geom_vline(xintercept = 26, linetype = 4) +
  geom_vline(xintercept = 27, linetype = 4) +
  geom_label(aes(x = 25, y = 55), label = "MK's Ruling\nDPR's Defiance\nPublic Protests") + 
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(subtitle = "Satisfied or Very Satisfied with Jokowi's Performance")
p_approval

df
# ggsave("approval_trend.pdf", p_approval, width = 8, height = 6)
```

## Favorability

```{r, echo=FALSE}

dat$feelpol_president_jokowi2pt = car::recode(dat$feelpol_president_jokowi,
                                              "1:3 = 0;
                                              4:5 = 1")
dat$feelpol_president_prabowo2pt = car::recode(dat$feelpol_president_prabowo,
                                              "1:3 = 0;
                                              4:5 = 1")


m_jokowi = lm(feelpol_president_jokowi2pt  ~ 
         femalefac + age_group + edulab3pt + religionlab3pt + region + javanese + weeklab,
       data = dat)
df_jokowi = ggeffect(m_jokowi, terms = "weeklab")
df_jokowi$predicted = df_jokowi$predicted*100
df_jokowi$conf.low = df_jokowi$conf.low*100
df_jokowi$conf.high = df_jokowi$conf.high*100
df_jokowi$figure = "Jokowi"

m_prabowo = lm(feelpol_president_prabowo2pt  ~ 
         femalefac + age_group + edulab3pt + religionlab3pt + region + javanese + weeklab,
       data = dat)
df_prabowo = ggeffect(m_prabowo, terms = "weeklab")
df_prabowo$predicted = df_prabowo$predicted*100
df_prabowo$conf.low = df_prabowo$conf.low*100
df_prabowo$conf.high = df_prabowo$conf.high*100
df_prabowo$figure = "Prabowo"

df = bind_rows(df_jokowi, df_prabowo)

p_fav <- ggplot(df, aes(x = x, y = predicted, group = figure)) +
  geom_line(aes(linetype = figure)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, group = figure), alpha = .1) +
  scale_x_discrete(name = "Survey Week") +
  scale_y_continuous(name = "Percent Respondents",
                     limits = c(40, 100),
                     breaks = seq(-15,100,5),
                     labels = seq(-15,100,5)) +
  geom_hline(yintercept = 50) +
  geom_vline(xintercept = 10, linetype = 4) +
  geom_vline(xintercept = 11, linetype = 4) +
  geom_label(aes(x = 8, y = 45), label = "MK's Ruling\nDPR's Defiance\nPublic Protests") + 
  geom_vline(xintercept = 19, linetype = 4) +
  geom_vline(xintercept = 20, linetype = 4) +
  geom_label(aes(x = 18, y = 45), label = "Prabowo\nInauguration") + 
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position = "bottom") +
  labs(subtitle = "Positive or Very Positive Feelings toward Jokowi and Prabowo",
       linetype = "")
p_fav

df
# ggsave("favorability_trend.pdf", p_fav, width = 8, height = 6)

```


### Combined

```{r, echo=FALSE}

p <- p_approval + p_fav
ggsave(p, filename = "./output/figure_3.pdf", width = 12, height = 5)
ggsave(p, filename = "./output/figure_3.tiff", width = 12, height = 5, dpi = 300)

```


## By Java Island

### Approval

```{r, echo=FALSE}

m = lm(polpref_jokowisatis2pt  ~ 
         femalefac + age_group + edulab3pt + religionlab3pt + javanese + weeklab*javaislandlab,
       data = dat)

d = data.frame(ggeffect(m, terms = c("weeklab", "javaislandlab")))[,c(1,2,3,6)]
colnames(d) = c("week", "est", "se", "javaisland")
d$javaisland = factor(d$javaisland, levels = c("Java Island", "Outside Java Island"))

d$est = round(d$est*100,2)
d$se = round(d$se*100,2)
d$lo95 = d$est - qnorm(.975)*d$se
d$hi95 = d$est + qnorm(.975)*d$se

p_java1 <- ggplot(d, aes(x = week, y = est, group = javaisland)) +
  geom_line(aes(linetype = javaisland)) +
  # geom_ribbon(aes(ymin = lo95, ymax = hi95, group = javaisland), alpha = .1) + 
   scale_x_discrete(name = "Survey Week") +
  scale_y_continuous(name = "Percent Respondents",
                     limits = c(50, 100),
                     breaks = seq(-15,100,5),
                     labels = seq(-15,100,5)) +
  geom_hline(yintercept = 50) +
  geom_vline(xintercept = 26, linetype = 4) +
  geom_vline(xintercept = 27, linetype = 4) +
  geom_label(aes(x = 25, y = 55), label = "MK's Ruling\nDPR's Defiance\nPublic Protests") + 
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position = c(.15, .20)) +
  labs(subtitle = "Satisfied or Very Satisfied with Jokowi's Performance",
       linetype = "Residence")

p_java1

# ggsave("approval_byjava.pdf", p_java1, width = 8, height = 6)
```



### Favorability

```{r, echo=FALSE}

m = lm(feelpol_president_jokowi2pt  ~ 
         femalefac + age_group + edulab3pt + religionlab3pt + javanese + weeklab*javaislandlab,
       data = dat)

d = data.frame(ggeffect(m, terms = c("weeklab", "javaislandlab")))[,c(1,2,3,6)]
colnames(d) = c("week", "est", "se", "javaisland")
d$javaisland = factor(d$javaisland, levels = c("Java Island", "Outside Java Island"))


d$est = round(d$est*100,2)
d$se = round(d$se*100,2)
d$lo95 = d$est - qnorm(.975)*d$se
d$hi95 = d$est + qnorm(.975)*d$se

p_java2 <- ggplot(d, aes(x = week, y = est, group = javaisland)) +
  geom_line(aes(linetype = javaisland)) +
  # geom_ribbon(aes(ymin = lo95, ymax = hi95, group = javaisland), alpha = .1) + 
   scale_x_discrete(name = "Survey Week") +
  scale_y_continuous(name = "Percent Respondents",
                     limits = c(40, 100),
                     breaks = seq(-15,100,5),
                     labels = seq(-15,100,5)) +
  geom_hline(yintercept = 50) +
  geom_vline(xintercept = 10, linetype = 4) +
  geom_vline(xintercept = 11, linetype = 4) +
  geom_label(aes(x = 7, y = 45), label = "MK's Ruling\nDPR's Defiance\nPublic Protests") + 
  geom_vline(xintercept = 19, linetype = 4) +
  geom_vline(xintercept = 20, linetype = 4) +
  geom_label(aes(x = 18, y = 45), label = "Prabowo\nInauguration") + 
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position = c(.85, .15)) +
  labs(subtitle = "Positive or Very Positive Feeling toward Jokowi",
       linetype = "Residence")

p_java2

# ggsave("favorability_byjava.pdf", p_java2, width = 8, height = 6)
```



## By Media

### Approval

```{r, echo=FALSE}

dat$socmed3pt = car::recode(dat$infosource_socmed04,
                           "0:1 = '0 to 1';
                           2:3 = '2 to 3';
                           4 = '4'")
dat$socmed3pt = factor(dat$socmed3pt,
                       levels = c("0 to 1", "2 to 3", "4"))

m_socmed = lm(polpref_jokowisatis2pt  ~ 
         femalefac + age_group + edulab3pt + religionlab3pt + region + javanese + 
           weeklab*socmed3pt,
       data = dat)

d_socmed = data.frame(ggeffect(m_socmed, terms = c("weeklab", "socmed3pt")))[,c(1,2,3,6)]
colnames(d_socmed) = c("week", "est", "se", "daily")


d = d_socmed

d$est = round(d$est*100,2)
d$se = round(d$se*100,2)
d$lo95 = d$est - qnorm(.975)*d$se
d$hi95 = d$est + qnorm(.975)*d$se

p_media1 <- ggplot(d, aes(x = week, y = est, group = daily)) +
  geom_line(aes(linetype = daily)) +
  # geom_ribbon(aes(ymin = lo95, ymax = hi95, group = source), alpha = .1) +
   scale_x_discrete(name = "Survey Week") +
  scale_y_continuous(name = "Percent Respondents",
                     limits = c(50, 100),
                     breaks = seq(-15,100,5),
                     labels = seq(-15,100,5)) +
  geom_hline(yintercept = 50) +
  geom_vline(xintercept = 26, linetype = 4) +
  geom_vline(xintercept = 27, linetype = 4) +
  geom_label(aes(x = 25, y = 55), label = "MK's Ruling\nDPR's Defiance\nPublic Protests") + 
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position = c(.15, .20),
        legend.key.width = unit(1.5, "cm")) +
  labs(subtitle = "Satisfied or Very Satisfied with Jokowi's Performance",
       linetype = "Number of Platforms")

p_media1

```

### Favorability

```{r, echo=FALSE}

m_socmed = lm(feelpol_president_jokowi2pt  ~ 
         femalefac + age_group + edulab3pt + religionlab3pt + region + javanese + 
           weeklab*socmed3pt,
       data = dat)


d_socmed = data.frame(ggeffect(m_socmed, terms = c("weeklab", "socmed3pt")))[,c(1,2,3,6)]
colnames(d_socmed) = c("week", "est", "se", "daily")
d_socmed$source = "Social Media"


d = d_socmed
d$est = round(d$est*100,2)
d$se = round(d$se*100,2)
d$lo95 = d$est - qnorm(.975)*d$se
d$hi95 = d$est + qnorm(.975)*d$se

p_media2 <- ggplot(d, aes(x = week, y = est, group = daily)) +
  geom_line(aes(linetype = daily)) +
  # geom_ribbon(aes(ymin = lo95, ymax = hi95, group = source), alpha = .1) +
   scale_x_discrete(name = "Survey Week") +
  scale_y_continuous(name = "Percent Respondents",
                     limits = c(30, 100),
                     breaks = seq(-15,100,5),
                     labels = seq(-15,100,5)) +
  geom_hline(yintercept = 50) +
  geom_vline(xintercept = 10, linetype = 4) +
  geom_vline(xintercept = 11, linetype = 4) +
  geom_label(aes(x = 7, y = 35), label = "MK's Ruling\nDPR's Defiance\nPublic Protests") + 
  geom_vline(xintercept = 19, linetype = 4) +
  geom_vline(xintercept = 20, linetype = 4) +
  geom_label(aes(x = 18, y = 35), label = "Prabowo\nInauguration") + 
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position = c(.85, .20),
        legend.key.width = unit(1.5, "cm")) +
  labs(subtitle = "Positive or Very Positive Feeling toward Jokowi",
       linetype = "Number of Platforms")

p_media2

```



```{r, echo=FALSE}
p <- ((p_java1 + p_java2) + plot_layout(guides = "collect") & theme(legend.position = 'bottom',
                                                              legend.text=element_text(size=12))) /
  ((p_media1 + p_media2) + plot_layout(guides = "collect") & theme(legend.position = 'bottom',
                                                              legend.text=element_text(size=12)))

ggsave(p, filename = "./output/figure_4.pdf", width = 12, height = 9)
ggsave(p, filename = "./output/figure_4.tiff", width = 12, height = 9, dpi = 300)

```



# Democratic Suitability

### General

```{r, echo=FALSE}

dat <- dat %>%
  mutate(
    start_date = str_sub(startdate, 1, 10),
    date_plot  = as.Date(start_date, "%Y-%m-%d")
  ) %>%
  group_by(week) %>%
  mutate(min_date = min(date_plot)) %>%
  ungroup() %>%                                   
  mutate(min_date_f = factor(min_date))

m <- lm(
  democ_suitable2pt ~ femalefac + age_group + edulab3pt + religionlab3pt +
    region + javanese + min_date_f,
  data = dat
)

df <- as.data.frame(ggeffect(m, terms = "min_date_f")) %>%
  mutate(
    x_date   = as.Date(as.character(x)),  # turn factor/chr like "2024-01-01" into Date
    predicted = predicted * 100,
    conf.low  = conf.low  * 100,
    conf.high = conf.high * 100
  ) %>%
  arrange(x_date)

p_dem_approval <- ggplot(df, aes(x = x_date, y = predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1) +
  scale_x_date(
    name = "Month",
    date_breaks = "1 month",
    date_labels = "%b %Y",
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  scale_y_continuous(
    name = "Percent Respondents",
    limits = c(50, 100),
    breaks = seq(50, 100, 5),
    labels = seq(50, 100, 5)
  ) +
  geom_vline(xintercept = as.Date("2024-02-14"), linetype = 4) +
  annotate("label", x = as.Date("2024-02-14"), y = 55,
           label = "General Election\n(February 14, 2024)") +
  labs(subtitle = "Democracy Suitable for Indonesia (%)") +
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    # remove the next line if you actually want the "Month" axis title to show
    axis.title.x = element_blank(),
    strip.background = element_blank(),
    legend.title = element_blank(),
    axis.line = element_line(colour = "black"),
    panel.border = element_blank(),
    legend.position = "bottom"
  )

p_dem_approval

# ggsave("approval_trend.pdf", p_approval, width = 8, height = 6)
```

### By Vote

```{r, echo=FALSE}
# make the factor levels once across ALL data (you already do this)
dat <- dat %>%
  mutate(
    polpref_pres_coll = coalesce(polpref_votepres, polpref_votepresb),
    polpref_pres_coll = case_when(
      polpref_pres_coll == 2 ~ "Prabowo Voter",
      polpref_pres_coll %in% c(1, 3) ~ "Non-Prabowo Voter",
      TRUE ~ NA_character_
    ),
    start_date = str_sub(startdate, 1, 10),
    date_plot  = as.Date(start_date, "%Y-%m-%d")
  ) %>%
  group_by(week) %>%
  mutate(min_date = min(date_plot)) %>%
  ungroup() %>%
  mutate(min_date_f = factor(min_date, levels = sort(unique(min_date))))

# models
m_prab_voter <- lm(
  democ_suitable2pt ~ femalefac + age_group + edulab3pt + religionlab3pt +
    region + javanese + min_date_f,
  data = dat %>% filter(polpref_pres_coll == "Prabowo Voter")
)

m_non_prab <- lm(
  democ_suitable2pt ~ femalefac + age_group + edulab3pt + religionlab3pt +
    region + javanese + min_date_f,
  data = dat %>% filter(polpref_pres_coll == "Non-Prabowo Voter")
)

# effects → data frames with a real Date column
df_prab <- as.data.frame(ggeffect(m_prab_voter, terms = "min_date_f")) %>%
  mutate(
    x_date   = as.Date(as.character(x)),
    predicted = predicted * 100,
    conf.low  = conf.low  * 100,
    conf.high = conf.high * 100,
    figure = "Prabowo Voter"
  )

df_non_prab <- as.data.frame(ggeffect(m_non_prab, terms = "min_date_f")) %>%
  mutate(
    x_date   = as.Date(as.character(x)),
    predicted = predicted * 100,
    conf.low  = conf.low  * 100,
    conf.high = conf.high * 100,
    figure = "Non-Prabowo Voter"
  )

df <- bind_rows(df_prab, df_non_prab) %>% arrange(x_date)

# plot
p_suit_split <- ggplot(df, aes(x = x_date, y = predicted)) +
  geom_line(aes(linetype = figure)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, group = figure), alpha = 0.12) +
  scale_x_date(
    name = "Month",
    date_breaks = "1 month",
    date_labels = "%b %Y",
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  scale_y_continuous(
    name = "Percent Respondents",
    limits = c(40, 100),
    breaks = seq(40, 100, 5),
    labels = seq(40, 100, 5)
  ) +
  geom_vline(xintercept = as.Date("2024-02-14"), linetype = 4) +
  annotate("label", x = as.Date("2024-02-14"), y = 42,
           label = "General Election\n(February 14, 2024)") +
  labs(subtitle = "Democracy Suitable for Indonesia (%)") +
  theme_bw() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    axis.title.x = element_blank(),
    strip.background = element_blank(),
    legend.title = element_blank(),
    axis.line = element_line(colour = "black"),
    panel.border = element_blank(),
    legend.position = "bottom"
  )

p_suit_split
# ggsave("favorability_trend.pdf", p_fav, width = 8, height = 6)

```


### Combined

```{r, echo=FALSE}

p <- p_dem_approval + p_suit_split

ggsave(p, filename = "./output/figure_5.pdf", width = 12, height = 5)
ggsave(p, filename = "./output/figure_5.tiff", width = 12, height = 5, dpi = 300)
```

