Analyses_robustness_ols_log.Rmd

---
title: "Robustness checks: linear models with log-transformed variables"
author: "Ana Macanovic"
date: "2024-07-24"
---

Various linear model specifications, but including log-transformed variables of interest.
All specifications include standard errors clustered at the professor level.


Load the packages:
```{r message=  F, warning = F, eval = T}
source("helper_functions.R")

packages_to_load <- c("readr", "dplyr", "tidyr", 
                      "ggplot2", "cowplot",
                      "tidyverse", "RPostgres", 
                      "lubridate", "lmtest", 
                      "sandwich", "ggpubr", 
                      "knitr", "scales", 
                      "ggeffects", "flextable", 
                      "officer", "DescTools",
                      "gglorenz", "corrplot")

fpackage_check(packages_to_load)


# For full reproducibility, load the packages with groundhog using the code below instead
# of the fpackage_check function

# library(groundhog)
# groundhog.library(packages_to_load, date = "2024-4-23")
```


```{r include=FALSE}
opts_chunk$set(echo = TRUE)
opts_chunk$set(eval = TRUE)
opts_chunk$set(warning = FALSE)
opts_chunk$set(message = FALSE)
```


Load the panel dataset:
```{r warning = F, message = F}
prof_panel_filter <- read_csv("panel_datasets/prof_panel_final_26_7.csv")
```


Make a list to save all the results into:
```{r}
all_reg_rob_list <- list()
```


## Main models - log-transformed variables

### Printed news attention
```{r warning = F}
news_formula_log_model <- "news_all_log ~ inferred_gender + news_all_l_log + cited_by_total_all_l_log + alt_online_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

news_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                       lm_formula_list = news_formula_log_model,
                                       year_cutoff_upper = 2023,
                                       year_cutoff_lower = 2012)
```

### Online news attention 

```{r}
online_news_formula_log_model <- "alt_online_all_log ~ inferred_gender + alt_online_all_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

online_news_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                              lm_formula = online_news_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```

### Twitter attention models

```{r}
twitter_formula_log_model <- "alt_twitter_log ~ inferred_gender + alt_twitter_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_online_all_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

twitter_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                              lm_formula_list = twitter_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```


### Combined coefficients

Print out selected coefficients in a table:
```{r}
table_models_log <- neat_regression_table_scopus(news_log_model[[1]],
                                          online_news_log_model[[1]],
                                          twitter_log_model[[1]])

# add to the list
all_reg_rob_list[['linear_log']] <- table_models_log

(table_models_save <- table_models_log %>%
    regulartable() %>% 
    set_caption("Main model - log transformed")%>%
  autofit())


word_document_name <-
    read_docx() %>%
    body_add_flextable(table_models_save) %>%
    print(target = "results/supplement_tables/main_model_table_log.docx")
```


## Main models - log-transformed variables, but last year instead of totals

### Printed news attention
```{r warning = F}
news_formula_log_model <- "news_all_log ~ inferred_gender + news_all_l_log + cited_by_l_log + alt_online_all_l_log + alt_twitter_l_log +coa_tot_cited_by_l_log + coa_tot_online_all_l_log + coa_tot_twitter_l_log+years_since_first_pub + as.factor(year)"

news_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                       lm_formula_list = news_formula_log_model,
                                       year_cutoff_upper = 2023,
                                       year_cutoff_lower = 2012)
```

### Online news attention 

```{r}
online_news_formula_log_model <- "alt_online_all_log ~ inferred_gender + alt_online_all_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

online_news_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                              lm_formula = online_news_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```

### Twitter attention models

```{r}
twitter_formula_log_model <- "alt_twitter_log ~ inferred_gender + alt_twitter_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_online_all_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

twitter_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                              lm_formula_list = twitter_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```


### Combined coefficients

