---
editor_options:
  chunk_output_type: console
fontsize: 12pt
geometry: margin=1in
header-includes: \usepackage{float} \usepackage{lineno} \usepackage{setspace}
  \usepackage[round]{natbib} \bibpunct[; ]{(}{)}{,}{a}{}{,} \usepackage{color} \usepackage{totcount} \usepackage[normalem]{ulem}
  \newtotcounter{citenum} \def\oldcite{} \let\oldcite=\bibcite \def\bibcite{\stepcounter{citenum}\oldcite}
  \usepackage{fancyhdr} \pagestyle{fancy} \fancyhf{} \fancyfoot[LE,LO]{\textcolor{red}{Preprint
  - This work has not yet been peer-reviewed}} \fancyfoot[RE,RO]{\thepage} \renewcommand{\headrulewidth}{0pt}
output:
  pdf_document:
    fig_caption: true
    keep_tex: yes
    number_sections: no
  word_document: default
  html_document:
    df_print: paged
---




\begin{center}
\textbf{\Large If this title is funny, will you cite me? \\
Citation impacts of humour\sout{, cultural references,} and other features of article titles in ecology and evolution
}
\vspace{5 mm}
	
\textsc{Stephen B. Heard$^{1,2}$, Chloe A. Cull$^{1}$, and Easton R. White$^{4}$}
\vspace{3 mm}

\small{$^{1}$Dept. of Biology, University of New Brunswick, Fredericton, NB Canada E3B 5A3 \\

$^{2}$Corresponding author: sheard@unb.ca \\

$^{3}$Permanent address: **** \\

$^{4}$Department of Biological Sciences, University of New Hampshire, Durham, NH, USA, 03824
}

\end{center}

\vspace{3 mm}

<!--\tableofcontents-->



```{r load_packages,echo=F,warning=F,message=F}
if (!require("pacman",character.only = TRUE))
  {
    install.packages("pacman",dep=TRUE)
    if(!require("pacman",character.only = TRUE)) stop("Package not found")
  }

# Keeping below source for github package. Ask Easton whether pacman works for github packages or not.
#devtools::install_github("rensa/stickylabeller")
pacman::p_load(dplyr, tidyr, ggplot2,MASS,stargazer,stringr,irr, corrplot,rstatix,ggpubr)
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE,message = FALSE,warning = FALSE,dev = c("png"),
                      dpi = 300,
                      echo = FALSE,
                      cache = TRUE)
```



# Methods - Data Analysis 

We use a series of generalized linear models (each with a Poisson error structure) to assess the effect of various predictors on both total citation count and self-citation count. Specifically, we examined the effect of article type (PrimaryReviewOther), whether the title was phrased as a question (Question), whether the title was assertive (Assertive), the presence of a colon in the title (Colon), the presence of any acronyms (Acronyms), whether the political region was noted in the title (Location), the presence of a taxonomic name (Taxonomic_name), and the average humour score for each title (avg_humour). We assessed each combination of these predictor variables and ranked models according to AIC criteria. We present only the best fitting model for each response variable in the results section. We did not include highly-correlated (>0.7) predictor variables in the same model. We examined residual plots to verify model assumptions were met.

\pagebreak

# Results 

