##########################################################
##########################################################
# Preprocessing of print media corpus
### SNAP MEDIA ANALYSIS
### LAST UPDATED: 29-Oct-2019
### CONTACT: BENJAMIN.CHRISINGER@SPI.OX.AC.UK 
##########################################################
##########################################################

#load required packages
library(tm)
library(wordcloud)
library(RColorBrewer)
library(plyr)
library(ggplot2)
library(RTextTools)
library(e1071)
library(lda)
library(tm)
library(stm)
library(stringr)
library(urltools)
library(parallel)

#read in chunks of print corpus and assemble into one large text dataset
ln0A = read_file("LNmedia_SearchTerms_PRE2009")
s0A= strsplit(ln0A, split="DOCUMENTS\\\n\\\n\\\n")
df0A = as.data.frame(s0A)
colnames(df0A) = "text"

ln0B = read_file("LNmedia_SearchTerms_PRE2010")
s0B= strsplit(ln0B, split="DOCUMENTS\\\n\\\n\\\n")
df0B = as.data.frame(s0B)
colnames(df0B) = "text"

ln0C = read_file("LNmedia_SearchTerms")
s0C= strsplit(ln0C, split="DOCUMENTS\\\n\\\n\\\n")
df0C = as.data.frame(s0C)
colnames(df0C) = "text"

ln1 = read_file("LexisNexis_SNAPSearchTerms1")
s1= strsplit(ln1, split="DOCUMENTS\\\n\\\n\\\n")
df1 = as.data.frame(s1)
colnames(df1) = "text"

ln2A = read_file("LexisNexis_SNAPSearchTerms2A")
s2A= strsplit(ln2A, split="DOCUMENTS\\\n\\\n\\\n\\\n")
df2A = as.data.frame(s2A)
colnames(df2A) = "text"

ln2B = read_file("LexisNexis_SNAPSearchTerms2A")
s2B= strsplit(ln2B, split="DOCUMENTS\\\n\\\n\\\n\\\n")
df2B = as.data.frame(s2B)
colnames(df2B) = "text"

ln3 = read_file("LexisNexis_SNAPSearchTerms3")
s3= strsplit(ln3, split="DOCUMENTS\\\n\\\n\\\n")
df3 = as.data.frame(s3)
colnames(df3) = "text"

ln4A = read_file("LexisNexis_SNAPSearchTerms4A")
s4A= strsplit(ln4A, split="DOCUMENTS\\\n\\\n\\\n\\\n")
df4A = as.data.frame(s4A)
colnames(df4A) = "text"

ln4B = read_file("LexisNexis_SNAPSearchTerms4B")
s4B= strsplit(ln4B, split="DOCUMENTS\\\n\\\n\\\n")
df4B = as.data.frame(s4B)
colnames(df4B) = "text"

ln4C = read_file("LexisNexis_SNAPSearchTerms4C")
s4C= strsplit(ln4C, split="DOCUMENTS\\\n\\\n\\\n")
df4C = as.data.frame(s4C)
colnames(df4C) = "text"

ln4D = read_file("LexisNexis_SNAPSearchTerms4D")
s4D= strsplit(ln4D, split="DOCUMENTS\\\n\\\n\\\n")
df4D = as.data.frame(s4D)
colnames(df4D) = "text"

ln4E = read_file("LexisNexis_SNAPSearchTerms4E")
s4E= strsplit(ln4E, split="DOCUMENTS\\\n\\\n\\\n")
df4E = as.data.frame(s4E)
colnames(df4E) = "text"

ln4F = read_file("LexisNexis_SNAPSearchTerms4F")
s4F= strsplit(ln4F, split="DOCUMENTS\\\n\\\n\\\n")
df4F = as.data.frame(s4F)
colnames(df4F) = "text"

ln4G = read_file("LexisNexis_SNAPSearchTerms4G")
s4G= strsplit(ln4G, split="DOCUMENTS\\\n\\\n\\\n")
df4G = as.data.frame(s4G)
colnames(df4G) = "text"

ln4H = read_file("LexisNexis_SNAPSearchTerms4H")
s4H= strsplit(ln4H, split="DOCUMENTS\\\n\\\n\\\n")
df4H = as.data.frame(s4H)
colnames(df4H) = "text"

ln4I = read_file("LexisNexis_SNAPSearchTerms4I")
s4I= strsplit(ln4I, split="DOCUMENTS\\\n\\\n\\\n")
df4I = as.data.frame(s4I)
colnames(df4I) = "text"

ln5 = read_file("LexisNexis_SNAPSearchTerms5")
s5= strsplit(ln5, split="DOCUMENTS\\\n\\\n\\\n")
df5 = as.data.frame(s5)
colnames(df5) = "text"

df = rbind(df0A,df0B,df0C,df1,df2A,df2B,df3,df4A,df4B,df4C,df4D,df4E,df4F,df4G,df4H,df4I,df5)
#rm(df1,df2A,df2B,df3,df4A,df4B,df4C,df4D,df4E,df4F,df4G,df4H,df4I,df5,s1,s2A,s2B,s3,s4A,s4B,s4C,s4D,s4E,s4F,s4G,s4H,s4I,s5,ln1,ln2A,ln2B,ln3,ln4A,ln4B,ln4C,ln4D,ln4E,ln4F,ln4G,ln4H,ln4I,ln5)

# start formatting the text as individual documents for analysis

df$text = as.character(df$text)
df$ID = 1:nrow(df)
df$text = gsub("\r?\n|\r", " ", df$text)

# remove any urls from lexis
df$text = ifelse(str_detect(df$text,"http://www-lexisnexis-com.stanford.")==TRUE,
                 str_sub(df$text, str_locate(df$text,"inline=y")[,2]+2),
                 df$text)

# extract date
df$date = str_sub(df$text,str_locate(df$text,paste(format(ISOdate(2000, 1:12, 1), "%B"),collapse = "|"))[,1],
                  str_locate(df$text,"[1-2][0-9][0-9][0-9]")[,2])
df$year = as.integer(str_sub(df$date,-4,-1))
df$month_abb = str_sub(df$date,
                       str_locate(df$date,paste(format(ISOdate(2000, 1:12, 1), "%b"),collapse = "|"))[,1],
                       str_locate(df$date,paste(format(ISOdate(2000, 1:12, 1), "%b"),collapse = "|"))[,2])
df$month = match(df$month_abb,month.abb)
df$nmonth = (2018-df$year)*12 + df$month

