domain <- df_fields %>%
    dplyr::select("domain", "id") %>%
    filter(!is.na(domain)) %>%
    mutate(domain2 = as.numeric(as.factor(domain))) %>%
    # make numeric of domain variable
group_by(id) %>%
    # number of publications domain is based domain as proportions
mutate(npubs_domain = n(), health_domain = sum(domain2 == 1)/npubs_domain, life_domain = sum(domain2 ==
    2)/npubs_domain, physical_domain = sum(domain2 == 3)/npubs_domain, social_domain = sum(domain2 ==
    4)/npubs_domain) %>%
    # only keep one record per id
ungroup() %>%
    mutate(id.dup = ifelse(duplicated(id), 1, 0)) %>%
    filter(id.dup == 0) %>%
    dplyr::select("id", "npubs_domain", "health_domain", "life_domain", "physical_domain", "social_domain")

# describe(domain)

df_ppf <- df_ppf %>%
    left_join(domain)

df_ppf$npubs_mid_a <- df_ppf$npubs_a - df_ppf$npubs_first_a - df_ppf$npubs_last_a

df_mid <- df_ppf %>%
    group_by(id) %>%
    summarize(npubs_mid_a_t = sum(npubs_mid_a))

df_ppf <- df_ppf %>%
    left_join(df_mid)

df_gender <- df_pubs %>%
    ungroup() %>%
    mutate(id.dup = ifelse(duplicated(id), 1, 0)) %>%
    filter(id.dup == 0) %>%
    dplyr::select("id", "gender_2", "firstname", "lastname") %>%
    mutate(gender2 = ifelse(is.na(firstname) & is.na(gender_2), 3, gender_2))

df_gender$gender2[is.na(df_gender$gender2)] <- 4
df_gender$gender2 <- as.factor(df_gender$gender2)
levels(df_gender$gender2) <- c("men", "women", "no first name", "missing")

# table(df_gender$gender_2, useNA = 'always') table(df_gender$gender2, useNA = 'always')

tdf_gender <- df_gender %>%
    filter(is.na(gender_2)) %>%
    dplyr::select("id", "gender2", "firstname", "lastname")

fshowdf(tdf_gender)

df_ppf <- df_ppf %>%
    left_join(df_gender)

df_ppf$time <- df_ppf$pub_year - df_ppf$phd_year

df_ppf$covid <- as.numeric(df_ppf$pub_year == 2020 | df_ppf$pub_year == 2021)
df_ppf$covid2020 <- as.numeric(df_ppf$pub_year == 2020)
df_ppf$covid2021 <- as.numeric(df_ppf$pub_year == 2021)

df_ppf$c2005 <- as.numeric(df_ppf$phd_year > 2005)
df_ppf$c2010 <- as.numeric(df_ppf$phd_year > 2010)
df_ppf$c2015 <- as.numeric(df_ppf$phd_year > 2015)

df_ppf$phd_year0 <- df_ppf$phd_year - 1990

df_ppf$stem <- ifelse((df_ppf$field == "Biological and Health Sciences" | df_ppf$field == "Physical and Mathematical Sciences" |
    df_ppf$field == "Engineering" | df_ppf$field == "Agricultural Sciences"), "STEM-M", ifelse(df_ppf$field ==
    "missing", "missing", "NON-STEM-M"))

df_ppf$ethnicity <- as.factor(df_ppf$ethnicity)
df_ppf$ethnicity <- relevel(df_ppf$ethnicity, 2)
# levels(df_ppf$ethnicity)

df_ppf$multidisc <- 1 - (df_ppf$health_domain^2 + df_ppf$life_domain^2 + df_ppf$physical_domain^2 + df_ppf$social_domain^2)

table(df_ppf$time, df_ppf$covid)

childdf <- data.frame(cohort = 1990:2018, child5 = NA, N = NA)