```{r article_type,echo=F,message=F,warnings=F,fig.cap='Total citations and self citations versus article type (Other, Primary, or Review). The article types were significantly different (generalized linear model with a Poisson error term for counts and a Gaussian error term for average humor) from one another.   \\label{fig:article_type}',echo=F,fig.height=8}

svg("figure1.svg",width = 6.5,height=8)


titles_minimal <- read.csv(file='Data/Cleaned_Title_Database_no_scorers.csv',header=T)

par(mfrow=c(2,2),mar=c(3,5,1,1),oma=c(6,1,1,1))
plot(titles_minimal$Total_citations~as.factor(titles_minimal$PrimaryReviewOther),pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Total citations',xlab=' ',cex.lab=1.5)
#mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(humour_total_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(humour_total_glm ))[2,4],2)) ), side = 3)

plot(titles_minimal$Self_citations~as.factor(titles_minimal$PrimaryReviewOther),pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Self citations',xlab=' ',cex.lab=1.5)

#mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(humour_self_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(humour_self_glm ))[2,4],2)) ), side = 3)

titles_long <- read.csv(file='Data/Cleaned_Title_Database_with_scorers.csv',header=T)

humour <- titles_long %>%
  filter(metric == 'humour') %>%
  mutate(score = as.numeric(score)) %>%
  group_by(Index) %>%
  mutate(avg_humour = mean(score,na.rm=T)) %>%
  distinct(Index, .keep_all=TRUE)

titles2 <- left_join(titles_minimal,humour)


plot(titles2$avg_humour~as.factor(titles_minimal$PrimaryReviewOther),pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Average humor',xlab=' ',cex.lab=1.5)
#mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(humour_total_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(humour_total_glm ))[2,4],2)) ), side = 3)

require(stringr)
titles2$Journal <- str_trim(titles2$Journal, side = c("both"))
titles2$Journal[titles2$Journal=='OIkos'] <- 'Oikos'


par(mar=c(10,5,1,1))

plot(titles2$avg_humour~as.factor(titles2$Journal),pch=16,col=rgb(0.5,0.5,0.5,0.3),las=2,ylab='Average humor',xlab=' ',cex.lab=1.5, pars  =  list(xaxt = "n"))
axis(1, at=c(1:9), labels = FALSE)
text(c(1:9 +1), par("usr")[3] - 1, labels = names(table(as.factor(titles2$Journal))), srt = 45, pos = 2, xpd = TRUE)

m1 <- glm(titles2$Total_citations ~ titles2$PrimaryReviewOther, family = 'poisson')
m2 <- glm(titles2$Self_citations ~ titles2$PrimaryReviewOther, family = 'poisson')
m3 <- lm(titles2$avg_humour ~ titles2$PrimaryReviewOther)
m4 <- lm(titles2$avg_humour ~ titles2$Journal)

#GLM.GAMMA <- glm(avg_humour ~ Journal,data=titles2, family="Gamma"(link='log'))
#mtext('Article type', side=1, line=-2,outer=T,cex = 1.5)

dev.off()
```

```{r}
require(knitr)
kable(data.frame(table(titles2$Journal,titles2$PrimaryReviewOther)))

```


Citation counts for the papers we tracked were extremely variable, ranging from (sadly) zero to `r max(titles_minimal$Total_citations)`. Unsurprisingly, review papers were cited more heavily, on average, than primary research papers (Fig. \ref{fig:article_type}). Total citations increased by an average of  `r 100*round(abs(1-exp(as.numeric(coef(m1)[2]))),2)` and `r 100*round(abs(1-exp(as.numeric(coef(m1)[3]))),2)` for Primary and Review articles compared to Other article types. Simiarly, self citations increased by an average of  `r 100*round(abs(1-exp(as.numeric(coef(m2)[2]))),2)` and `r 100*round(abs(1-exp(as.numeric(coef(m2)[3]))),2)` for Primary and Review articles compared to Other article types. All of these differences were significant (p<0.0001). Other articles types were rated significantly more humourous than Review and Primary articles. Humour did not vary significantly by journal except for Trends in Ecology and Evolution, which had a `r as.numeric(round(coef(m4)[length(coef(m4))],2))` increase in humour compared to other journals. Our best fitting models and parameter estimates were similar whether we analyzed all articles or just primary research papers (Tables \ref{tab:best_fit_model_total}-\ref{tab:best_fit_model_self_primary}).


\pagebreak

```{r cites_vs_title_detail,echo=F,message=F,warnings=F,fig.cap='Total and self citations for article where the title included an acronym or taxonomic name. The exp($\\beta$) above each plot represents the change in citations for the presence of an acronym or taxonomic name. For example, exp($\\beta$) = 0.55 represents a 45% decline in citations for the presence of an acronym. This is for a generalized linear model with a Poisson error term with acronym or taxonomic name as the only predictor variable. Multivariate results are available in tables S1-S4. \\label{fig:presence}',echo=F,fig.height=6}

titles_minimal <- read.csv(file='Data/Cleaned_Title_Database_no_scorers.csv',header=T)



#Models
glm_total_Acronyms <- glm(formula = Total_citations~Acronyms,data = titles_minimal,family = 'poisson')
glm_total_Taxonomic_name <- glm(formula = Total_citations~Taxonomic_name,data = titles_minimal,family = 'poisson')
#glm_total_Location <- glm(formula = Total_citations~Location,data = titles_minimal,family = 'poisson')

#Plots
par(mfrow=c(2,2),mar=c(4,4.5,2.5,1))

boxplot(titles_minimal$Total_citations~titles_minimal$Acronyms,ylim=c(0,500),xlab='Acronyms',ylab='Total citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_total_Acronyms)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_total_Acronyms))[2,4],2)) ), side = 3)

boxplot(titles_minimal$Total_citations~titles_minimal$Taxonomic_name,ylim=c(0,500),xlab='Taxonomic name',ylab='Total citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_total_Taxonomic_name)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_total_Taxonomic_name))[2,4],2)) ), side = 3)

#boxplot(titles_minimal$Total_citations~titles_minimal$Location,ylim=c(0,500),xlab='Political region noted',ylab='Total citations',cex.lab=1.5,las=1)
#mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_total_Location)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_total_Location))[2,4],2)) ), side = 3)



# Models
glm_self_Acronyms <- glm(formula = Self_citations~Acronyms,data = titles_minimal,family = 'poisson')
glm_self_Taxonomic_name <- glm(formula = Self_citations~Taxonomic_name,data = titles_minimal,family = 'poisson')
#glm_self_Location <- glm(formula = Self_citations~Location,data = titles_minimal,family = 'poisson')

# Plots
boxplot(titles_minimal$Self_citations~titles_minimal$Acronyms,ylim=c(0,60),xlab='Acronyms',ylab='Self citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_self_Acronyms)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_self_Acronyms))[2,4],2)) ), side = 3)

boxplot(titles_minimal$Self_citations~titles_minimal$Taxonomic_name,ylim=c(0,60),xlab='Taxonomic name',ylab='Self citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_self_Taxonomic_name)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_self_Taxonomic_name))[2,4],2)) ), side = 3)

#boxplot(titles_minimal$Self_citations~titles_minimal$Location,ylim=c(0,60),xlab='Political region noted',ylab='Self citations',cex.lab=1.5,las=1)
#mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_self_Location)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_self_Location))[2,4],2)) ), side = 3)


```