# extract location and source
df$source0 = str_sub(df$text,1,str_locate(df$text,"[1-2][0-9][0-9][0-9]")[,2])
df$source1 = str_sub(df$text,1,str_locate(df$text,paste(format(ISOdate(2000, 1:12, 1), "%B"),collapse = "|"))[,1]-1)
df$source1 = trimws(df$source1)

df$state = str_sub(df$source1,str_locate(df$source1,"\\(")[,1]+1,str_locate(df$source1,"\\)")[,2]-1)

df$source = ifelse(str_detect(df$source1,"\\(|\\s\\s\\s")==TRUE,
                   str_sub(df$source1,1,sapply(c(str_locate(df$source1,"\\(")[,1]-1,str_locate(df$source1,"\\s\\s\\s")[,1]),min)),
                   df$source1)
df$source = trimws(df$source)

# extract publication type and section
df$type0 = str_sub(df$text,str_locate(df$text,"TYPE: ")[,1]+6,str_locate(df$text,"TYPE: ")[,2]+30)
df$type = str_sub(df$type0,1,str_locate(df$type0,"\\s\\s\\s|Copy")[,1]-1)

df$section = str_sub(df$text, str_locate(df$text,"SECTION: ")[,2]+1,str_locate(df$text,";")[,2]-1)

# extract document text
df$document = str_sub(df$text,(regexpr("LENGTH:", df$text)+18),
                      (regexpr("LOAD-DATE:", df$text)))

dat = data.frame("ID" = as.factor(df$ID),
                 "source" = df$source,
                 "type" = df$type,
                 "section" = df$section,
                 "state" = df$state,
                 "date" = df$date,
                 "year" = df$year,
                 "month" = df$month,
                 "nmonth" = df$nmonth,
                 "documents" = as.character(df$document),
                 stringsAsFactors = F) #n=79,418

# remove duplicates (n=64,203)
dat = unique(dat[,-1])

dat = dat[!is.na(dat$nmonth),] #n=63622
dat = dat[!is.na(dat$source),] #n=62381
dat = dat[!is.na(dat$documents),] #n=62381

### START REDUCING DATABASE
# Make sure to include data that specifically reference the program (n=62,381)
dat$check0 = ifelse(
  (regexpr("food stamp", dat$documents, ignore.case = T) != -1 |
     regexpr("Supplemental Nutrition", dat$documents, ignore.case = T) != -1 |
     regexpr("Supplement Nutrition Assistance Program", dat$documents, ignore.case = T) != -1 |
     regexpr("Supplemental Assistance Nutrition Program", dat$documents, ignore.case = T) != -1 |
     regexpr("SNAP program", dat$documents, ignore.case = F) != -1 |
     regexpr("SNAP Program", dat$documents, ignore.case = F) != -1 |
     regexpr("CalFresh", dat$documents, ignore.case = T) != -1 |
     regexpr("SNAP benefits", dat$documents, ignore.case = F) != -1 |
     regexpr("SNAP-Ed", dat$documents, ignore.case = F) != -1 |
     regexpr("3SquaresVT", dat$documents, ignore.case = T) != -1) == "TRUE", 1, 0)

# Identify data that use the word "SNAP" that are actually about SNAP
dat$check1 = ifelse(
  ((regexpr("SNAP",dat$documents, ignore.case = F) != -1 | regexpr("S.N.A.P.",dat$documents, ignore.case = F) != -1) &
     (regexpr("food assistance",dat$documents, ignore.case = T) != -1 |
        (regexpr("nutrition assistance",dat$documents, ignore.case = T) != -1 |
           (regexpr("USDA",dat$documents, ignore.case = F) != -1 |
              (regexpr("EBT",dat$documents, ignore.case = F) != -1 |
                 (regexpr("E.B.T.",dat$documents, ignore.case = F) != -1 |
                    (regexpr("U.S.D.A.",dat$documents, ignore.case = F) != -1 |
                       dat$check0 == 1))))))), 1, 0)

dat$check2 = ifelse((
  dat$check0 == 1 | 
    dat$check1 == 1) == "TRUE",1,0) 