Print out selected coefficients in a table:
```{r}
table_models_log <- neat_regression_table_scopus(news_log_model[[1]],
                                          online_news_log_model[[1]],
                                          twitter_log_model[[1]])

# add to the list
all_reg_rob_list[['linear_log']] <- table_models_log

(table_models_save <- table_models_log %>%
    regulartable() %>% 
    set_caption("Main model - log transformed")%>%
  autofit())


word_document_name <-
    read_docx() %>%
    body_add_flextable(table_models_save) %>%
    print(target = "results/supplement_tables/main_model_table_log_last_year.docx")
```


## Main models - log-transformed variables and no first lag

### Printed news attention
```{r warning = F}
news_formula_log_model <- "news_all_log ~ inferred_gender + cited_by_total_all_l_log + alt_online_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

news_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                       lm_formula_list = news_formula_log_model,
                                       year_cutoff_upper = 2023,
                                       year_cutoff_lower = 2012)
```

### Online news attention 

```{r}
online_news_formula_log_model <- "alt_online_all_log ~ inferred_gender  + cited_by_total_all_l_log + news_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

online_news_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                              lm_formula = online_news_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```

### Twitter attention models

```{r}
twitter_formula_log_model <- "alt_twitter_log ~ inferred_gender  + cited_by_total_all_l_log + news_all_total_l_log + alt_online_all_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

twitter_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter,
                                              lm_formula_list = twitter_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```


### Combined coefficients

Print out selected coefficients in a table:
```{r}
table_models_log <- neat_regression_table_scopus(news_log_model[[1]],
                                          online_news_log_model[[1]],
                                          twitter_log_model[[1]])

# add to the list
all_reg_rob_list[['linear_log_nolag']] <- table_models_log

(table_models_save <- table_models_log %>%
    regulartable() %>% 
    set_caption("Main model - log transformed")%>%
  autofit())


word_document_name <-
    read_docx() %>%
    body_add_flextable(table_models_save) %>%
    print(target = "results/supplement_tables/main_model_table_log_nolag.docx")
```



## Main models - log-transformed variables excluding 2020 and 2021

### Printed news attention
```{r warning = F}
news_formula_log_model <- "news_all_log ~ inferred_gender + news_all_l_log + cited_by_total_all_l_log + alt_online_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

news_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_2020,
                                       lm_formula_list = news_formula_log_model,
                                       year_cutoff_upper = 2023,
                                       year_cutoff_lower = 2012)
```

### Online news attention 

```{r}
online_news_formula_log_model <- "alt_online_all_log ~ inferred_gender + alt_online_all_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

online_news_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_2020,
                                              lm_formula = online_news_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```

### Twitter attention models

```{r}
twitter_formula_log_model <- "alt_twitter_log ~ inferred_gender + alt_twitter_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_online_all_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

twitter_log_model <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_2020,
                                              lm_formula_list = twitter_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```


### Combined coefficients

Print out selected coefficients in a table:
```{r}
table_models_log <- neat_regression_table_scopus(news_log_model[[1]],
                                          online_news_log_model[[1]],
                                          twitter_log_model[[1]])

# add to the list
all_reg_rob_list[['linear_log_2019']] <- table_models_log

(table_models_save <- table_models_log %>%
    regulartable() %>% 
    set_caption("Main model - log transformed until 2020")%>%
  autofit())


word_document_name <-
    read_docx() %>%
    body_add_flextable(table_models_save) %>%
    print(target = "results/supplement_tables/main_model_table_log_2019.docx")
```


## Main models - log-transformed excluding those who held a public post

Load data on profesors who had a public position:
```{r}
# Connect to the database:
# fill in own credentials
port <- 5432
user <- "postgres"
password <- "dutchmediaprofssql"
database_name <- "postgres"


con <- dbConnect(Postgres(),
                 dbname= database_name,
                 port = port,
                 user = user, 
                 password = password)

con # Checks connection is working

public_roles <- dbReadTable(con, "govt_positions")

prof_panel_filter_public <- filter(prof_panel_filter,
                                   ! profile_id %in% public_roles$profile_id)
```


### Printed news attention