We found that several characteristics of a paper title were significant predictors of the paper's citation count (Figs. \ref{fig:presence}, Tables \ref{tab:best_fit_model_total}, \ref{tab:best_fit_model_self}). However, most of these predictors had relatively small effects on both total and self citation counts (Figs. \ref{fig:presence},\ref{fig:basic_metrics}). For example, we found that the presence of an acronym or taxonomic name decreased the total citation count by `r 100-100*as.numeric(round(exp(coef(glm_total_Acronyms)[2]),2))`% and `r 100-100*as.numeric(round(exp(coef(glm_total_Taxonomic_name)[2]),2))`%, respectively (Fig. \ref{fig:presence}; values were also similar in the multivariate analyses in Tables \ref{tab:best_fit_model_total}, \ref{tab:best_fit_model_self}).

We found that the presence of a colon or political region increased total citations (Fig. \ref{fig:basic_metrics}, Tables \ref{tab:best_fit_model_total}, \ref{tab:best_fit_model_self}). We also found that a title phased as an assertive statement or as a question also increased the total citation counts (Fig. \ref{fig:basic_metrics}, Tables \ref{tab:best_fit_model_total}, \ref{tab:best_fit_model_self}). Self citations generally followed a similar pattern to total citations (Tables \ref{tab:best_fit_model_total}, \ref{tab:best_fit_model_self}). However, for some variables, such as the presence of an acronym or taxonomic name, they had a significant effect for total citations, but not self citations (Tables \ref{tab:best_fit_model_total}, \ref{tab:best_fit_model_self}).



\pagebreak

## Humor

```{r humor,echo=F,message=F,warnings=F,fig.cap='Both total citations and self citations signficantly decreased with higher humour scores. However, the effect size for total citations was very small compared to self citations.\\label{fig:humour}',echo=F}

svg("figure3.svg",width = 8,height=4)


titles_long <- read.csv(file='Data/Cleaned_Title_Database_with_scorers.csv',header=T)

humour <- titles_long %>%
  filter(metric == 'humour') %>%
  mutate(score = as.numeric(score)) %>%
  group_by(Index) %>%
  mutate(avg_humour = mean(score,na.rm=T)) %>%
  distinct(Index, .keep_all=TRUE)

par(mfrow=c(1,3))
plot(humour$Total_citations~humour$avg_humour,pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Total citations',xlab='Average humour score',cex.lab=1.5)
humour_total_glm <- glm(formula = Total_citations~avg_humour,data = humour,family = 'poisson')
newdata <- expand.grid(avg_humour=seq(min(humour$avg_humour), max(humour$avg_humour), length.out=200))
newdata <- cbind(newdata, citations=predict(humour_total_glm, newdata, type='response'))
points(newdata,type='l',lwd=4)
#mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(humour_total_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(humour_total_glm ))[2,4],2)) ), side = 3)

plot(humour$Self_citations~humour$avg_humour,pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Self citations',xlab='Average humour score',cex.lab=1.5)
humour_self_glm <- glm(formula = Self_citations~avg_humour,data = humour,family = 'poisson')
newdata <- expand.grid(avg_humour=seq(min(humour$avg_humour), max(humour$avg_humour), length.out=200))
newdata <- cbind(newdata, citations=predict(humour_self_glm, newdata, type='response'))
points(newdata,type='l',lwd=4)
#mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(humour_self_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(humour_self_glm ))[2,4],2)) ), side = 3)


humour$total_vs_self  <-  round((humour$Total_citations - humour$Self_citations)/humour$Self_citations)
humour <- humour[-which(humour$total_vs_self=='Inf'),]
plot(humour$total_vs_self~humour$avg_humour,pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Total/ self',xlab='Average humour score',cex.lab=1.5)
humour_total_vs_self_glm <- glm(formula = total_vs_self~avg_humour,data = humour,family = 'poisson')
newdata <- expand.grid(avg_humour=seq(min(humour$avg_humour), max(humour$avg_humour), length.out=200))
newdata <- cbind(newdata, citations=predict(humour_total_vs_self_glm , newdata, type='response'))
points(newdata,type='l',lwd=4)
#mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(humour_total_vs_self_glm  )[2])),2)) ~ "," ~ p == .(round(coef(summary(humour_total_vs_self_glm))[2,4],2)) ), side = 3)
dev.off()
```