# If an ambiguous state-specific program name also occurs with the state's name, include it. (n=72,713)
dat$check3 = ifelse(regexpr("Basic Food Program", dat$documents, ignore.case = T) != -1 &
                      (regexpr("Washington", dat$documents, ignore.case = T) != -1 |
                         regexpr("WA", dat$documents, ignore.case = F) != -1 |
                         dat$check2 == 1), 1,
                    ifelse(regexpr("Food and Nutrition Services", dat$documents, ignore.case = T) != -1 &
                             (regexpr("North Carolina", dat$documents, ignore.case = T) != -1 |
                                regexpr("NC", dat$documents, ignore.case = F) != -1 |
                                dat$check2 == 1), 1,
                           ifelse(regexpr("Food Assistance Program", dat$documents, ignore.case = T) != -1 &
                                    (regexpr("Alabama", dat$documents, ignore.case = T) != -1 |
                                       regexpr("AL", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Colorado", dat$documents, ignore.case = T) != -1 |
                                       regexpr("CO", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Florida", dat$documents, ignore.case = T) != -1 |
                                       regexpr("FL", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Iowa", dat$documents, ignore.case = T) != -1 |
                                       regexpr("IA", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Kansas", dat$documents, ignore.case = T) != -1 |
                                       regexpr("KS", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Ohio", dat$documents, ignore.case = T) != -1 |
                                       regexpr("OH", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1), 1,
                                  ifelse(regexpr("Food Supplement Program", dat$documents, ignore.case = T) != -1 &
                                           (regexpr("Delaware", dat$documents, ignore.case = T) != -1 |
                                              regexpr("DE", dat$documents, ignore.case = F) != -1 |
                                              dat$check2 == 1) |
                                           (regexpr("Maine", dat$documents, ignore.case = T) != -1 |
                                              regexpr("ME", dat$documents, ignore.case = F) != -1 |
                                              dat$check2 == 1) |
                                           (regexpr("Maryland", dat$documents, ignore.case = T) != -1 |
                                              regexpr("MD", dat$documents, ignore.case = F) != -1 |
                                              dat$check2 == 1) |
                                           (regexpr("Michigan", dat$documents, ignore.case = T) != -1 |
                                              regexpr("MI", dat$documents, ignore.case = F) != -1 |
                                              dat$check2 == 1), 1,
                                         ifelse(regexpr("FoodShare", dat$documents, ignore.case = T) != -1 &
                                                  (regexpr("Wisconsin", dat$documents, ignore.case = T) != -1 |
                                                     regexpr("WI", dat$documents, ignore.case = F) != -1 |
                                                     dat$check2 == 1),1,
                                                ifelse(regexpr("Nutrition Assistance", dat$documents, ignore.case = T) != -1 &
                                                         (regexpr("Arizona", dat$documents, ignore.case = T) != -1 |
                                                            regexpr("AZ", dat$documents, ignore.case = F) != -1 |
                                                            dat$check2 == 1),1,0))))))
dat$check = ifelse(dat$check0 + dat$check1 + dat$check2 + dat$check3 > 0, 1, 0)
table(dat$check0) #n=57,112
table(dat$check1) #n=6,661
table(dat$check2) #n=57,432
dat = data.frame(dat[dat$check == "1",]) #n=58,782

#Save database
save(dat,file=paste0("PrintCorpusFull_Database","_",Sys.Date(),".RData"))

dat2=dat[!is.na(dat$nmonth) & dat$nmonth<350 & dat$nmonth>0,] #58,771

# identify urls from the top news outlets
# https://assets.pewresearch.org/wp-content/uploads/sites/13/2016/06/30143308/state-of-the-news-media-report-2016-final.pdf
newspapers <- read_csv("majorpaper_avgalign.csv")

dat2$source2=trimws(toupper(dat2$source))
dat2$source2= str_replace(dat2$source2, "THE ","")
dat2$source2 = str_replace_all(dat2$source2, "[[:punct:]]", " ")

dat2$top_papers=dat2$source2 %in% newspapers$name

dat2$majorYN = ifelse(dat2$top_papers==TRUE,"major_paper", "other_paper")
dat2$majorYN = as.factor(dat2$majorYN)

#PREPROCESSING - full print corpus dataset
set.seed(1234)
PrintCorpusFull_processed <- textProcessor(dat2$documents, 
                                           metadata = dat2,
                                           verbose = T)
PrintCorpusFull_out <- prepDocuments(PrintCorpusFull_processed$documents,
                                     PrintCorpusFull_processed$vocab,
                                     PrintCorpusFull_processed$meta, 
                                     lower.thresh = nrow(dat2) * 0.01,
                                     upper.thresh = nrow(dat2) * 0.90)
PrintCorpusFull_docs <- PrintCorpusFull_out$documents
PrintCorpusFull_vocab <- PrintCorpusFull_out$vocab
PrintCorpusFull_meta <-PrintCorpusFull_out$meta
save(PrintCorpusFull_processed, file=paste("PrintCorpusFull_processed",Sys.Date(),sep="_"))
save(PrintCorpusFull_out, file=paste("PrintCorpusFull_out",Sys.Date(),sep="_"))

# Assign bias scores to subset of Print Corpus and save
top_bias = dat2[dat2$source2 %in% newspapers$name,]
top_bias = merge(dat2,newspapers,by.x="source2",by.y="name",all=F)
save(top_bias,file=paste0("PrintCorpusBias_Database","_",Sys.Date(),".RData"))

# PREPROCESSING - subset of print corpus 
set.seed(1234)
PrintCorpusBias_processed <- textProcessor(top_bias$documents, 
                                           metadata = top_bias,
                                           verbose = T)
PrintCorpusBias_out <- prepDocuments(PrintCorpusBias_processed$documents,
                                     PrintCorpusBias_processed$vocab,
                                     PrintCorpusBias_processed$meta, 
                                     lower.thresh = nrow(top_bias) * 0.01,
                                     upper.thresh = nrow(top_bias) * 0.90)
PrintCorpusBias_docs <- PrintCorpusBias_out$documents
PrintCorpusBias_vocab <- PrintCorpusBias_out$vocab
PrintCorpusBias_meta <-PrintCorpusBias_out$meta
save(PrintCorpusBias_processed, file=paste("PrintCorpusBias_processed",Sys.Date(),sep="_"))
save(PrintCorpusBias_out, file=paste("PrintCorpusBias_out",Sys.Date(),sep="_"))


##########################################################
##########################################################
# Preprocessing of online media corpus
### SNAP MEDIA ANALYSIS
### LAST UPDATED: 29-Oct-2019
### CONTACT: BENJAMIN.CHRISINGER@SPI.OX.AC.UK 
##########################################################
##########################################################
library(tm)
library(wordcloud)
library(RColorBrewer)
library(plyr)
library(ggplot2)
library(RTextTools)
library(e1071)
library(lda)
library(tm)
library(stm)
library(stringr)
library(urltools)
library(parallel)

#read in web scrape database
articles = read.csv("articles-with-topics.csv")
# n=13,987

#create a variable for the article texts and make sure they'll be recognized as text
articles$documents = as.character(articles$Text)

#remove documents with missing data
articles = remove_missing(articles)
# n=13,940

dat <- data.frame("ID" = as.factor(articles$ID),
                  "year" = as.integer(str_pad(substr(articles$Date,1,4),2,"left","0")),
                  "month" = as.integer(str_pad(substr(articles$Date,5,6),2,"left","0")),
                  "day" = as.integer(str_pad(substr(articles$Date,7,8),2,"left","0")),
                  "Headline" = as.character(articles$Headline),
                  "Url" = as.character(trimws(articles$Url)),
                  "documents" = as.character(articles$Text))

# remove duplicates 
dat = unique(dat[,-1])

dat$Date = paste(dat$year,dat$month,dat$day,sep="-")
dat$Date = as.Date(dat$Date, "%Y-%m-%d")
dat$numdate = as.numeric(dat$Date - min(dat$Date,na.rm = T))
dat$numweek = as.integer(dat$numdate / 7)
dat$nummonth = as.integer(dat$numdate / 30)

#remove documents outside the search timeline
dat = dat[dat$year > 1989 & dat$year < 2018,]

dat$domain = domain(dat$Url)
dat$domain0 = suffix_extract(dat$domain)$domain

### START REDUCING DATABASE

# Make sure to include dat that specifically reference the program (n=13,126)
dat$check0 = ifelse(
  (regexpr("food stamp", dat$documents, ignore.case = T) != -1 |
     regexpr("Supplemental Nutrition", dat$documents, ignore.case = T) != -1 |
     regexpr("Supplement Nutrition Assistance Program", dat$documents, ignore.case = T) != -1 |
     regexpr("Supplemental Assistance Nutrition Program", dat$documents, ignore.case = T) != -1 |
     regexpr("SNAP program", dat$documents, ignore.case = F) != -1 | 
     regexpr("SNAP Program", dat$documents, ignore.case = F) != -1 | 
     regexpr("CalFresh", dat$documents, ignore.case = T) != -1 |
     regexpr("SNAP benefits", dat$documents, ignore.case = F) != -1 |
     regexpr("SNAP-Ed", dat$documents, ignore.case = F) != -1 |
     regexpr("3SquaresVT", dat$documents, ignore.case = T) != -1) == "TRUE", 1, 0)

# Identify dat that use the word "SNAP" that are actually about SNAP.
dat$check1 = ifelse(
  ((regexpr("SNAP",dat$documents, ignore.case = F) != -1 | regexpr("S.N.A.P.",dat$documents, ignore.case = F) != -1) &
     (regexpr("food assistance",dat$documents, ignore.case = T) != -1 |
        (regexpr("nutrition assistance",dat$documents, ignore.case = T) != -1 |
           (regexpr("USDA",dat$documents, ignore.case = F) != -1 |
              (regexpr("EBT",dat$documents, ignore.case = F) != -1 |
                 (regexpr("U.S.D.A.",dat$documents, ignore.case = F) != -1 |
                    dat$check0 == 1)))))), 1, 0)

dat$check2 = ifelse((
  dat$check0 == 1 | 
    dat$check1 == 1) == "TRUE",1,0)

# If an ambiguous state-specific program name also occurs with the state's name, include it.
dat$check3 = ifelse(regexpr("Basic Food Program", dat$documents, ignore.case = T) != -1 &
                      (regexpr("Washington", dat$documents, ignore.case = T) != -1 |
                         regexpr("WA", dat$documents, ignore.case = F) != -1 |
                         dat$check2 == 1), 1,
                    ifelse(regexpr("Food and Nutrition Services", dat$documents, ignore.case = T) != -1 &
                             (regexpr("North Carolina", dat$documents, ignore.case = T) != -1 |
                                regexpr("NC", dat$documents, ignore.case = F) != -1 |
                                dat$check2 == 1), 1,
                           ifelse(regexpr("Food Assistance Program", dat$documents, ignore.case = T) != -1 &
                                    (regexpr("Alabama", dat$documents, ignore.case = T) != -1 |
                                       regexpr("AL", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Colorado", dat$documents, ignore.case = T) != -1 |
                                       regexpr("CO", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Florida", dat$documents, ignore.case = T) != -1 |
                                       regexpr("FL", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Iowa", dat$documents, ignore.case = T) != -1 |
                                       regexpr("IA", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Kansas", dat$documents, ignore.case = T) != -1 |
                                       regexpr("KS", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1) |
                                    (regexpr("Ohio", dat$documents, ignore.case = T) != -1 |
                                       regexpr("OH", dat$documents, ignore.case = F) != -1 |
                                       dat$check2 == 1), 1,
                                  ifelse(regexpr("Food Supplement Program", dat$documents, ignore.case = T) != -1 &
                                           (regexpr("Delaware", dat$documents, ignore.case = T) != -1 |
                                              regexpr("DE", dat$documents, ignore.case = F) != -1 |
                                              dat$check2 == 1) |
                                           (regexpr("Maine", dat$documents, ignore.case = T) != -1 |
                                              regexpr("ME", dat$documents, ignore.case = F) != -1 |
                                              dat$check2 == 1) |
                                           (regexpr("Maryland", dat$documents, ignore.case = T) != -1 |
                                              regexpr("MD", dat$documents, ignore.case = F) != -1 |
                                              dat$check2 == 1) |
                                           (regexpr("Michigan", dat$documents, ignore.case = T) != -1 |
                                              regexpr("MI", dat$documents, ignore.case = F) != -1 |
                                              dat$check2 == 1), 1,
                                         ifelse(regexpr("FoodShare", dat$documents, ignore.case = T) != -1 &
                                                  (regexpr("Wisconsin", dat$documents, ignore.case = T) != -1 |
                                                     regexpr("WI", dat$documents, ignore.case = F) != -1 |
                                                     dat$check2 == 1),1,
                                                ifelse(regexpr("Nutrition Assistance", dat$documents, ignore.case = T) != -1 &
                                                         (regexpr("Arizona", dat$documents, ignore.case = T) != -1 |
                                                            regexpr("AZ", dat$documents, ignore.case = F) != -1 |
                                                            dat$check2 == 1),1,0))))))

#summarize numbers of articles at each check
dat$check = ifelse(dat$check0 + dat$check1 + dat$check2 + dat$check3 > 0, 1, 0)
table(dat$check0) #n=10,110
table(dat$check1) #n=5,508
table(dat$check2) #n=10,171
table(dat$check) #n=11,105
dat = data.frame(dat[dat$check == "1",])
# n = 11,105

# identify urls from the top news outlets
# https://assets.pewresearch.org/wp-content/uploads/sites/13/2016/06/30143308/state-of-the-news-media-report-2016-final.pdf
top_print = read.csv("top_print.csv")
top_print$domain = as.character(trimws(top_print$url))
top_print_domain = suffix_extract(top_print$domain)
dat_domain = suffix_extract(dat$domain)
top_print_domain$link = tolower(ifelse(is.na(top_print_domain$subdomain),paste(top_print_domain$domain,top_print_domain$suffix,sep="."),
                                       paste(top_print_domain$subdomain,top_print_domain$domain,top_print_domain$suffix,sep=".")))
dat$link = tolower(ifelse(is.na(dat_domain$subdomain),paste(dat_domain$domain,dat_domain$suffix,sep="."),
                          paste(dat_domain$subdomain,dat_domain$domain,dat_domain$suffix,sep=".")))
dat$top_print_source = as.character(str_match(dat$link,
                                              paste(top_print_domain$link,collapse="|")))
dat$majorYN = ifelse(is.na(dat$top_print_source), 0, 1)
dat$majorYN = as.factor(dat$majorYN)


#Identify top "hard content" url sources, assign partisanship scores
# From: Bakshy, Eytan; Messing, Solomon; Adamic, Lada, 2015, "Replication Data for: Exposure to Ideologically Diverse News and Opinion on Facebook", doi:10.7910/DVN/LDJ7MS, Harvard Dataverse, V2, UNF:6:5LxloFlZg1/Wqd5XntqsgQ==
# https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/LDJ7MS 

hardcontent = read.csv("top500.csv")
hardcontent$domain = as.character(trimws(hardcontent$domain))
hardcontent$domain = domain(hardcontent$domain)
test = suffix_extract(hardcontent$domain)
test2 = suffix_extract(dat$domain)
hardcontent$link = ifelse(is.na(test$subdomain),paste(test$domain,test$suffix,sep="."),
                          paste(test$subdomain,test$domain,test$suffix,sep="."))
dat$link = ifelse(is.na(test2$subdomain),paste(test2$domain,test$suffix,sep="."),
                  paste(test2$subdomain,test2$domain,test2$suffix,sep="."))
dat$hcsource = as.character(str_match(dat$link,
                                      paste(hardcontent$link,collapse="|")))
d1 = merge(dat,hardcontent, by="link", all.x = T)

#save dataset for later. documents without partisanship scores (n=7,255)
dat2 = d1[is.na(d1$avg_align)==TRUE,]
# n = 7,255

#remove documents without partisanship scores
dat = d1[is.na(d1$avg_align)==FALSE,]
# n = 3,850

# FOCUS ON 2010 ONWARDS (n=3799)
OnlineCorpus = dat[dat$Date >= "2010-01-01",]
# n = 3,799

#save both datasets
write.csv(OnlineCorpus,file="OnlineCorpus_bias-score.csv")
write.csv(dat2, file="OnlineCorpus_no-bias-score.csv")
save(OnlineCorpus,file=paste0("OnlineCorpus","_",Sys.Date(),".RData"))

#clean up workspace
rm(d1,df,hardcontent)

#PREPROCESSING
set.seed(1234)
OnlineCorpus_processed <- textProcessor(dat$documents, 
                                        metadata = dat,
                                        verbose = T)

OnlineCorpus_out <- prepDocuments(OnlineCorpus_processed$documents,
                                  OnlineCorpus_processed$vocab,
                                  OnlineCorpus_processed$meta, 
                                  lower.thresh = nrow(dat) * 0.01,
                                  upper.thresh = nrow(dat) * 0.90)
OnlineCorpus_docs <- OnlineCorpus_out$documents
OnlineCorpus_vocab <- OnlineCorpus_out$vocab
OnlineCorpus_meta <-OnlineCorpus_out$meta
save(OnlineCorpus_processed, file=paste("OnlineCorpus_processed",Sys.Date(),sep="_"))
save(OnlineCorpus_out, file=paste("OnlineCorpus_out",Sys.Date(),sep="_"))

write.csv(OnlineCorpus_meta,"OnlineCorpus_meta.csv")





##########################################################
##########################################################
# Run structural topic models and investigate model outputs
### SNAP MEDIA ANALYSIS
### LAST UPDATED: 29-Oct-2019
### CONTACT: BENJAMIN.CHRISINGER@SPI.OX.AC.UK 
##########################################################
##########################################################
#load required packages
library(readr)
library(tm)
library(stm)
library(stringr)
library(gsubfn)
library(qdapRegex)
library(parallel)
install.packages("quartzFonts",dependencies = T)
quartzFonts(avenir = c("Avenir Book", "Avenir Black", "Avenir Book Oblique", 
                       "Avenir Black Oblique"))

####
### RUN STRUCTURAL TOPIC MODELS
####

### 
# ONLINE CORPUS
load("OnlineCorpus_processed_2019-10-29")
load("OnlineCorpus_out_2019-10-29")

OnlineCorpus_docs <- OnlineCorpus_out$documents
OnlineCorpus_vocab <- OnlineCorpus_out$vocab
OnlineCorpus_meta <-OnlineCorpus_out$meta
#run selectModel for different numbers of topics (3-10)
#covariates: time (by week since 01/01/1990), whether or not a major print newspaper, and political alignment (from Top 500)
OnlineCorpus_snapStorage <- manyTopics(OnlineCorpus_out$documents, OnlineCorpus_out$vocab, K=c(3:10),
                                       prevalence=~numweek + majorYN + avg_align,
                                       data=OnlineCorpus_meta, seed=1234)
save(OnlineCorpus_snapStorage, file=paste("OnlineCorpus_selectModel_3to10topics_",Sys.Date(),sep=""))
#looks like 8-topic model has best average semcoh/exclusivity

###
## PRINT CORPUS - FULL
load("PrintCorpusFull_out_2019-10-29")

PrintCorpusFull_out_docs <- PrintCorpusFull_out$documents
PrintCorpusFull_out_vocab <- PrintCorpusFull_out$vocab
PrintCorpusFull_out_meta <-PrintCorpusFull_out$meta
#run selectModel for different numbers of topics (3-10)
#covariates: time (by month since 01/01/1990), source
PrintCorpusFull_snapStorage <- manyTopics(PrintCorpusFull_out$documents, PrintCorpusFull_out$vocab, K=c(3:10),
                                          prevalence=~nmonth + majorYN + source2,
                                          data=PrintCorpusFull_out_meta,seed=1234)
save(PrintCorpusFull_snapStorage, file=paste("PrintCorpusFull_selectModel_3to10topics",Sys.Date(),sep="_"))


###
## PRINT CORPUS - WITH BIAS SCORES
load("PrintCorpusBias_out_2019-10-28")

PrintCorpusBias_docs <- PrintCorpusBias_out$documents
PrintCorpusBias_vocab <- PrintCorpusBias_out$vocab
PrintCorpusBias_meta <-PrintCorpusBias_out$meta

#run selectModel for different numbers of topics (3-10)
#covariates: time (by month since 01/01/1990), source, political alignment
PrintCorpusBias_snapStorage <- manyTopics(PrintCorpusBias_out$documents, PrintCorpusBias_out$vocab, K=c(3:10),
                                          prevalence=~nmonth + majorYN + source2 + avg_align,
                                          data=PrintCorpusBias_meta,seed=1234)
save(PrintCorpusBias_snapStorage, file=paste("PrintCorpusBias_selectModel_3to10topics",Sys.Date(),sep=""))


####
### INVESTIGATE MODEL OUTPUTS
####

###
## examine the outputs of PrintCorpusFull_selectModel_3to10topics
###
load("stm_outputs/PrintCorpusFull_selectModel_3to10topics_2019-10-11")
load("dataset_creation/PrintCorpusFull_out_2019-10-30")
load("dataset_creation/PrintCorpusFull_processed_2019-10-30")

PrintCorpusFull_snapStorage$exclusivity
PrintCorpusFull_snapStorage$semcoh
summary(PrintCorpusFull_snapStorage$out[[6]])

# 8-topic model appears to have best average semantic coherence and exclusivity 
# plotted in PrintCorpFull_comparing k_3 to 10 topics_29-Oct-2019.xlsx

full_snap8 = PrintCorpusFull_snapStorage$out[[6]]

par(bty="n",lwd=2,bg="transparent",family="avenir")
plot.STM(full_snap8,n=8,labeltype="frex",family="avenir",
         main="Topics within Print Corpus, \n 1990-2017")

full_thoughts8.1= findThoughts(full_snap8, texts=PrintCorpusFull_out$meta$documents,n=100,topics=1)$docs[[1]] #administrative notices
full_thoughts8.2= findThoughts(full_snap8, texts=PrintCorpusFull_out$meta$documents,n=100,topics=2)$docs[[1]] #poverty rates and welfare
full_thoughts8.3= findThoughts(full_snap8, texts=PrintCorpusFull_out$meta$documents,n=100,topics=3)$docs[[1]] #community notices
full_thoughts8.4= findThoughts(full_snap8, texts=PrintCorpusFull_out$meta$documents,n=100,topics=4)$docs[[1]] #federal budget
full_thoughts8.5= findThoughts(full_snap8, texts=PrintCorpusFull_out$meta$documents,n=100,topics=5)$docs[[1]] #crime reports
full_thoughts8.6= findThoughts(full_snap8, texts=PrintCorpusFull_out$meta$documents,n=100,topics=6)$docs[[1]] #national poverty and welfare
full_thoughts8.7= findThoughts(full_snap8, texts=PrintCorpusFull_out$meta$documents,n=100,topics=7)$docs[[1]] #experiences "on food stamps"
full_thoughts8.8= findThoughts(full_snap8, texts=PrintCorpusFull_out$meta$documents,n=100,topics=8)$docs[[1]] #federal budget and cuts

PrintCorpusFull_out$meta$yearmonth = 2018 - (PrintCorpusFull_out$meta$nmonth)/12
full_snap8_effect0<-estimateEffect(formula=1:8~1,
                                   stmobj= full_snap8,metadata=PrintCorpusFull_out$meta )
full_snap8_effect1<-estimateEffect(formula=1:8~s(yearmonth),
                                   stmobj= full_snap8,metadata=PrintCorpusFull_out$meta )
full_snap8_effect2<-estimateEffect(formula=1:8~s(yearmonth) + majorYN ,
                                   stmobj=full_snap8,metadata=PrintCorpusFull_out$meta)
full_snap8_effect3<-estimateEffect(formula=1:8~s(yearmonth) + source2,
                                   stmobj=full_snap8,metadata=PrintCorpusFull_out$meta)

###
## examine the outputs of PrintCorpusBias_selectModel_3to10topics
###
load("stm_outputs/PrintCorpusBias_selectModel_3to10topics2019-10-29")
load("dataset_creation/PrintCorpusBias_out_2019-10-30")
load("dataset_creation/PrintCorpus_bias_2019-10-30.RData")

PrintCorpusBias_snapStorage$exclusivity
PrintCorpusBias_snapStorage$semcoh
summary(PrintCorpusBias_snapStorage$out[[6]])

# 8-topic model appears to have best average semantic coherence and exclusivity 
# plotted in PrintCorpBias_comparing k_3 to 10 topics_29-Oct-2019.xlsx

bias_snap8 = PrintCorpusBias_snapStorage$out[[6]]

par(bty="n",lwd=2,bg="transparent",family="avenir")
plot.STM(bias_snap8,n=8,labeltype="frex",family="avenir",
         main="Topics within Print Corpus Subset, \n 1990-2017")

bias_thoughts8.1= findThoughts(bias_snap8, texts=PrintCorpusBias_out$meta$documents,n=100,topics=1)$docs[[1]] #politics
bias_thoughts8.2= findThoughts(bias_snap8, texts=PrintCorpusBias_out$meta$documents,n=100,topics=2)$docs[[1]] #benefit administration
bias_thoughts8.3= findThoughts(bias_snap8, texts=PrintCorpusBias_out$meta$documents,n=100,topics=3)$docs[[1]] #welfare reform
bias_thoughts8.4= findThoughts(bias_snap8, texts=PrintCorpusBias_out$meta$documents,n=100,topics=4)$docs[[1]] #poverty rates
bias_thoughts8.5= findThoughts(bias_snap8, texts=PrintCorpusBias_out$meta$documents,n=100,topics=5)$docs[[1]] #crime reports
bias_thoughts8.6= findThoughts(bias_snap8, texts=PrintCorpusBias_out$meta$documents,n=100,topics=6)$docs[[1]] #poverty experiences or human interest 
bias_thoughts8.7= findThoughts(bias_snap8, texts=PrintCorpusBias_out$meta$documents,n=100,topics=7)$docs[[1]] #experiences "on food stamps"
bias_thoughts8.8= findThoughts(bias_snap8, texts=PrintCorpusBias_out$meta$documents,n=100,topics=8)$docs[[1]] #federal budget and cuts

PrintCorpusBias_out$meta$yearmonth = 2018 - (PrintCorpusBias_out$meta$nmonth)/12
bias_snap8_effect0<-estimateEffect(formula=1:8~1,
                                   stmobj= bias_snap8,metadata=PrintCorpusBias_out$meta )
bias_snap8_effect1<-estimateEffect(formula=1:8~s(yearmonth),
                                   stmobj= bias_snap8,metadata=PrintCorpusBias_out$meta )
bias_snap8_effect2<-estimateEffect(formula=1:8~s(yearmonth) + source2 ,
                                   stmobj=bias_snap8,metadata=PrintCorpusBias_out$meta)
bias_snap8_effect3<-estimateEffect(formula=1:8~s(yearmonth) + avg_align,
                                   stmobj=bias_snap8,metadata=PrintCorpusBias_out$meta)

plot.estimateEffect(bias_snap8_effect3, covariate = "avg_align", topics = c(1:8),
                    model = bias_snap8, method = "continuous",
                    ci.level = 0.95, printlegend = F, axes =T,
                    ylim=c(0,0.25),
                    main="SNAP Topic Variation by Partisanship Score of Media Source \n (Print Corpus Subset, 1990-2017)",
                    xlab = "Time")
summary(bias_snap8_effect3)

###
## examine the outputs of OnlineCorpus_selectModel_3to10topics
###
load("stm_outputs/OnlineCorpus_selectModel_3to10topics_2019-10-30")
load("dataset_creation/OnlineCorpus_out_2019-10-29")
load("dataset_creation/OnlineCorpus_2019-10-30.RData")

OnlineCorpus_snapStorage$exclusivity
OnlineCorpus_snapStorage$semcoh
summary(OnlineCorpus_snapStorage$out[[6]])

# 8-topic model appears to have best average semantic coherence and exclusivity 
# plotted in PrintCorpBias_comparing k_3 to 10 topics_29-Oct-2019.xlsx

online_snap8 = OnlineCorpus_snapStorage$out[[6]]

par(bty="n",lwd=2,bg="transparent",family="avenir")
plot.STM(online_snap8,n=8,labeltype="frex",family="avenir",
         main="Topics within Print Corpus Subset, \n 1990-2017")

online_thoughts8.1= findThoughts(online_snap8, texts=OnlineCorpus_out$meta$documents,n=100,topics=1)$docs[[1]] #fraud
online_thoughts8.2= findThoughts(online_snap8, texts=OnlineCorpus_out$meta$documents,n=100,topics=2)$docs[[1]] #anti-hunger programs
online_thoughts8.3= findThoughts(online_snap8, texts=OnlineCorpus_out$meta$documents,n=100,topics=3)$docs[[1]] #SNAP shopping and retailers
online_thoughts8.4= findThoughts(online_snap8, texts=OnlineCorpus_out$meta$documents,n=100,topics=4)$docs[[1]] #work requirements
online_thoughts8.5= findThoughts(online_snap8, texts=OnlineCorpus_out$meta$documents,n=100,topics=5)$docs[[1]] #federal budget
online_thoughts8.6= findThoughts(online_snap8, texts=OnlineCorpus_out$meta$documents,n=100,topics=6)$docs[[1]] #poverty rates
online_thoughts8.7= findThoughts(online_snap8, texts=OnlineCorpus_out$meta$documents,n=100,topics=7)$docs[[1]] #eligibility and immigrants
online_thoughts8.8= findThoughts(online_snap8, texts=OnlineCorpus_out$meta$documents,n=100,topics=8)$docs[[1]] #experiences on food stamps

OnlineCorpus_out$meta$yearmonth = lubridate::year(OnlineCorpus_out$meta$Date)+(lubridate::month(OnlineCorpus_out$meta$Date)-1)/12
online_snap8_effect0<-estimateEffect(formula=1:8~1,
                                     stmobj= online_snap8,metadata=OnlineCorpus_out$meta )
online_snap8_effect1<-estimateEffect(formula=1:8~s(yearmonth),
                                     stmobj= online_snap8,metadata=OnlineCorpus_out$meta )
online_snap8_effect2<-estimateEffect(formula=1:8~s(yearmonth) + majorYN ,
                                     stmobj=online_snap8,metadata=OnlineCorpus_out$meta)
online_snap8_effect3<-estimateEffect(formula=1:8~s(yearmonth) + link,
                                     stmobj=online_snap8,metadata=OnlineCorpus_out$meta)
online_snap8_effect4<-estimateEffect(formula=1:8~s(yearmonth) + avg_align,
                                     stmobj=online_snap8,metadata=OnlineCorpus_out$meta)
online_snap8_effect5<-estimateEffect(formula=1:8~s(yearmonth) + majorYN + avg_align,
                                     stmobj=online_snap8,metadata=OnlineCorpus_out$meta)

plot.estimateEffect(online_snap8_effect4, covariate = "avg_align", topics = c(1:8),
                    model = online_snap8, method = "continuous",
                    ci.level = 0.95, printlegend = F, axes =T,#xlim=c(25,325),
                    ylim=c(0,0.25),
                    main="SNAP Topic Variation \n by Partisanship Score of Media Source (Online Corpus, 2010-2017)",
                    xlab = "Time")

###
## generate plots of topic proportions over time 
## used to generate Fig2
###
par(bty="n",lwd=2,bg="transparent",family="avenir")
plot.estimateEffect(full_snap8_effect1, covariate = "yearmonth", topics = c(1:8),
                    model = full_snap8, method = "continuous",
                    ci.level = 0, printlegend = F, axes =T,xlim=c(2012,2016),
                    ylim=c(0,0.25),
                    labeltype = "frex",
                    main="SNAP Topic Variation \n (Print Corpus, 1990-2017)",
                    xlab = "Time")
plot.estimateEffect(full_snap8_effect1, covariate = "yearmonth", topics = c(1:8),
                    model = full_snap8, method = "continuous",
                    ci.level = 0, printlegend = T, axes =T,xlim=c(2012,2016),
                    ylim=c(0,0.4),
                    labeltype = "frex",
                    main="SNAP Topic Variation \n (Print Corpus, 1990-2017)",
                    xlab = "Time")

par(bty="n",lwd=2,bg="transparent",family="avenir")
plot.estimateEffect(bias_snap8_effect2, covariate = "yearmonth", topics = c(1:8),
                    model = bias_snap8, method = "continuous",
                    ci.level = 0, printlegend = F, axes =T,xlim=c(2012,2016),
                    ylim=c(0,0.25),
                    main="SNAP Topic Variation \n (Print Corpus Subset, 1990-2017)",
                    xlab = "Time")
plot.estimateEffect(bias_snap8_effect2, covariate = "yearmonth", topics = c(1:8),
                    model = bias_snap8, method = "continuous",
                    ci.level = 0, printlegend = T, axes =T,xlim=c(2012,2016),
                    ylim=c(0,0.4),labeltype = "frex",
                    main="SNAP Topic Variation \n (Print Corpus Subset, 1990-2017)",
                    xlab = "Time")

par(bty="n",lwd=2,bg="transparent",family="avenir")
plot.estimateEffect(online_snap8_effect1, covariate = "yearmonth", topics = c(1:8),
                    model = online_snap8, method = "continuous",
                    ci.level = 0, printlegend = F, axes =T,xlim=c(2012,2016),
                    ylim=c(0,0.25),
                    main="SNAP Topic Variation \n (Online Corpus, 2010-2017)",
                    xlab = "Time")
plot.estimateEffect(online_snap8_effect1, covariate = "yearmonth", topics = c(1:8),
                    model = online_snap8, method = "continuous",
                    ci.level = 0, printlegend = T, axes =T,xlim=c(2012,2016),
                    ylim=c(0,0.4),labeltype="frex",
                    main="SNAP Topic Variation \n (Online Corpus, 2010-2017)",
                    xlab = "Time")

###
## examine the outputs of estimateEffect models for print subset and online corpora, which have political alignment scores
## data used to generate Fig3
###
summary(full_snap8_effect2)
summary(bias_snap8_effect3)
summary(online_snap8_effect5)


### 
## generate tables of average maximum topic prevalence per document
## data used to generate Table 3
###
require(dplyr)

PrintFull_scored8 = cbind(full_snap8$theta[,1],full_snap8$theta[,2], full_snap8$theta[,3],full_snap8$theta[,4],full_snap8$theta[,5],full_snap8$theta[,6],full_snap8$theta[,7],full_snap8$theta[,8])
PrintFull_scored8 = as.data.frame(PrintFull_scored8)
PrintFull_scored8$max_topic = max.col(PrintFull_scored8)
PrintFull_scored8$max_topic_theta <- do.call(pmax, PrintFull_scored8[1:8])
PrintFull_scored8 = cbind(PrintCorpusFull_out$meta,PrintFull_scored8)

PrintBias_scored8 = cbind(bias_snap8$theta[,1],bias_snap8$theta[,2], bias_snap8$theta[,3],bias_snap8$theta[,4],bias_snap8$theta[,5],bias_snap8$theta[,6],bias_snap8$theta[,7],bias_snap8$theta[,8])
PrintBias_scored8 = as.data.frame(PrintBias_scored8)
PrintBias_scored8$max_topic = max.col(PrintBias_scored8)
PrintBias_scored8$max_topic_theta <- do.call(pmax, PrintBias_scored8[1:8])
PrintBias_scored8 = cbind(PrintCorpusBias_out$meta,PrintBias_scored8)

Online_scored8 = cbind(online_snap8$theta[,1],online_snap8$theta[,2], online_snap8$theta[,3],online_snap8$theta[,4],online_snap8$theta[,5],online_snap8$theta[,6],online_snap8$theta[,7],online_snap8$theta[,8])
Online_scored8 = as.data.frame(Online_scored8)
Online_scored8$max_topic = max.col(Online_scored8)
Online_scored8$max_topic_theta <- do.call(pmax, Online_scored8[1:8])
Online_scored8 = cbind(OnlineCorpus_out$meta,Online_scored8)

require(dplyr)
maxtable=
  rbind(
    cbind(PrintFull_scored8 %>% group_by(max_topic) %>% summarise(n(),mean(max_topic_theta),sd(max_topic_theta)), title=c(rep("printFull_max",8))),
    cbind(PrintFull_scored8[PrintFull_scored8$max_topic_theta>0.50,] %>% group_by(max_topic) %>% summarise(n(),mean(max_topic_theta),sd(max_topic_theta)), title=c(rep("printFull_maxGT50",8))),
    cbind(PrintBias_scored8 %>% group_by(max_topic) %>% summarise(n(),mean(max_topic_theta),sd(max_topic_theta)), title=c(rep("printBias_max",8))),
    cbind(PrintBias_scored8[PrintBias_scored8$max_topic_theta>0.50,] %>% group_by(max_topic) %>% summarise(n(),mean(max_topic_theta),sd(max_topic_theta)), title=c(rep("printBias_maxGT50",8))),
    cbind(Online_scored8 %>% group_by(max_topic) %>% summarise(n(),mean(max_topic_theta),sd(max_topic_theta)), title=c(rep("online_max",8))),
    cbind(Online_scored8[Online_scored8$max_topic_theta>0.50,] %>% group_by(max_topic) %>% summarise(n(),mean(max_topic_theta),sd(max_topic_theta)), title=c(rep("online_maxGT50",8)))
  )

#write output to csv for Table 3
write.csv(maxtable,"Table3_data.csv")

#redact print media corpora 
PrintFull_scored8$documents = stringr::str_sub(PrintFull_scored8$documents,1,50)
PrintBias_scored8$documents = stringr::str_sub(PrintBias_scored8$documents,1,50)

#write scored datasets to csv for replication dataset
write.csv(Online_scored8,"OnlineCorpus_Scored.csv")
write.csv(PrintFull_scored8,"PrintCorpus_Scored_redacted.csv")
write.csv(PrintBias_scored8,"PrintSubsetCorpus_Scored_redacted.csv")


###
## generate outputs for Table 1
###
summary(full_snap8)
summary(bias_snap8)
summary(online_snap8)

PrintFull_scored8$major_paper = ifelse(PrintFull_scored8$majorYN=="major_paper",1,0)
PrintBias_scored8$major_paper = ifelse(PrintBias_scored8$majorYN=="major_paper",1,0)
Online_scored8$major_paper = ifelse(Online_scored8$majorYN=="1",1,0)

table1=rbind(
  cbind(
    rbind(
      cbind(PrintFull_scored8 %>% group_by(max_topic) %>% summarise(n(),mean(yearmonth),sd(yearmonth),sum(major_paper)),
            "mean(avg_align)"=NA,"sd(avg_align)"=NA),
      cbind(max_topic="Overall",
            PrintFull_scored8 %>% summarise(n(),mean(yearmonth),sd(yearmonth),sum(major_paper)),
            "mean(avg_align)"=NA,"sd(avg_align)"=NA)),
    corpus="Print"),
  cbind(
    rbind(PrintBias_scored8 %>% group_by(max_topic) %>% summarise(n(),mean(yearmonth),sd(yearmonth),sum(major_paper),mean(avg_align),sd(avg_align)),
          cbind(max_topic="Overall",
                PrintBias_scored8 %>% summarise(n(),mean(yearmonth),sd(yearmonth),sum(major_paper),mean(avg_align),sd(avg_align)))),
    corpus="Print Subset"),
  cbind(
    rbind(Online_scored8 %>% group_by(max_topic) %>% summarise(n(),mean(yearmonth),sd(yearmonth),sum(major_paper),mean(avg_align),sd(avg_align)),
          cbind(max_topic="Overall",
                Online_scored8 %>% summarise(n(),mean(yearmonth),sd(yearmonth),sum(major_paper),mean(avg_align),sd(avg_align)))),
    corpus="Online")
)
write.csv(table1,"Table1_data.csv",row.names = FALSE)




