getwd()
## [1] "/Users/idlhy/Library/CloudStorage/OneDrive-개인/R FILE"
setwd("/Users/idlhy/Library/CloudStorage/OneDrive-개인/바탕 화면/Research/Puberty Delinquency(HRB)")
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(readr)
library(haven)
library(readxl)

# 1. Intercept
intercept <- read_excel("./4. sas print files/220929_female_w2_w6_only drinking.xlsx", sheet = "intercept_prevalence")
colors <- c("late" = "#0000CC",
            "early" = "#CC0000",
            "on time" = "#1B9E77")
linetypes <- c("late" = 5,
               "early" = 1,
               "on time" = 6)

ggplot(intercept) +
  aes(x = Grade) +
  geom_line(aes(y = early, color = "early", lty = "early"), size = 1) +
  geom_line(aes(y = normal, color = "on time", lty = "on time"), size = 1) +
  geom_line(aes(y = late, color = "late", lty = "late"), size = 1) +
  geom_ribbon(aes(ymin = e_l, ymax = e_u), color = "#CC0000", alpha = 0.05, lty = "dotted") +
  geom_ribbon(aes(ymin = n_l, ymax = n_u), color = "#1B9E77", alpha = 0.05, lty = "dotted") +
  geom_ribbon(aes(ymin = l_l, ymax = l_u), color = "#0000CC", alpha = 0.05, lty = "dotted") +
  theme_bw() +
  scale_color_manual(values = colors) +
  scale_linetype_manual(values = linetypes) + 
  labs(x = "Grade"
       , y = "Prevalence"
       , color = "Puberty Group"
       , lty = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))

## Only E and L
ggplot(intercept) +
  aes(x = Grade) +
  geom_line(aes(y = early, color = "early", lty = "early"), size = 1) +
  geom_line(aes(y = late, color = "late", lty = "late"), size = 1) +
  geom_ribbon(aes(ymin = e_l, ymax = e_u), color = "#CC0000", alpha = 0.1, lty = "dotted") +
  geom_ribbon(aes(ymin = l_l, ymax = l_u), color = "#0000CC", alpha = 0.1, lty = "dotted") +
  theme_bw() +
  scale_color_manual(values = colors
                     , limits = c("early", "late")) +
  labs(x = "Grade"
       , y = "Prevalence"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))

# 2. Peer Substance Abuse
PeerSubstanceUse <- read_excel("./4. sas print files/220929_female_w2_w6_only drinking.xlsx", sheet = "Peer Substance Use")

ggplot(PeerSubstanceUse) +
  aes(x = Grade) +
  geom_line(aes(y = early, color = "early", lty = "early"), size = 1) +
  geom_line(aes(y = normal, color = "on time", lty = "on time"), size = 1) +
  geom_line(aes(y = late, color = "late", lty = "late"), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = e_l, ymax = e_u), color = "#CC0000", alpha = 0.04, lty = "dotted") +
  geom_ribbon(aes(ymin = n_l, ymax = n_u), color = "#1B9E77", alpha = 0.04, lty = "dotted") +
  geom_ribbon(aes(ymin = l_l, ymax = l_u), color = "#0000CC", alpha = 0.04, lty = "dotted") +
  theme_bw() +
  scale_color_manual(values = colors
                     , limits = c("early", "late", "on time")) +
  labs(x = "Grade"
       , y = "Coefficient"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(-3.5, 3.5)

## E and L only
ggplot(PeerSubstanceUse) +
  aes(x = Grade) +
  geom_line(aes(y = early, color = "early", lty = "early"), size = 1) +
  geom_line(aes(y = late, color = "late", lty = "late"), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = e_l, ymax = e_u), color = "#CC0000", alpha = 0.1, lty = "dotted") +
  geom_ribbon(aes(ymin = l_l, ymax = l_u), color = "#0000CC", alpha = 0.1, lty = "dotted") +
  theme_bw() +
  scale_color_manual(values = colors
                     , limits = c("early", "late")) +
  labs(x = "Grade"
       , y = "Coefficient"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16))) + 
  ylim(-3.5, 3.5)