After controlling for other predictor variables (Tables \ref{tab:best_fit_model_total}, \ref{tab:best_fit_model_self}), the average title humor score was still a significant predictor of both total and self citations. The effect was relatively small with a decrease of `r 100-100*as.numeric(round(exp(coef(humour_total_glm)[2]),2))`% and `r 100-100*as.numeric(round(exp(coef(humour_self_glm)[2]),2))`%, respectively. percent in total citations for each 1 point increase in the average humor score. However, this equates to a difference of `r 100*(1-round(exp(5*as.numeric(coef(humour_total_glm)[2])),4))`% between the least and most humorous titles. The effect was larger for self citations with a `r 100*(1-round(exp(5*as.numeric(coef(humour_self_glm)[2])),4))`% decrease for the most humorous titles. These results were similar between these univariate analyses (Fig. \ref{fig:humour}) and after our model selection process (Tables \ref{tab:best_fit_model_total}, \ref{tab:best_fit_model_self}). We also see that when we corrected for "author perceived self-importance", by dividing total citations by self, we see humour actually increases citations by `r 100*abs(1-as.numeric(round(exp(coef(humour_total_vs_self_glm)[2]),2)))`%.





## Papers with the highest rated funny titles (average humour scores >3) include:

- `r paste(humour$Title[humour$avg_humour>3][1])`  
- `r paste(humour$Title[humour$avg_humour>3][2])`  
- `r paste(humour$Title[humour$avg_humour>3][3])`  
- `r paste(humour$Title[humour$avg_humour>3][4])`  
- `r paste(humour$Title[humour$avg_humour>3][5])`  
- `r paste(humour$Title[humour$avg_humour>3][6])`  
- `r paste(humour$Title[humour$avg_humour>3][7])`  
- `r paste(humour$Title[humour$avg_humour>3][8])`  
- `r paste(humour$Title[humour$avg_humour>3][9])`  





\clearpage



<!--

# Main findings

- With a proper generalized linear model (and a Poisson error structure) there are several significant predicators of total and self citations. *This is a big change from earlier analyses.*
- Total citations increased significantly with a colon, a political area, assertive titles, or question titles (Figs. \ref{fig:basic_metrics},\ref{fig:presence}). However, total citations significantly decreased with the presence of acronyms or scientific names (Fig. \ref{fig:presence}).  
- Self citations increased significantly with the inclusion of a colon, scientific name, or political region (Figs. \ref{fig:basic_metrics},\ref{fig:presence}). However, total citations significantly decreased with the presence of an question title, (Fig. \ref{fig:basic_metrics}).  
- Both total citations and self citations significantly decreased with higher humour scores. However, the effect size for total citations was very small compared to self citations (Fig. \ref{fig:humour}).  
- Total citations significantly increased with cultural score while self citations significantly decreased with average cultural score (Fig. \ref{fig:cultural}).  
- Both total and self citations descreased significantly with higher offensive scores, but the effect sizes were small (Fig. \ref{fig:offensive}).  - There was not a lot of concordance between reviewers (Fig. \ref{fig:concordance})   

-->






\pagebreak





\pagebreak


# SM

\begin{center}
\textbf{\Large Supplemental material for: \\ If this title is funny, will you cite me? Citation impacts of humour, cultural references, and other features of article titles in ecology and evolution
}
\vspace{5 mm}
	
\textsc{Stephen B. Heard$^{1,2}$, Chloe A. Cull$^{1}$, and Easton R. White$^{4}$}
\vspace{3 mm}

\small{$^{1}$Dept. of Biology, University of New Brunswick, Fredericton, NB Canada E3B 5A3 \\

$^{2}$Corresponding author: sheard@unb.ca \\

$^{3}$Permanent address: **** \\

$^{4}$Department of Biological Sciences, University of New Hampshire, Durham, NH, USA, 03824
}