childdf$child5 = c(0, 0.022, 0.008, 0.014, 0.022, 0.024, 0, 0.008, 0.016, 0.037, 0.032, 0.044, 0.06,
    0.052, 0.083, 0.07, 0.16, 0.233, 0.271, 0.279, 0.352, 0.426, 0.461, 0.497, 0.53, 0.397, 0.5, 0.413,
    0.336)
childdf$N = c(30, 89, 119, 148, 136, 127, 115, 132, 186, 162, 158, 205, 233, 210, 254, 227, 256, 288,
    284, 337, 341, 319, 343, 296, 166, 174, 198, 172, 110)

# mean(childdf$N)
fit.lo <- loess(childdf$child5 ~ childdf$cohort, family = "symmetric", weigths = testdf$N, degree = 2,
    span = 0.6)
childdf$smooth5 <- fitted(fit.lo)

# fshowdf(childdf)

ggplot(data = childdf, aes(x = cohort, y = child5, group = 1)) + geom_line(aes(x = cohort, y = smooth5,
    group = 1)) + geom_point(color = "red", size = childdf$N/100) + labs(x = "PhD-cohort", y = "proportion with children under 5",
    caption = "Source: Phd Survey 2014 and 2019 linked to Register data of Statistics Netherlands \nNotes: Size of points is proportional to number of observations. Mean N per cohort is 200; Smooth line created with loess.") +
    scale_x_continuous(breaks = seq(1990, 2020, by = 5)) + theme(plot.caption = element_text(hjust = 0))

df_phd_des <- df_phd %>%
    dplyr::select("npubs_t", "npubs_a_t", "npubs_first_a_t", "npubs_mid_a_t", "npubs_last_a_t", "phd_year",
        "gender2", "c2005", "ethnicity", "uni", "npubs_domain", "health_domain", "life_domain", "physical_domain",
        "social_domain", "multidisc") %>%
    dummy_cols(select_columns = c("gender2", "ethnicity", "uni")) %>%
    dplyr::select(!c("gender2", "ethnicity", "uni")) %>%
    na.omit()

Table1 <- df_phd_des %>%
    describe()


Table1 <- as.data.frame(Table1)[, c(3, 4, 5, 8, 9)]
rownames(Table1) <- c("all output", "all articles", "first authored articles", "mid position authored articles",
    "last authored articles", "PhD-cohort", "PhD after 2005", "output per domain", "domain: Health",
    "domain: Life", "domain: Physical", "domain: Social", "multidisciplinarity", "gender: men", "gender: women",
    "gender: no first name", "gender: missing", "ethnicity: native Dutch", "ethnicity: 'traditional' minority",
    "ethnicity: other minority", "ethnicity: missing", "university: EUR", "university: LU", "university: RU",
    "university: RUG", "university: TUD", "university: TUE", "university: TI", "university: UM", "university: UT",
    "university: UU", "university: UvA", "university: VU", "university: WUR")
# fsave(Table1, 'Table1_raw')

Table1 <- Table1[c(1:7, 14:21, 9:12, 22:34), ]

Table1 <- cbind(rownames(Table1), Table1)

tableft <- flextable(Table1) %>%
    set_header_labels(values = c(" ", "mean/prop.", "sd", "median", "min", "max")) %>%
    # set_caption( 'Descriptive statistics at Phd-level (N = 8.303') %>%
width(j = 1, width = 2) %>%
    align(i = 1, align = "center", part = "header") %>%
    colformat_double(j = c("min", "max"), digits = c(0), big.mark = "") %>%
    colformat_double(j = c("mean", "sd", "median"), digits = c(3), big.mark = "") %>%
    colformat_double(i = c("PhD-cohort"), j = c("mean", "median"), digits = c(0), big.mark = "") %>%
    set_table_properties(layout = "autofit") %>%
    footnote(i = 6, j = 1, value = as_paragraph(c("The smallest PhD-cohort is 1992 and consists of 70 men and 23 women. For a complete overview of the size of each cohort see our replication website.")),
        ref_symbols = c("a"), part = "body")

tableft