```{r warning = F}
news_formula_log_model <- "news_all_log ~ inferred_gender + news_all_l_log + cited_by_total_all_l_log + alt_online_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

news_model_no_public <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_public,
                                       lm_formula = news_formula_log_model,
                                       year_cutoff_upper = 2023,
                                       year_cutoff_lower = 2012)
```

### Online news attention 
```{r}
online_news_formula_log_model <- "alt_online_all_log ~ inferred_gender + alt_online_all_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

online_news_model_no_public <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_public,
                                              lm_formula = online_news_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```

### Twitter/X
```{r}
twitter_formula_log_model <- "alt_twitter_log ~ inferred_gender + alt_twitter_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_online_all_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

twitter_model_no_public <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_public,
                                              lm_formula = twitter_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```


### Combined coefficients

Combine this all to plot the coefficients in one plot.
```{r}
options(scipen=999)

table_models_no_public_log <- neat_regression_table_scopus(news_model_no_public[[1]],
                                                    online_news_model_no_public[[1]],
                                                    twitter_model_no_public[[1]])

# add to the list
all_reg_rob_list[['linear_no_public_log']] <- table_models_no_public_log

(table_models_save <- table_models_no_public_log %>%
    regulartable() %>% 
    set_caption("Main model - log transformed excluding professors with positions in government")%>%
  autofit())


word_document_name <-
    read_docx() %>%
    body_add_flextable(table_models_save) %>%
    print(target = "results/supplement_tables/main_model_table_no_public_log.docx")
```

## Main models - excluding those who were an emeritus professor

Load data on profesors who had a public position:
```{r}
emeritus_roles <- dbReadTable(con, "emeritus_info")

emeritus_roles <- filter(emeritus_roles, emeritus == TRUE)

prof_panel_filter_emeritus <- filter(prof_panel_filter,
                                   ! profile_id %in% emeritus_roles$profile_id)
```


### Printed news attention

```{r warning = F}
news_formula_log_model <- "news_all_log ~ inferred_gender + news_all_l_log + cited_by_total_all_l_log + alt_online_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

news_model_no_emer <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_emeritus,
                                       lm_formula = news_formula_log_model,
                                       year_cutoff_upper = 2023,
                                       year_cutoff_lower = 2012)
```

### Online news attention 
```{r}
online_news_formula_log_model <- "alt_online_all_log ~ inferred_gender + alt_online_all_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_twitter_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

online_news_model_no_emer <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_emeritus,
                                              lm_formula = online_news_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```

### Twitter/X
```{r}
twitter_formula_log_model <- "alt_twitter_log ~ inferred_gender + alt_twitter_l_log + cited_by_total_all_l_log + news_all_total_l_log + alt_online_all_total_l_log +coa_tot_cited_by_total_l_log + coa_tot_online_all_total_l_log + coa_tot_twitter_total_l_log+years_since_first_pub + as.factor(year)"

twitter_model_no_emer <- lm_fitter_cl_robust_scopus(panel_dataset = prof_panel_filter_emeritus,
                                              lm_formula = twitter_formula_log_model,
                                              year_cutoff_upper = 2023,
                                              year_cutoff_lower = 2012)
```


### Combined coefficients

Combine this all to plot the coefficients in one plot.
```{r}
options(scipen=999)

table_models_no_emer <- neat_regression_table_scopus(news_model_no_emer[[1]],
                                           online_news_model_no_emer[[1]],
                                           twitter_model_no_emer[[1]])

# add to the list
all_reg_rob_list[['linear_no_emer']] <- table_models_no_emer

(table_models_save <- table_models_no_emer %>%
    regulartable() %>% 
    set_caption("Main model - excluding professors with emeritus positions")%>%
  autofit())


word_document_name <-
    read_docx() %>%
    body_add_flextable(table_models_save) %>%
    print(target = "results/supplement_tables/main_model_table_no_emeritus_log.docx")
```

Save the list:
```{r}
saveRDS(all_reg_rob_list, "results/ols_log_models.RDS")
```