\end{center}

\vspace{3 mm}

\renewcommand{\thepage}{S\arabic{page}} 
\renewcommand{\thesection}{S\arabic{section}}  
\renewcommand{\thetable}{S\arabic{table}}  
\renewcommand{\thefigure}{S\arabic{figure}}
\setcounter{table}{0}
\setcounter{figure}{0}
\setcounter{section}{0}
\setcounter{page}{0}

\pagebreak

# Extra findings

```{r cites_vs_title_metrics,echo=F,message=F,warnings=F,fig.cap='See after figure for findings \\label{fig:basic_metrics}',echo=F,fig.width=11,fig.height=8,out.extra='angle=90'}

titles_minimal <- read.csv(file='Data/Cleaned_Title_Database_no_scorers.csv',header=T)
  
# Tests
glm_total_colon <- glm(formula = Total_citations~Colon,data = titles_minimal,family = 'poisson')
glm_total_question <- glm(formula = Total_citations~Question,data = titles_minimal,family = 'poisson')
glm_total_assertive <- glm(formula = Total_citations~Assertive,data = titles_minimal,family = 'poisson')
glm_total_Location <- glm(formula = Total_citations~Location,data = titles_minimal,family = 'poisson')


# Plots
par(mfrow=c(2,4),mar=c(4,4.5,2.5,0.5))

boxplot(titles_minimal$Total_citations~titles_minimal$Colon,ylim=c(0,500),xlab='Colon',ylab='Total citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_total_colon)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_total_colon))[2,4],2)) ), side = 3)

boxplot(titles_minimal$Total_citations~titles_minimal$Question,ylim=c(0,500),xlab='Question',ylab='Total citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_total_question)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_total_question))[2,4],2)) ), side = 3)

boxplot(titles_minimal$Total_citations~titles_minimal$Assertive,ylim=c(0,500),xlab='Assertive',ylab='Total citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_total_assertive)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_total_assertive))[2,4],2)) ), side = 3)

boxplot(titles_minimal$Total_citations~titles_minimal$Location,ylim=c(0,500),xlab='Political region noted',ylab='Total citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_total_Location)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_total_Location))[2,4],2)) ), side = 3)

#Models
glm_self_colon <- glm(formula = Self_citations~Colon,data = titles_minimal,family = 'poisson')
glm_self_question <- glm(formula = Self_citations~Question,data = titles_minimal,family = 'poisson')
glm_self_assertive <- glm(formula = Self_citations~Assertive,data = titles_minimal,family = 'poisson')
glm_self_Location <- glm(formula = Self_citations~Location,data = titles_minimal,family = 'poisson')

# Plots
boxplot(titles_minimal$Self_citations~titles_minimal$Colon,ylim=c(0,60),xlab='Colon',ylab='Self citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_self_colon)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_self_colon))[2,4],2)) ), side = 3)
boxplot(titles_minimal$Self_citations~titles_minimal$Question,ylim=c(0,60),xlab='Question',ylab='Self citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_self_question)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_self_question))[2,4],2)) ), side = 3)
boxplot(titles_minimal$Self_citations~titles_minimal$Assertive,ylim=c(0,60),xlab='Assertive',ylab='Self citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_self_assertive)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_self_assertive))[2,4],2)) ), side = 3)

boxplot(titles_minimal$Self_citations~titles_minimal$Location,ylim=c(0,60),xlab='Political region noted',ylab='Self citations',cex.lab=1.5,las=1)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(glm_self_Location)[2])),2)) ~ "," ~ p == .(round(coef(summary(glm_self_Location))[2,4],2)) ), side = 3)
```


\pagebreak

# All article types

```{r full_model_total_citations,eval=T,echo=F,results='asis',warning=FALSE,message=FALSE,error=FALSE}

titles_minimal <- read.csv(file='Data/Cleaned_Title_Database_no_scorers.csv',header=T)

titles_long <- read.csv(file='Data/Cleaned_Title_Database_with_scorers.csv',header=T)

humour <- titles_long %>%
  filter(metric == 'humour') %>%
  mutate(score = as.numeric(score)) %>%
  group_by(Index) %>%
  mutate(avg_humour = mean(score,na.rm=T)) %>%
  distinct(Index, .keep_all=TRUE)

titles2 <- left_join(titles_minimal,humour)
#titles2 <- left_join(titles2,cultural)
#titles2 <- left_join(titles2,offensive)

# Should do model selection here...
mymod <- glm(formula = Total_citations ~ PrimaryReviewOther + Question + Assertive + Colon + Acronyms + Location + Taxonomic_name + avg_humour,data=titles2, family = 'poisson')

top_total_citations_mod <- stepAIC(mymod,direction='both',trace = 0)
#summary(mymod)

stargazer(top_total_citations_mod, type = "latex", title="Best fitting model after AIC model selection for total citations.\\label{tab:best_fit_model_total}",header=F)# dep.var.labels=c("Early doubling time","Overall doubling time"),covariate.labels=c("Restrict restaurants","log (Population density)","Vaccination rate","GNI per capita","Population percent in rural areas","Tightness score","Constant"),omit.stat=c("LL","ser","f"), no.space=TRUE,ci=TRUE,model.numbers = FALSE, notes = "",header=FALSE)

```