domain_gender <- df_phd %>%
    dplyr::select("gender2", "health_domain", "life_domain", "physical_domain", "social_domain", "multidisc",
        "c2005") %>%
    na.omit() %>%
    describeBy(group = "gender2")

domain_men <- domain_gender[[1]][2:7, 3:4]
domain_women <- domain_gender[[2]][2:7, 3:4]

domain_gender <- cbind(domain_men, domain_women)


# fsave(domain_gender, 'Appendix2_raw')
table(df_phd$gender2)

A2 <- domain_gender

rownames(A2) <- c("domain: Health", "domain: Life", "domain: Physical", "domain: Social", "multidisciplinarity",
    "PhD after 2005")

A2 <- cbind(rownames(A2), A2)
colnames(A2) <- c(" t", "mean", "sd", "mean1", "sd1")

A2ft <- flextable(A2, col_keys = c(" t", "mean", "sd", "col1", "mean1", "sd1")) %>%
    width(j = "col1", width = 0.2) %>%
    empty_blanks() %>%
    set_header_labels(values = c(" ", "mean", "sd", "", "mean", "sd")) %>%
    add_header_row(values = c("", "men \n(N=4.740)", "", "women \n(N=3.488)"), colwidths = c(1, 2, 1,
        2)) %>%
    width(j = 1, width = 2) %>%
    # set_caption( 'Research domain by gender') %>%
width(j = 1, width = 2) %>%
    align(i = 1, align = "center", part = "header") %>%
    colformat_double(digits = c(3)) %>%
    set_table_properties(layout = "autofit")
A2ft

cor.test(df_ppf$covid, df_ppf$time)  #career age

cor.test(df_ppf$covid, df_ppf$pub_year)  #period

df_ppf2 <- df_ppf[df_ppf$time > 1 & df_ppf$time < 32, ]
cor.test(df_ppf2$covid, df_ppf2$time)  #career age
cor.test(df_ppf2$covid, df_ppf2$pub_year)  #period

df_ppf %>%
    group_by(pub_year, phd_year) %>%
    summarize(fr = sum(npubs, na.rm = TRUE), naut = n(), nauthors = length(unique(id))) -> data_npubs_aggyear

df_ppf %>%
    filter(gender2 == "men") %>%
    group_by(pub_year, phd_year) %>%
    summarize(fr = sum(npubs, na.rm = TRUE), naut = n(), nauthors = length(unique(id))) -> data_npubs_aggyear_m

df_ppf %>%
    filter(gender2 == "women") %>%
    group_by(pub_year, phd_year) %>%
    summarize(fr = sum(npubs, na.rm = TRUE), naut = n(), nauthors = length(unique(id))) -> data_npubs_aggyear_w


df_ppf %>%
    group_by(pub_year) %>%
    summarize(fr = sum(npubs, na.rm = TRUE), naut = n(), nauthors = length(unique(id))) -> data_npubs_aggyear2