#3. Parent Attachment
ParentAttachment <- read_excel("./4. sas print files/220929_female_w2_w6_only drinking.xlsx", sheet = "Parent Attachment")
ggplot(ParentAttachment) +
  aes(x = Grade) +
  geom_line(aes(y = early, color = "early", lty = "early"), size = 1) +
  geom_line(aes(y = normal, color = "on time", lty = "on time"), size = 1) +
  geom_line(aes(y = late, color = "late", lty = "late"), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = e_l, ymax = e_u), color = "#CC0000", alpha = 0.04, lty = "dotted") +
  geom_ribbon(aes(ymin = n_l, ymax = n_u), color = "#1B9E77", alpha = 0.04, lty = "dotted") +
  geom_ribbon(aes(ymin = l_l, ymax = l_u), color = "#0000CC", alpha = 0.04, lty = "dotted") +
  theme_bw() +
  scale_color_manual(values = colors
                     , limits = c("early", "late", "on time")) +
  labs(x = "Grade"
       , y = "Coefficient"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(-3.5, 3.5)

## E and L only
ggplot(ParentAttachment) +
  aes(x = Grade) +
  geom_line(aes(y = early, color = "early", lty = "early"), size = 1) +
  geom_line(aes(y = late, color = "late", lty = "late"), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = e_l, ymax = e_u), color = "#CC0000", alpha = 0.1, lty = "dotted") +
  geom_ribbon(aes(ymin = l_l, ymax = l_u), color = "#0000CC", alpha = 0.1, lty = "dotted") +
  theme_bw() +
  scale_color_manual(values = colors
                     , limits = c("early", "late")) +
  labs(x = "Grade"
       , y = "Coefficient"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(-3.5, 3.5)

# 4. Coeducation
Coeducation <- read_excel("./4. sas print files/220929_female_w2_w6_only drinking.xlsx", sheet = "Coeducation")
ggplot(Coeducation) +
  aes(x = Grade) +
  geom_line(aes(y = early, color = "early", lty = "early"), size = 1) +
  geom_line(aes(y = normal, color = "on time", lty = "on time"), size = 1) +
  geom_line(aes(y = late, color = "late", lty = "late"), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = e_l, ymax = e_u), color = "#CC0000", alpha = 0.04, lty = "dotted") +
  geom_ribbon(aes(ymin = n_l, ymax = n_u), color = "#1B9E77", alpha = 0.04, lty = "dotted") +
  geom_ribbon(aes(ymin = l_l, ymax = l_u), color = "#0000CC", alpha = 0.04, lty = "dotted") +
  theme_bw() +
  scale_color_manual(values = colors
                     , limits = c("early", "late", "on time")) +
  labs(x = "Grade"
       , y = "Coefficient"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(-3.5, 3.5)

## E and L only
ggplot(Coeducation) +
  aes(x = Grade) +
  geom_line(aes(y = early, color = "early", lty = "early"), size = 1) +
  geom_line(aes(y = late, color = "late", lty = "late"), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = e_l, ymax = e_u), color = "#CC0000", alpha = 0.1, lty = "dotted") +
  geom_ribbon(aes(ymin = l_l, ymax = l_u), color = "#0000CC", alpha = 0.1, lty = "dotted") +
  theme_bw() +
  scale_color_manual(values = colors
                     , limits = c("early", "late")) +
  labs(x = "Grade"
       , y = "Coefficient"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(-3.5, 3.5)

# 5. Overall Group
# Intercept
T_int <- read_excel("./4. sas print files/220929_female_w2_w6_only drinking.xlsx", sheet = "T_int")
ggplot(T_int) +
  aes(x = Grade) +
  geom_line(aes(y = total), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = t_l, ymax = t_u), color = "#CC0000", alpha = 0.04, lty = "dotted") +
  theme_bw() +
  labs(x = "Grade"
       , y = "Prevalence"
       , title = "Intercept"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(0, 0.5)

# Parent attachment
T_parent <- read_excel("./4. sas print files/220929_female_w2_w6_only drinking.xlsx", sheet = "T_parent")
ggplot(T_parent) +
  aes(x = Grade) +
  geom_line(aes(y = total), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = t_l, ymax = t_u), color = "#CC0000", alpha = 0.04, lty = "dotted") +
  theme_bw() +
  labs(x = "Grade"
       , y = "Coefficient"
       , title = "Parent Attachment"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(-3.5, 3.5)

# Peer Substance Use
T_peer <- read_excel("./4. sas print files/220929_female_w2_w6_only drinking.xlsx", sheet = "T_peer")
ggplot(T_peer) +
  aes(x = Grade) +
  geom_line(aes(y = total), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = t_l, ymax = t_u), color = "#CC0000", alpha = 0.04, lty = "dotted") +
  theme_bw() +
  labs(x = "Grade"
       , y = "Coefficient"
       , title = "Peer Substance Abuse"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(-3.5, 3.5)

# Coeducation
T_coedu <- read_excel("./4. sas print files/220929_female_w2_w6_only drinking.xlsx", sheet = "T_coedu")
ggplot(T_coedu) +
  aes(x = Grade) +
  geom_line(aes(y = total), size = 1) +
  geom_line(aes(y = 0), size = 1) +
  geom_ribbon(aes(ymin = t_l, ymax = t_u), color = "#CC0000", alpha = 0.04, lty = "dotted") +
  theme_bw() +
  labs(x = "Grade"
       , y = "Coefficient"
       , title = "Coeducation"
       , lty = "Puberty Group"
       , color = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (16)),
        legend.text = element_text(family = "Times New Roman", size = (16)),
        axis.title = element_text(family = "Times New Roman", size = (16)),
        axis.text = element_text(family = "Times New Roman", size = (16)))+ 
  ylim(-3.5, 3.5)

# 6. Puberty Group Distribution
library(haven)
library(ggridges)
data <- read_dta("./2. stata_cleaning data/1. current_using_data/KCYPS2010m1_w2_w6_female_220910.dta")
data$n_puberty_f_3 <- data$n_puberty_f_3 %>%
  as.character() %>%
  fct_recode(
    "On time" = "0",
    "Early" = "1",
    "Late" = "2")
colors2 <- c("Late" = "#0000CC",
             "Early" = "#CC0000",
             "On time" = "#1B9E77")
data %>%
  filter(!is.na(n_puberty_f_3)) %>%
  ggplot(aes(x = n_puberty_f_2, 
             y = factor(n_puberty_f_3, levels = c("Early", "On time", "Late")), 
             height = ..density.., 
             fill = n_puberty_f_3)) +
  geom_density_ridges(stat = "density", 
                      alpha = 0.7) +
  scale_fill_manual(values = colors2) +
  theme_bw() +
  labs(x = "Age (month)"
       , y = ""
       , title = ""
       , fill = "Puberty Group") +
  theme(legend.position = "bottom",
        plot.title = element_text(family = "Times New Roman", size = (16)),
        legend.title = element_text(family = "Times New Roman", size = (15)),
        legend.text = element_text(family = "Times New Roman", size = (15)),
        axis.title = element_text(family = "Times New Roman", size = (15)),
        axis.text = element_text(family = "Times New Roman", size = (15)))