```{r full_model_self_citations,eval=T,echo=F,results='asis',warning=FALSE,message=FALSE,error=FALSE}

titles_minimal <- read.csv(file='Data/Cleaned_Title_Database_no_scorers.csv',header=T)

titles2 <- left_join(titles_minimal,humour)
#titles2 <- left_join(titles2,cultural)
#titles2 <- left_join(titles2,offensive)

# Should do model selection here...
mymod <- glm(formula = Self_citations ~ PrimaryReviewOther + Question + Assertive + Colon + Acronyms + Location + Taxonomic_name + avg_humour,data=titles2, family = 'poisson')

top_self_citations_mod <- stepAIC(mymod,direction='both',trace = 0)

#summary(mymod)

stargazer(top_self_citations_mod, type = "latex", title="Best fitting model after AIC model selection for self citations.\\label{tab:best_fit_model_self}",header = F)# dep.var.labels=c("Early doubling time","Overall doubling time"),covariate.labels=c("Restrict restaurants","log (Population density)","Vaccination rate","GNI per capita","Population percent in rural areas","Tightness score","Constant"),omit.stat=c("LL","ser","f"), no.space=TRUE,ci=TRUE,model.numbers = FALSE, notes = "",header=FALSE)

```

\pagebreak

# Only Primary types

```{r full_model_total_citations_Primary,eval=T,echo=F,results='asis',warning=FALSE,message=FALSE,error=FALSE}

titles_minimal <- read.csv(file='Data/Cleaned_Title_Database_no_scorers.csv',header=T)

titles_long <- read.csv(file='Data/Cleaned_Title_Database_with_scorers.csv',header=T)

humour <- titles_long %>%
  filter(metric == 'humour') %>%
  mutate(score = as.numeric(score)) %>%
  group_by(Index) %>%
  mutate(avg_humour = mean(score,na.rm=T)) %>%
  distinct(Index, .keep_all=TRUE)

titles2 <- left_join(titles_minimal,humour)
#titles2 <- left_join(titles2,cultural)
#titles2 <- left_join(titles2,offensive)

# Should do model selection here...
mymod <- glm(formula = Total_citations ~ Question + Assertive + Colon + Acronyms + Location + Taxonomic_name + avg_humour,data=titles2 %>%filter(PrimaryReviewOther =='Primary'), family = 'poisson')

top_total_citations_mod <- stepAIC(mymod,direction='both',trace = 0)
#summary(mymod)

stargazer(top_total_citations_mod, type = "latex", title="Best fitting model after AIC model selection for total citations for PRIMARY articles only.\\label{tab:best_fit_model_total_primary}",header=F)# dep.var.labels=c("Early doubling time","Overall doubling time"),covariate.labels=c("Restrict restaurants","log (Population density)","Vaccination rate","GNI per capita","Population percent in rural areas","Tightness score","Constant"),omit.stat=c("LL","ser","f"), no.space=TRUE,ci=TRUE,model.numbers = FALSE, notes = "",header=FALSE)

```


```{r full_model_self_citations_Primary,eval=T,echo=F,results='asis',warning=FALSE,message=FALSE,error=FALSE}

titles_minimal <- read.csv(file='Data/Cleaned_Title_Database_no_scorers.csv',header=T)

titles2 <- left_join(titles_minimal,humour)
#titles2 <- left_join(titles2,cultural)
#titles2 <- left_join(titles2,offensive)

# Should do model selection here...
mymod <- glm(formula = Self_citations ~  Question + Assertive + Colon + Acronyms + Location + Taxonomic_name + avg_humour,data=titles2 %>%filter(PrimaryReviewOther =='Primary'), family = 'poisson')

top_self_citations_mod <- stepAIC(mymod,direction='both',trace = 0)

#summary(mymod)

stargazer(top_self_citations_mod, type = "latex", title="Best fitting model after AIC model selection for self citations for PRIMARY articles only.\\label{tab:best_fit_model_self_primary}",header = F)# dep.var.labels=c("Early doubling time","Overall doubling time"),covariate.labels=c("Restrict restaurants","log (Population density)","Vaccination rate","GNI per capita","Population percent in rural areas","Tightness score","Constant"),omit.stat=c("LL","ser","f"), no.space=TRUE,ci=TRUE,model.numbers = FALSE, notes = "",header=FALSE)

```