Figure1 <- ggplot(data_npubs_aggyear, aes(x = pub_year, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + labs(color = "Phd Cohort", x = "year", y = "number of publications") + scale_y_continuous(limits = c(0,
    10)) + scale_x_continuous(limits = c(1990, 2024)) + theme_minimal() + theme(plot.background = element_rect(colour = "black",
    fill = NA, size = 1))
Figure1

Figure1_m <- ggplot(data_npubs_aggyear_m, aes(x = pub_year, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + labs(color = "Phd Cohort", x = "year", y = "number of publications") + scale_y_continuous(limits = c(0,
    10)) + scale_x_continuous(limits = c(1990, 2024)) + theme_minimal() + theme(plot.background = element_rect(colour = "black",
    fill = NA, size = 1))
Figure1_m

Figure1_w <- ggplot(data_npubs_aggyear_w, aes(x = pub_year, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + labs(color = "Phd Cohort", x = "year", y = "number of publications") + scale_y_continuous(limits = c(0,
    10)) + scale_x_continuous(limits = c(1990, 2024)) + theme_minimal() + theme(plot.background = element_rect(colour = "black",
    fill = NA, size = 1))
Figure1_w

Figure1 <- ggplot(data_npubs_aggyear, aes(x = pub_year, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + labs(color = "Phd Cohort", x = "year", y = "number of publications") + scale_y_continuous(limits = c(0,
    10)) + scale_x_continuous(limits = c(1990, 2024)) + theme_minimal() + theme(plot.margin = margin(t = 1,
    r = 1, b = 1, l = 1, unit = "cm"))

# theme(plot.background = element_rect(colour='black', fill=NA, size=1))

Figure1_m <- ggplot(data_npubs_aggyear_m, aes(x = pub_year, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + labs(color = "Phd Cohort", x = "year", y = "number of publications") + scale_y_continuous(limits = c(0,
    10)) + scale_x_continuous(limits = c(1990, 2024)) + theme_minimal() + theme(legend.position = "none") +
    theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "cm"))


Figure1_w <- ggplot(data_npubs_aggyear_w, aes(x = pub_year, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + labs(color = "Phd Cohort", x = "year", y = "number of publications") + scale_y_continuous(limits = c(0,
    10)) + scale_x_continuous(limits = c(1990, 2024)) + theme_minimal() + theme(legend.position = "none") +
    theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "cm"))

Figure1_revision <- ggarrange(Figure1, ggarrange(Figure1_m, Figure1_w, ncol = 2, labels = c("men", "women"),
    hjust = c(-2.5, -1.5)), nrow = 2, labels = "total", hjust = -2.5, common.legend = TRUE, legend = "right")

Figure1_revision

# ggexport(Figure1_revision, filename = './data/processed/20250531Figure1_rev.pdf') { png(file =
# './data/processed/20250531Figure1_rev.png', width = 1600, height = 1600, units = 'px', res = 200)
# Figure1_revision } dev.off() # # ggexport(Figure1_revision, filename =
# './data/processed/20250531figure1_rev.png', width = 1600, height=800, res=300, pointsize = 8)

df_ppf %>%
    group_by(time, phd_year) %>%
    summarize(fr = sum(npubs, na.rm = TRUE), naut = n(), nauthors = length(unique(id)), covid = first(covid),
        ) -> data_npubs_aggyear

data_covid <- data_npubs_aggyear[data_npubs_aggyear$covid == 1, ]

# ggplot(data_npubs_aggyear, aes(x=time, y=fr, group=phd_year, color=as.factor(phd_year))) +
# geom_line()

Figure2 <- ggplot(data_npubs_aggyear, aes(x = time, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + geom_point(data = data_covid, size = 1, show.legend = FALSE) + labs(color = "Phd Cohort",
    x = "career age", y = "number of publications", caption = "Note: Points refer to the COVID-19 years") +
    scale_y_continuous(limits = c(0, 10)) + theme_minimal() + theme(plot.background = element_rect(colour = "black",
    fill = NA, size = 1), plot.caption = element_text(hjust = 0))

# fsave(Figure2, 'Figure2')
Figure2

df_ppf %>%
    filter(gender2 == "men") %>%
    group_by(time, phd_year) %>%
    summarize(fr = sum(npubs, na.rm = TRUE), naut = n(), nauthors = length(unique(id)), covid = first(covid),
        ) -> data_npubs_aggyear_m

data_covid_m <- data_npubs_aggyear_m[data_npubs_aggyear_m$covid == 1, ]

# ggplot(data_npubs_aggyear, aes(x=time, y=fr, group=phd_year, color=as.factor(phd_year))) +
# geom_line()

Figure2_m <- ggplot(data_npubs_aggyear_m, aes(x = time, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + geom_point(data = data_covid_m, size = 1, show.legend = FALSE) + labs(color = "Phd Cohort",
    x = "career age", y = "number of publications", caption = "Note: Points refer to the COVID-19 years") +
    scale_y_continuous(limits = c(0, 10)) + theme_minimal() + theme(plot.background = element_rect(colour = "black",
    fill = NA, size = 1), plot.caption = element_text(hjust = 0))

# fsave(Figure2, 'Figure2')
Figure2_m

df_ppf %>%
    filter(gender2 == "women") %>%
    group_by(time, phd_year) %>%
    summarize(fr = sum(npubs, na.rm = TRUE), naut = n(), nauthors = length(unique(id)), covid = first(covid),
        ) -> data_npubs_aggyear_w

data_covid_w <- data_npubs_aggyear_w[data_npubs_aggyear_w$covid == 1, ]

# ggplot(data_npubs_aggyear, aes(x=time, y=fr, group=phd_year, color=as.factor(phd_year))) +
# geom_line()

Figure2_w <- ggplot(data_npubs_aggyear_w, aes(x = time, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + geom_point(data = data_covid_w, size = 1, show.legend = FALSE) + labs(color = "Phd Cohort",
    x = "career age", y = "number of publications", caption = "Note: Points refer to the COVID-19 years") +
    scale_y_continuous(limits = c(0, 10)) + theme_minimal() + theme(plot.background = element_rect(colour = "black",
    fill = NA, size = 1), plot.caption = element_text(hjust = 0))

# fsave(Figure2, 'Figure2')
Figure2_w

Figure2 <- ggplot(data_npubs_aggyear, aes(x = time, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + geom_point(data = data_covid, size = 1, show.legend = FALSE) + labs(color = "Phd Cohort",
    x = "career age", y = "number of publications") + scale_y_continuous(limits = c(0, 10)) + theme_minimal()
theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "cm"))

Figure2_m <- ggplot(data_npubs_aggyear_m, aes(x = time, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + geom_point(data = data_covid_m, size = 1, show.legend = FALSE) + labs(color = "Phd Cohort",
    x = "career age", y = "number of publications") + scale_y_continuous(limits = c(0, 10)) + theme_minimal() +
    theme(legend.position = "none") + theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "cm"))

Figure2_w <- ggplot(data_npubs_aggyear_w, aes(x = time, y = fr/nauthors, group = phd_year, color = as.factor(phd_year))) +
    geom_line() + geom_point(data = data_covid_w, size = 1, show.legend = FALSE) + labs(color = "Phd Cohort",
    x = "career age", y = "number of publications") + scale_y_continuous(limits = c(0, 10)) + theme_minimal() +
    theme(legend.position = "none") + theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "cm"))


Figure2_revision <- ggarrange(Figure2, ggarrange(Figure2_m, Figure2_w, ncol = 2, labels = c("men", "women"),
    hjust = c(-2.5, -1.5)), nrow = 2, labels = "total", hjust = -2.5, common.legend = TRUE, legend = "right")



Figure2_revision <- annotate_figure(Figure2_revision, bottom = text_grob("Note: Points refer to the COVID-19 years",
    color = "black", just = "left", x = 0.05, y = 1.4, size = 8, lineheight = 0.6))


Figure2_revision

# fsave(Figure2_revision, file='Figure2_revision.rda') ggexport(Figure2_revision, filename =
# './data/processed/20250531Figure2_rev.pdf') { png(file =
# './data/processed/20250531Figure2_rev.png', width = 1600, height = 1600, units = 'px', res = 200)
# Figure2_revision } dev.off()



# ggexport(Figure1_revision, filename = './data/processed/figure1.png', width = 1600, height=800,
# res=300, pointsize = 8)

df_ppf_sel <- df_ppf %>%
    filter(gender2 == "men" | gender2 == "women")
describe(df_ppf_sel$npubs)

describeBy(df_ppf_sel$npubs, group = df_ppf_sel$gender2)

df_ppf_sel[df_ppf_sel$npubs == max(df_ppf_sel$npubs), ]$id

tdf_pubs <- df_pubs[df_pubs$id == df_ppf_sel[df_ppf_sel$npubs == max(df_ppf_sel$npubs), ]$id, ]

fshowdf(tdf_pubs)