\pagebreak

# Concordance measures


```{r concordance,echo=F,eval=T,fig.width=5,fig.cap='Correlation coefficients (which are almost always higher than concordance values) between each pair of scorers.'}

titles_long <- read.csv(file='Data/Cleaned_Title_Database_with_scorers.csv',header=T)


titles2 <- titles_long %>%
  select(Index:score) %>%
  filter(score %in% c(0,1,2,3,4,5,6)) %>%
  filter(metric=='humour') %>%
  mutate(score = as.numeric(score))  %>%
  mutate(score = ifelse(score>0,1,0)) %>% # Optional line to make binary variables
  tidyr::spread(key=name,value=score) %>%
  dplyr::select(-c(Index:metric))

overall_humor_concordance <- kappam.light(titles2)
#kappam.fleiss(titles2,detail=T)
#titles2$Easton <- sample(0:3,size=nrow(titles2),replace = T)


par(mfrow=c(1,1),mar=c(0,0,0,0),oma=c(0.5,0.5,0,0))

corrplot(cor(na.omit(titles2)),type = 'lower',method='shade',diag=T,number.cex  = 0.8)
mtext('Humour metric',side=3,line=-3,adj= 0.75)

write.csv(x=cor(na.omit(titles2)),file='scorer_correlations.csv',quote=F)


```

It looks like there is not high concordance among the various scorers. There was a concordance value of `r overall_humor_concordance$value` on a scale from 0 to 1, with 1 being perfect concordance. This is even the case when scores are converted to binary variables, e.g., funny or not.

\clearpage


# Title length 

Not very interesting results in terms of title length. Average humour dropped with longer titles. Total and self citations were relatively flat with title length. 

```{r, echo=F}


titles_minimal <- read.csv(file='Data/Cleaned_Title_Database_no_scorers.csv',header=T)

titles_minimal$title_length <- str_count(titles_minimal$Title)

titles_long <- read.csv(file='Data/Cleaned_Title_Database_with_scorers.csv',header=T)

humour <- titles_long %>%
  filter(metric == 'humour') %>%
  mutate(score = as.numeric(score)) %>%
  group_by(Index) %>%
  mutate(avg_humour = mean(score,na.rm=T)) %>%
  distinct(Index, .keep_all=TRUE)

titles2 <- left_join(titles_minimal,humour)

ggplot(data= titles2, aes(x=title_length,y=avg_humour)) + geom_point(alpha=0.3) + geom_smooth(method = 'loess') + ylim(0,5)



ggplot(data= titles_minimal, aes(x=title_length,y=Total_citations)) + geom_point(alpha=0.3) + geom_smooth(method = 'loess') + ylim(0,500)


ggplot(data= titles_minimal, aes(x=title_length,y=Self_citations)) + geom_point(alpha=0.3) + geom_smooth(method = 'loess')
```


\clearpage

# Cultural scores

```{r,echo=F,eval=T,message=F,warnings=F,fig.cap='Total citations significantly increased with cultural score while self citations significanly decreased with average cultural score..\\label{fig:cultural}'}

titles_long <- read.csv(file='Data/Cleaned_Title_Database_with_scorers.csv',header=T)

cultural <- titles_long %>%
  filter(score != 'racist wordings') %>%
  filter(metric == 'cultural') %>%
  mutate(score = as.numeric(score)) %>%
  group_by(Index) %>%
  mutate(avg_cultural = mean(score,na.rm=T)) %>%
  distinct(Index, .keep_all=TRUE)


par(mfrow=c(1,2))
plot(cultural$Total_citations~cultural$avg_cultural,pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Total citations',xlab='Average cultural score',cex.lab=1.5)
cultural_total_glm <- glm(formula = Total_citations~avg_cultural,data = cultural,family = 'poisson')
newdata <- expand.grid(avg_cultural=seq(min(cultural$avg_cultural), max(cultural$avg_cultural), length.out=200))
newdata <- cbind(newdata, citations=predict(cultural_total_glm, newdata, type='response'))
points(newdata,type='l',lwd=4)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(cultural_total_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(cultural_total_glm ))[2,4],2)) ), side = 3)

plot(cultural$Self_citations~cultural$avg_cultural,pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Self citations',xlab='Average cultural score',cex.lab=1.5)
cultural_self_glm <- glm(formula = Self_citations~avg_cultural,data = cultural,family = 'poisson')
newdata <- expand.grid(avg_cultural=seq(min(cultural$avg_cultural), max(cultural$avg_cultural), length.out=200))
newdata <- cbind(newdata, citations=predict(cultural_self_glm, newdata, type='response'))
points(newdata,type='l',lwd=4)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(cultural_self_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(cultural_self_glm ))[2,4],2)) ), side = 3)



```

## Papers with the highest rated culturally-relevant titles (average humour scores >0.4) include:

- `r paste(cultural$Title[cultural$avg_cultural>0.4][1])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][2])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][3])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][4])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][5])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][6])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][7])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][8])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][9])`
- `r paste(cultural$Title[cultural$avg_cultural>0.4][10])`

\pagebreak




# Offensive titles

```{r,echo=F,eval=T,message=F,warnings=F,fig.cap='Both total and self citations descreased significantly with higher offensive scores, but the effect sizes were small.\\label{fig:offensive}',echo=F}

titles_long <- read.csv(file='Data/Cleaned_Title_Database_with_scorers.csv',header=T)

offensive <- titles_long %>%
  filter(metric == 'offensive') %>%
  mutate(score = as.numeric(score)) %>%
  group_by(Index) %>%
  mutate(avg_offensive = mean(score,na.rm=T)) %>%
  distinct(Index, .keep_all=TRUE)

par(mfrow=c(1,2))
plot(offensive$Total_citations~offensive$avg_offensive,pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Total citations',xlab='Average offensive score',cex.lab=1.5)
offensive_total_glm <- glm(formula = Total_citations~avg_offensive,data = offensive,family = 'poisson')
newdata <- expand.grid(avg_offensive=seq(min(offensive$avg_offensive), max(offensive$avg_offensive), length.out=200))
newdata <- cbind(newdata, citations=predict(offensive_total_glm, newdata, type='response'))
points(newdata,type='l',lwd=4)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(offensive_total_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(offensive_total_glm ))[2,4],2)) ), side = 3)


plot(offensive$Self_citations~offensive$avg_offensive,pch=16,col=rgb(0.5,0.5,0.5,0.3),las=1,ylab='Self citations',xlab='Average offensive score',cex.lab=1.5)
offensive_self_glm <- glm(formula = Self_citations~avg_offensive,data = offensive,family = 'poisson')
newdata <- expand.grid(avg_offensive=seq(min(offensive$avg_offensive), max(offensive$avg_offensive), length.out=200))
newdata <- cbind(newdata, citations=predict(offensive_self_glm, newdata, type='response'))
points(newdata,type='l',lwd=4)
mtext(text= bquote(exp(beta) ==.(round(exp(as.numeric(coef(offensive_self_glm )[2])),2)) ~ "," ~ p == .(round(coef(summary(offensive_self_glm ))[2,4],2)) ), side = 3)



```


## Papers with the highest rated offensive titles (average offensive scores >0.4) include:

- `r paste(offensive$Title[offensive$avg_offensive>0.10][1])`
- `r paste(offensive$Title[offensive$avg_offensive>0.10][2])`
- `r paste(offensive$Title[offensive$avg_offensive>0.10][3])`
- `r paste(offensive$Title[offensive$avg_offensive>0.10][4])`
- `r paste(offensive$Title[offensive$avg_offensive>0.10][5])`
- `r paste(offensive$Title[offensive$avg_offensive>0.10][6])`
- `r paste(offensive$Title[offensive$avg_offensive>0.10][7])`
- `r paste(offensive$Title[offensive$avg_offensive>0.10][8])`
- `r paste(offensive$Title[offensive$avg_offensive>0.10][9])`








```{r,echo=F,eval=F}


pwc <- titles_minimal %>%  pairwise_t_test(Total_citations ~ Taxonomic_name, p.adjust.method = "bonferroni")
ggboxplot(titles_minimal, x = 'Taxonomic_name', y = 'Total_citations',  ylab = "Total citations", xlab = "Taxonomic name")  +  labs(    subtitle = get_test_label(pwc, detailed = TRUE))#,    caption = get_pwc_label(pwc)  )

pwc <- titles_minimal %>%  pairwise_t_test(Total_citations ~ Colon, p.adjust.method = "bonferroni")
ggboxplot(titles_minimal, x = 'Colon', y = 'Total_citations',  ylab = "Total citations", xlab = "Colon")  +  labs(    subtitle = get_test_label(pwc, detailed = TRUE))

ggboxplot(titles_minimal, x = "Taxonomic_name", y = "Total_citations",  ylab = "Total citations", xlab = "Taxonomic name")+ stat_compare_means(method = "t.test",  label.x = 1.5, label.y = 1800)+ stat_compare_means( aes(label = ..p.signif..), label.x = 1.5, label.y = 2000)
```








#


\clearpage

# Random notes that have been deleted