metametazipf_notebook.Rmd

---
title: "The MetaMetaZipf Notebook"
output: html_notebook
---

This is an R Notebook which reproduces the analysis, tables and figures of the working paper "MetaMetaZipf. What do analyses of city size distributions have in common?", as of August 2020. The data is stored on the [MetaZipf](https://github.com/ClementineCttn/MetaZipf) GitHub repository.

```{r setup, include = FALSE}
library(rmarkdown)
library(igraph)
library(scales)
library(wesanderson)
library(ggplot2)
library(ggrepel)
library(ggalt)
library(gridExtra)
library(data.table)
library(stringr)
library(tm)
library(cluster)
library(SnowballC)
library(stargazer)
library(reshape2)


'%!in%' <- function(x,y)!('%in%'(x,y))
norm_vec <- function(x) sqrt(sumNum(x^2))

sumNum = function(x) { 
  s = sum(x, na.rm = T)
  return(s)
}
meanNum = function(x) { 
  s = mean(x, na.rm = T)
  return(s)
}

sdNum = function(x) { 
  s = sd(x, na.rm = T)
  return(s)
}

toSpace <- content_transformer(function(x, pattern) {
  return (gsub(pattern, " ", x))
})

stat.comp<-  function( x,y){
  K <-length(unique(y))
  n <-length(x)
  m <-mean(x)
  TSS <-sum((x-m)^2)
  nk<-table(y)
  mk<-tapply(x,y,mean)
  BSS <-sum(nk* (mk-m)^2)
  result<-c(mk,100.0*BSS/TSS)
  names(result) <-c( paste("G",1:K),"% epl.")
  return(result)
}

```


# Importing Data from the metaZipf project: 
* the information about articles in the corpus (refs)
* the empirical estimations they report along with their specifications (data)
* the citation matrix (cites)
* the cleaned full-texts (docs)
* the lookup table from journal to discipline (J2D)

```{r import}
refs <- read.csv2("data/zipf_refs.csv", sep=';')
cites <- read.csv2("data/zipf_cites.csv", sep=',')
data <- read.csv2("data/zipf_meta.csv", sep=",", dec=".")
cname <- "data/FullText/"
docs <- Corpus(DirSource(cname))   
J2D <- read.csv2("data/journals2Disciplines.csv", sep=';')
```


## Cleaning citation data
```{r cleanCitation, warning=FALSE}
cites_out <- cites[67:1221,]
cites_out$YEAR <- as.numeric(gsub("\\D+", "", cites_out$REFID))
cites_out$JOURNAL <- gsub(".*\\d.", "", cites_out$REFID)
cites_out$AUTHOR <- gsub("\\d.*", "", cites_out$REFID)
cites_out$JOURNAL <- gsub("[[:punct:]]", " ", cites_out$JOURNAL)
cites_out$AUTHOR <- gsub("[[:punct:]]", " ", cites_out$AUTHOR)
cites_out <- data.frame(cites_out,J2D[match(cites_out$JOURNAL, J2D$JOURNAL),])
cites_out$DISCIPLINE <- cites_out$DISCPLINE
cites_out$DISCPLINE <- NULL
cites_out$JOURNAL.1 <- NULL
cites <- rbind(cites[1:66,],cites_out)
cites[is.na(cites)] <- 0
cites_within <- cites[1:66,c(1,6:71)]
head(cites)
```

## Corpus analysis

```{r corpusAnalysis}

j <- cites[1:66,1:5]
j$n <- 1
journals <- aggregate(j[,"n"], by=list(j[,"JOURNAL"]), FUN = sum)
disciplines <- aggregate(j[,"n"], by=list(j[,"DISCIPLINE"]), FUN = sum)
q <- ggplot(journals, aes(x= Group.1, y = x))
 q + geom_lollipop(aes(reorder( Group.1, -x)),color = "orangered", cex=1) +
  coord_flip() +
  labs(x="Journal", y="Number of articles in the corpus")

rownames(cites_within) <- cites_within$REFID
cites_within$REFID <- NULL
cmat <- t(as.matrix(cites_within))
c <- colSums(cmat)
co <- c[order(-c)]
ref <- names(co)
cope <- data.frame(ref, co)
cope_full_names <- data.frame(cope, cites[match(cope$ref, cites$REFID),c("AUTHOR", "YEAR", "JOURNAL")])
cope_full_names$ref <- paste(cope_full_names$AUTHOR, cope_full_names$YEAR, sep=" ")

par(mar = c(2,2,2,2))
q <- ggplot(cope_full_names, aes(x=ref, y = co))
q + geom_lollipop(aes(reorder(ref, -co)),color = "seagreen3", cex=1) +
  coord_flip() +
  labs(x="Reference", y="Number of in-citations from the corpus")


age_journals <- data.frame(aggregate(j[,"YEAR"], by=list(j[,"JOURNAL"]), FUN = mean), journals$x)
age_disciplines <- data.frame(aggregate(j[,"YEAR"], by=list(j[,"DISCIPLINE"]), FUN = mean), disciplines$x)
rownames(age_disciplines) = age_disciplines$Group.1
rownames(age_journals) = age_journals$Group.1
age_disciplines$Group.1 <- NULL
age_journals$Group.1 <- NULL
colnames(age_disciplines) <- c("mean age", "n articles")
colnames(age_journals) <- c("mean age", "n articles")
round(mean(j$YEAR), digit=0)
round(age_disciplines, digit=0)
corpussummary <- data.frame(j, cope[match(j$REFID, cope$ref),])
colnames(corpussummary)[8] <- "in_citations"
mean(corpussummary[corpussummary$in_citations > 10, "YEAR"])
```

## External citation analysis

```{r citationAnalysis}
cites_out <- cites[67:1221,]
cites_out$n_cites <- rowSums(cites_out[,6:71])
single_outcites <- cites_out[order(-cites_out$n_cites),c(1:5,72)] 
s_outcites <- single_outcites[,c("REFID", "n_cites")]
colnames(s_outcites) <- c("ref", "n")
s_outcites <- subset(s_outcites[-100,], n>=5)
q <- ggplot(s_outcites, aes(x=ref, y = n))
q + geom_lollipop(aes(reorder(ref, -n)),color = "coral3", cex=1) +
  coord_flip() +
  labs(x="Reference", y="Number of citations from the corpus (>=5)")
```

Who cited Nitsch (2005)?
```{r exampleCitation}
who <- cites_out[cites_out$REFID == "Nitsch_2005_Journal_Urban_Economics",6:71]
who[1,]
```

Who did not cite G. K. Zipf himself?
```{r noZipfCitation}
nozipf<- as.data.frame(colSums(cites_out[cites_out$REFID %in% c("Zipf_1941_Unity_disunity",
                                                                "Zipf_1949_Human_Behavior_Principle_Least_Effort"),6:71]))
colnames(nozipf) <- "ref_zipf"
nozipf$ref_zipf <- ifelse(nozipf$ref_zipf == 0 , "no ref to zipf", "ref to zipf")
nozipf$REFID <- rownames(nozipf)
cites_for_ref_zipfs <- cites[1:66,]
referencing_zipf <- data.frame(cites_for_ref_zipfs, nozipf[match(cites_for_ref_zipfs$REFID, nozipf$REFID),])
rfczpf <- referencing_zipf[,c("YEAR", "DISCIPLINE", "ref_zipf")]

q <- ggplot(rfczpf, aes(x=YEAR))
q + geom_histogram(aes(y = stat(density), color = ref_zipf, fill = ref_zipf), 
                   alpha = 0.4, position = "identity", binwidth = 1) +
  geom_density(aes(color = ref_zipf), size = 1) +
  scale_color_manual(values = c("Coral2", "dodgerblue3")) +
scale_fill_manual(values = c("Coral2", "dodgerblue3"))  + labs(color = "") + guides(fill = FALSE, size = FALSE)+ theme(legend.position = "top")

table(rfczpf$ref_zipf, rfczpf$DISCIPLINE)
```

Which are the authors cited from Urban Studies?
```{r exCitationFromJournal}
cites_out[cites_out$JOURNAL == "Urban Studies","AUTHOR"]
```

Size of bibliography of corpus articles
```{r sizeBiblio}
n_citations <- as.data.frame(colSums(cites_out[,6:71], na.rm=T))
colnames(n_citations) <- "n_citations"
n_citations$REFID <- rownames(n_citations)
n_citations_corpus <- data.frame(cites_for_ref_zipfs, n_citations[match(cites_for_ref_zipfs$REFID, n_citations$REFID),])
n_citations_corpus$ref <- paste(n_citations_corpus$AUTHOR, n_citations_corpus$YEAR, sep=" ")

q <- ggplot(n_citations_corpus, aes(x=ref, y = n_citations))
q + geom_lollipop(aes(reorder(ref, -n_citations)),color = "dodgerblue3", cex=1) +
  coord_flip() +
  labs(x="Reference", y="Size of bibliography by corpus article")

```

Most cited external aticles
```{r MostCitedArticles}

jou_out <- as.data.frame(table(cites_out$JOURNAL))
jou_out <- jou_out[order(-jou_out$Freq),] 
colnames(jou_out) <- c("ref", "n")
jou <- subset(jou_out[-100,], n>=5)
q <- ggplot(jou, aes(x=ref, y = n))
q + geom_lollipop(aes(reorder(ref, -n)),color = "goldenrod3", cex=1) +
  coord_flip() +
  labs(x="Journal", y="Number of citations from the corpus (>=5)")
```

# Construction of the citation similarity network
```{r citationSimilarityNetwork}

outcitemat <- as.matrix(cites_out[,6:71])
toutcitemat <- t(outcitemat)
refSim <- rownames(toutcitemat)
cosSim <- data.frame()
k=0
ilist <- c()
for(i in refSim){
  ilist <- c(ilist, i)
  for(j in refSim){
    if (j %!in% ilist){
    k <- k + 1
    cosSim[k,1] <- i
    cosSim[k,2] <- j
    vi <- toutcitemat[i,]
    vj <- toutcitemat[j,]
    cosSim[k,3] <- (sumNum(vi * vj)) / (norm_vec(vi) *  norm_vec(vj)) 
    }
  }
}
colnames(cosSim) <- c('i', 'j', 'cosSim')
write.csv(cosSim[order(-cosSim$cosSim),], "simNets/SimilarCitations.csv")

summary(cosSim$cosSim)
cs.cit <- cosSim[cosSim$cosSim >= 0.25,]
g.cit <- graph_from_data_frame(cs.cit, directed=F)

citingNs <- as.data.frame(apply(toutcitemat[rownames(toutcitemat) %in% V(g.cit)$name,],1, FUN = norm_vec))
colnames(citingNs) <- "citingN"
citingNs$ref <- rownames(citingNs)
orderedName <- data.frame(V(g.cit)$name)
citingN <- data.frame(orderedName, citingNs[match(orderedName$V.g.cit..name, citingNs$ref),])[,"citingN"]
citingN <- citingN[!is.na(citingN)] 

par(mar = c(0,0,0,0))
clln.cit <- cluster_louvain(g.cit)
layout <- layout_nicely(g.cit,2)
g.cit$layout <- layout
plot(g.cit, edge.width = sqrt(cs.cit$cosSim) * 2, vertex.size = citingN,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = membership(clln.cit))
```

# Construction of the discipline similarity network
```{r disciplineNetwork}
discipline2ref <- aggregate(cites_out[,c(6:71)], by=list(cites_out$DISCIPLINE), FUN = sumNum)
rownames(discipline2ref) <- discipline2ref$Group.1
discipline2ref$Group.1 <- NULL
cols <- colnames(discipline2ref)
for(i in cols){
  absolute <- discipline2ref[,i]
  total <- sum(absolute)
  relative <- absolute/total
  discipline2ref[,paste0("freq_",i)] <- relative
}

disc2ref <- as.matrix(discipline2ref[,67:132])
colnames(disc2ref) <- cols
scaled_disc2ref <- disc2ref / colSums(disc2ref)
tdisc2ref <- t(disc2ref)

refSim <- rownames(tdisc2ref)
cosSim <- data.frame()
k=0
ilist <- c()
for(i in refSim){
  ilist <- c(ilist, i)
  for(j in refSim){
    if (j %!in% ilist){
      k <- k + 1
      cosSim[k,1] <- i
      cosSim[k,2] <- j
      vi <- tdisc2ref[i,]
      vj <- tdisc2ref[j,]
      cosSim[k,3] <- (sumNum(vi * vj)) / (norm_vec(vi) *  norm_vec(vj)) 
    }
  }
}
colnames(cosSim) <- c('i', 'j', 'cosSim')
write.csv(cosSim[order(-cosSim$cosSim),], "simNets/SimilarDisciplinesCited_rel.csv")
summary(cosSim$cosSim)
cs.cit <- cosSim[cosSim$cosSim >= 0.9,]
g.cit <- graph_from_data_frame(cs.cit, directed=F)

disCitedNs <- as.data.frame(apply(toutcitemat[rownames(tdisc2ref) %in% V(g.cit)$name,],1, FUN = norm_vec))
colnames(disCitedNs) <- "disCitedNs"
disCitedNs$ref <- rownames(disCitedNs)
orderedName <- data.frame(V(g.cit)$name)
disCitedNs <- data.frame(orderedName, disCitedNs[match(orderedName$V.g.cit..name, disCitedNs$ref),])[,"disCitedNs"]
disCitedNs <- disCitedNs[!is.na(disCitedNs)] 

colnames(orderedName) <- "REFID"
refinfo <- cites[1:66,c("REFID", "DISCIPLINE")]
own_discipline <- data.frame(orderedName, refinfo[match(orderedName$REFID, refinfo$REFID),])
odisc <- own_discipline$DISCIPLINE

layout <- layout_nicely(g.cit,2)
g.cit$layout <- layout
par(mar = c(0,0,0,0))
plot(g.cit, edge.width = sqrt(cs.cit$cosSim) * 2, vertex.size = disCitedNs,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = odisc)


```

# Construction of full-text similarity network
```{r fullTextSim}
docs <- tm_map(docs, toSpace, "-")
docs <- tm_map(docs, toSpace, ":")
docs <- tm_map(docs, toSpace, "'")
docs <- tm_map(docs, toSpace, "`")
docs <- tm_map(docs, toSpace, "‘")
docs <- tm_map(docs, toSpace, "_")
docs <- tm_map(docs, toSpace, "–")
docs <- tm_map(docs, toSpace, "/")
docs <- tm_map(docs, removePunctuation)
docs <- tm_map(docs, tolower)   
docs <- tm_map(docs, removeNumbers)   
docs <- tm_map(docs, removeWords, stopwords("english"))   
dtm <- DocumentTermMatrix(docs)   
dtm.m = as.matrix(dtm)
tdm <- TermDocumentMatrix(docs)   
tdm.m = as.matrix(tdm)

 tdm.df = as.data.frame(tdm.m)
 tdm.df.f <- tdm.df / colSums(tdm.df)

n_words <- as.data.frame(colSums(tdm.df, na.rm=T))
 colnames(n_words) <- "n_words"
 n_words$REFID <- substr(rownames(n_words), 1,8)
 n_words_corpus <- data.frame(cites_for_ref_zipfs, n_words[match(cites_for_ref_zipfs$REFID, n_words$REFID),])
 n_words_corpus$ref <- paste(n_words_corpus$AUTHOR, n_words_corpus$YEAR, sep=" ")
 n_words_corpus <- n_words_corpus[n_words_corpus$AUTHOR != "Krakover S.",]
 q <- ggplot(n_words_corpus, aes(x=ref, y = n_words))
 q + geom_lollipop(aes(reorder(ref, -n_words)),color = "seagreen4", cex=1) +
   coord_flip() +
   labs(x="Reference", y="Number of words by corpus article")
 
 termmat <- t(as.matrix(tdm.df.f))
 
 refSimTerm <- rownames(termmat)
 cosSimTerm <- data.frame()
 k=0
 ilist <- c()
 for(i in refSimTerm){
   ilist <- c(ilist, i)
   for(j in refSimTerm){
     if (j %!in% ilist){
       k <- k + 1
       cosSimTerm[k,1] <- i
       cosSimTerm[k,2] <- j
       vi <- termmat[i,]
       vj <- termmat[j,]
       cosSimTerm[k,3] <- (sumNum(vi * vj)) / (norm_vec(vi) *  norm_vec(vj)) 
     }
   }
 }
 colnames(cosSimTerm) <- c('i', 'j', 'cosSimTerm')
 write.csv(cosSimTerm[order(-cosSimTerm$cosSimTerm),], "simNets/SimilarWording.csv")
 summary(cosSimTerm)

 cs <- cosSimTerm[cosSimTerm$cosSimTerm > 0.7,]
 g.term <- graph_from_data_frame(cs, directed=F)
 totalTerm <- as.data.frame(apply(termmat[rownames(termmat) %in% V(g.term)$name,],1, FUN = norm_vec))
 colnames(totalTerm) <- "totalTerms"
 totalTerm$ref <- rownames(totalTerm)
 orderedName <- data.frame(V(g.term)$name)
 totalTerms <- data.frame(orderedName, totalTerm[match(orderedName$V.g.term..name, totalTerm$ref),])[,"totalTerms"]
 totalTerms <- totalTerms[!is.na(totalTerms)] 
 clln.term <- cluster_louvain(g.term)
 layout <- layout_nicely(g.term,2)
 g.term$layout <- layout
 par(mar = c(0,0,0,0))
 plot(g.term, edge.width = sqrt(cs$cosSimTerm) * 2, vertex.size = totalTerms *50,
      vertex.label.cex = 0.7,  
      edge.curved=.2, vertex.color = membership(clln.term))
 
```

# Construction country similarity network
```{r countrySim}

  paperList <- colnames(cites_out[,6:71])
  countryData <- data[data$TERRITORY_TYPE == "Country" & data$REFID %in% paperList,]
country2Ref <- table(countryData$CNTR_ID, countryData$REFID)[-1,]
country2Ref.mat <- as.matrix(country2Ref[rowSums(country2Ref)>0,colnames(country2Ref) %in% paperList])
country2Ref.mat[country2Ref.mat>0] <- 1
countrymat <- t(country2Ref.mat)

refSim <- rownames(countrymat)
cosSimC <- data.frame()
k=0
ilist <- c()
for(i in refSim){
  ilist <- c(ilist, i)
  for(j in refSim){
    if (j %!in% ilist){
      k <- k + 1
      cosSimC[k,1] <- i
      cosSimC[k,2] <- j
      vi <- countrymat[i,]
      vj <- countrymat[j,]
      s <- (sumNum(vi * vj)) / (norm_vec(vi) *  norm_vec(vj)) 
      cosSimC[k,3] <- ifelse(!is.na(s), s, 0)
    }
  }
}
colnames(cosSimC) <- c('i', 'j', 'cosSimC')
write.csv(cosSimC[order(-cosSimC$cosSimC),], "simNets/SimilarCountries.csv")
summary(cosSimC)
cs.cntr <- cosSimC[cosSimC$cosSimC >= 0.2,]
g.cntr <- graph_from_data_frame(cs.cntr, directed=F)
countryN <- as.data.frame(apply(countrymat[rownames(countrymat) %in% V(g.cntr)$name,],1, FUN = norm_vec))
colnames(countryN) <- "n_countries"
countryN$ref <- rownames(countryN)
orderedName <- data.frame(V(g.cntr)$name)
orderedCountryN <- data.frame(orderedName, countryN[match(orderedName$V.g.cntr..name, countryN$ref),])[,"n_countries"]
orderedCountryN <- orderedCountryN[!is.na(orderedCountryN)] 
clln.cntr <- cluster_louvain(g.cntr)
layout <- layout_nicely(g.cntr,2)
g.cntr$layout <- layout
 par(mar = c(0,0,0,0))
 plot(g.cntr, edge.width = sqrt(cs.cntr$cosSimC) * 2, vertex.size = orderedCountryN * 2,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = membership(clln.cntr))

```

# Construction city definition similarity network
```{r cityDefSim}

cityData <- data[!is.na(data$URBANSCALE) & data$REFID %in% paperList,]
city2Ref <- table(cityData$URBANSCALE, cityData$REFID)[,]
city2Ref.mat <- as.matrix(city2Ref[rowSums(city2Ref)>0,colnames(city2Ref) %in% paperList])
city2Ref.mat[city2Ref.mat>0] <- 1
citymat <- t(city2Ref.mat)

refSim <- rownames(citymat)
cosSimCt <- data.frame()
k=0
ilist <- c()
for(i in refSim){
  ilist <- c(ilist, i)
  for(j in refSim){
    if (j %!in% ilist){
      k <- k + 1
      cosSimCt[k,1] <- i
      cosSimCt[k,2] <- j
      vi <- citymat[i,]
      vj <- citymat[j,]
      s <- (sumNum(vi * vj)) / (norm_vec(vi) *  norm_vec(vj)) 
      cosSimCt[k,3] <- ifelse(!is.na(s), s, 0)
    }
  }
}
colnames(cosSimCt) <- c('i', 'j', 'cosSimCt')
write.csv(cosSimCt[order(-cosSimCt$cosSimCt),], "simNets/SimilarCities.csv")
summary(cosSimCt)
cs.city <- cosSimCt[cosSimCt$cosSimCt >= 0.1,]

g.city <- graph_from_data_frame(cs.city, directed=F)
cityN <- as.data.frame(apply(citymat[rownames(citymat) %in% V(g.city)$name,],1, FUN = norm_vec))
colnames(cityN) <- "n_cities"
cityN$ref <- rownames(cityN)
orderedName <- data.frame(V(g.city)$name)
orderedcityN <- data.frame(orderedName, cityN[match(orderedName$V.g.city..name, cityN$ref),])[,"n_cities"]
orderedcityN <- orderedcityN[!is.na(orderedcityN)] 

clln.city <- cluster_louvain(g.city)
layout <- layout_nicely(g.city,2)
g.city$layout <- layout
 par(mar = c(0,0,0,0))
 plot(g.city, edge.width = sqrt(cs.city$cosSimCt) * 2, vertex.size = orderedcityN * 2,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = membership(clln.city))

```


# Construction decade similarity network
```{r decadeSim}

decadeData <- data[!is.na(data$DATE) & data$REFID %in% paperList,]
decadeData$DECADE <- paste0(substr(as.character(decadeData$DATE), 1, 3), 0, "s")
decade2Ref <- table(decadeData$DECADE, decadeData$REFID)[,]
decade2Ref.mat <- as.matrix(decade2Ref[rowSums(decade2Ref)>0,colnames(decade2Ref) %in% paperList])
decade2Ref.mat[decade2Ref.mat>0] <- 1
decademat <- t(decade2Ref.mat)

refSim <- rownames(decademat)
cosSimCt <- data.frame()
k=0
ilist <- c()
for(i in refSim){
  ilist <- c(ilist, i)
  for(j in refSim){
    if (j %!in% ilist){
      k <- k + 1
      cosSimCt[k,1] <- i
      cosSimCt[k,2] <- j
      vi <- decademat[i,]
      vj <- decademat[j,]
      s <- (sumNum(vi * vj)) / (norm_vec(vi) *  norm_vec(vj))
      cosSimCt[k,3] <- ifelse(!is.na(s), s, 0)
    }
  }
}
colnames(cosSimCt) <- c('i', 'j', 'cosSimCt')
write.csv(cosSimCt[order(-cosSimCt$cosSimCt),], "simNets/SimilarDecades.csv")
summary(cosSimCt$cosSimCt)
cs.decade <- cosSimCt[cosSimCt$cosSimCt >= 0.65,]

g.decade <- graph_from_data_frame(cs.decade, directed=F)
decadeN <- as.data.frame(apply(decademat[rownames(decademat) %in% V(g.decade)$name,],1, FUN = norm_vec))
colnames(decadeN) <- "n_decades"
decadeN$ref <- rownames(decadeN)
orderedName <- data.frame(V(g.decade)$name)
ordereddecadeN <- data.frame(orderedName, decadeN[match(orderedName$V.g.decade..name, decadeN$ref),])[,"n_decades"]
ordereddecadeN <- ordereddecadeN[!is.na(ordereddecadeN)]

clln.decade <- cluster_louvain(g.decade)
layout <- layout_nicely(g.decade,2)
g.decade$layout <- layout
par(mar = c(0,0,0,0))
plot(g.decade, edge.width = sqrt(cs.decade$cosSimCt) * 2, vertex.size = ordereddecadeN * 2,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = membership(clln.decade))
 
```

# Construction of the networks of alpha similarity
```{r alphaSim}
alphaData <- data[!is.na(data$ALPHALOTKA) & data$REFID %in% paperList,]
alphaData$n <- 1
meanAlphaPerRef <- aggregate(alphaData[,"ALPHALOTKA"], by=list(alphaData$REFID), FUN=meanNum)
sdAlphaPerRef <- aggregate(alphaData[,"ALPHALOTKA"], by=list(alphaData$REFID), FUN=sdNum)
sdAlphaPerRef[is.na(sdAlphaPerRef$x),"x"] <- 0
nAlphaPerRef <- aggregate(alphaData[,"n"], by=list(alphaData$REFID), FUN=sumNum)

similarAlpha <- cbind(meanAlphaPerRef, sdAlphaPerRef, nAlphaPerRef)
rownames(similarAlpha) <- similarAlpha$Group.1
similarAlpha[,c(1,3, 5)] <- NULL
colnames(similarAlpha) <- c("meanAlpha", "sdAlpha", "nAlpha")

alphamat <- as.matrix(similarAlpha)

meanData <- mean(alphaData[,"ALPHALOTKA"])
sdData <- sd(alphaData[,"ALPHALOTKA"])
meanN <- mean(similarAlpha$nAlpha)

refSim <- rownames(alphamat)
diff_mean <- data.frame()
diff_sd <- data.frame()
diff_n <- data.frame()
k=0
ilist <- c()
for(i in refSim){
  ilist <- c(ilist, i)
  for(j in refSim){
    if (j %!in% ilist){
      k <- k + 1
      diff_mean[k,1] <- i
      diff_mean[k,2] <- j
      mi <- similarAlpha[i,1]
      mj <- similarAlpha[j,1]
      m <- abs((mi - mj) / meanData)
      diff_mean[k,3] <- m
      
      diff_sd[k,1] <- i
      diff_sd[k,2] <- j
      sdi <-similarAlpha[i,2]
      sdj <- similarAlpha[j,2]
      sd <- abs((sdi - sdj) / sdData)
      diff_sd[k,3] <- sd
      
      diff_n[k,1] <- i
      diff_n[k,2] <- j
      ni <-similarAlpha[i,3]
      nj <- similarAlpha[j,3]
      n <- abs((ni - nj) / meanN)
      diff_n[k,3] <- n
    }
  }
}

colnames(diff_mean) <- c('i', 'j', 'diff_mean')
write.csv(diff_mean[order(diff_mean$diff_mean),], "simNets/diff_mean.csv")
summary(diff_mean$diff_mean)
colnames(diff_sd) <- c('i', 'j', 'diff_sd')
write.csv(diff_sd[order(diff_sd$diff_sd),], "simNets/diff_sd.csv")
summary(diff_sd$diff_sd)
colnames(diff_n) <- c('i', 'j', 'diff_n')
write.csv(diff_n[order(diff_n$diff_n),], "simNets/diff_n.csv")
summary(diff_n$diff_n)

```

## Similarity of mean alpha
```{r meanAlphaSim}
cs.alpha <- diff_mean[diff_mean$diff_mean < 0.025,]
g.alpha <- graph_from_data_frame(cs.alpha, directed=F)
alphaM <- as.data.frame(apply(alphamat[rownames(alphamat) %in% V(g.alpha)$name,],1, FUN = norm_vec))
colnames(alphaM) <- "mean_alphas"
alphaM$ref <- rownames(alphaM)
orderedName <- data.frame(V(g.alpha)$name)
orderedalphaM <- data.frame(orderedName, alphaM[match(orderedName$V.g.alpha..name, alphaM$ref),])[,"mean_alphas"]
orderedalphaM <- orderedalphaM[!is.na(orderedalphaM)]
similarAlpha$REFID <- rownames(similarAlpha)
clln.alpha <- cluster_louvain(g.alpha)
layout <- layout_nicely(g.alpha,2)
g.alpha$layout <- layout

memb <- as.list(membership(clln.alpha))

vertexData <- data.frame(orderedREFID= names(membership(clln.alpha)))
vertexData <- data.frame(vertexData, similarAlpha[match(vertexData$orderedREFID, similarAlpha$REFID),])

for (r in vertexData$REFID){
  vertexData[vertexData$REFID == r,"memb_mean"] <- memb[r]
}

par(mar = c(0,0,0,0))
plot(g.alpha, edge.width = sqrt(cs.alpha$diff_mean) * 2, vertex.size = vertexData$meanAlpha^2 * 5,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = membership(clln.alpha))

# mean average value of alpha reported by group
aggregate(vertexData[,"meanAlpha"], by=list(vertexData$memb_mean), FUN=meanNum)

```

## Similarity of sd alpha
```{r sdAlphaSim}
cs.alpha <- diff_sd[diff_sd$diff_sd < 0.1,]
g.alpha <- graph_from_data_frame(cs.alpha, directed=F)
alphaSD <- as.data.frame(apply(alphamat[rownames(alphamat) %in% V(g.alpha)$name,],1, FUN = norm_vec))
colnames(alphaSD) <- "mean_sds"
alphaSD$ref <- rownames(alphaSD)
orderedName <- data.frame(V(g.alpha)$name)
orderedalphaSD <- data.frame(orderedName, alphaSD[match(orderedName$V.g.alpha..name, alphaSD$ref),])[,"mean_sds"]
orderedalphaSD <- orderedalphaSD[!is.na(orderedalphaSD)]

clln.alpha <- cluster_louvain(g.alpha)
layout <- layout_nicely(g.alpha,2)
g.alpha$layout <- layout

memb <- as.list(membership(clln.alpha))

vertexData <- data.frame(orderedREFID= names(membership(clln.alpha)))
vertexData <- data.frame(vertexData, similarAlpha[match(vertexData$orderedREFID, similarAlpha$REFID),])

for (r in vertexData$REFID){
  vertexData[vertexData$REFID == r,"memb_sd"] <- memb[r]
}

par(mar = c(0,0,0,0))
plot(g.alpha, edge.width = sqrt(cs.alpha$diff_sd) * 2, vertex.size = vertexData$sdAlpha * 50,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = membership(clln.alpha))

# mean average value of alpha reported by group
aggregate(vertexData[,"sdAlpha"], by=list(vertexData$memb_sd), FUN=meanNum)

```


## Similarity of n alpha
```{r nAlphaSim}
cs.alpha <- diff_n[diff_n$diff_n < 0.1,]
g.alpha <- graph_from_data_frame(cs.alpha, directed=F)
alphaN <- as.data.frame(apply(alphamat[rownames(alphamat) %in% V(g.alpha)$name,],1, FUN = norm_vec))
colnames(alphaN) <- "mean_ns"
alphaN$ref <- rownames(alphaN)
orderedName <- data.frame(V(g.alpha)$name)
orderedalphaN <- data.frame(orderedName, alphaN[match(orderedName$V.g.alpha..name, alphaN$ref),])[,"mean_ns"]
orderedalphaN <- orderedalphaN[!is.na(orderedalphaN)]

clln.alpha <- cluster_louvain(g.alpha)
layout <- layout_nicely(g.alpha,2)
g.alpha$layout <- layout

memb <- as.list(membership(clln.alpha))

vertexData <- data.frame(orderedREFID= names(membership(clln.alpha)))
vertexData <- data.frame(vertexData, similarAlpha[match(vertexData$orderedREFID, similarAlpha$REFID),])

for (r in vertexData$REFID){
  vertexData[vertexData$REFID == r,"memb_n"] <- memb[r]
}

par(mar = c(0,0,0,0))
plot(g.alpha, edge.width = sqrt(cs.alpha$diff_n) * 2, vertex.size = vertexData$nAlpha * 0.4,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = membership(clln.alpha))

# mean average value of alpha reported by group
aggregate(vertexData[,"nAlpha"], by=list(vertexData$memb_n), FUN=meanNum)


```

# Regression models of similarity in mean and sd alpha
## Preparation of data 
```{r prepmodelling}

# import diff in alpha means
meanAlphaDyads <- read.csv2("simNets/diff_mean.csv", sep=",", stringsAsFactors = F)[,-1]
meanAlphaDyads$dyadID <- paste(meanAlphaDyads$i, meanAlphaDyads$j, sep = "_")
meanAlphaDyads$meanAlphaSim <- -as.numeric(meanAlphaDyads$diff_mean)
meanAlphaDyads[,1:3] <- NULL
DYADS <- meanAlphaDyads
DYAD_ID_order <- DYADS$dyadID

# import diff in alpha sd
sdAlphaDyads <- read.csv2("simNets/diff_sd.csv", sep=",", stringsAsFactors = F)[,-1]
sdAlphaDyads$dyadIJ <- paste(substr(sdAlphaDyads$i, 1, 8), 
                             substr(sdAlphaDyads$j, 1, 8), sep = "_")
sdAlphaDyads$dyadJI <- paste(substr(sdAlphaDyads$i, 1, 8), 
                             substr(sdAlphaDyads$j, 1, 8), sep = "_")
sdAlphaDyads$dyadID <- ifelse(sdAlphaDyads$dyadIJ %in% DYAD_ID_order, sdAlphaDyads$dyadIJ, sdAlphaDyads$dyadJI)

sdAlphaDyads$sdAlphaSim <- -as.numeric(sdAlphaDyads$diff_sd)
sdAlphaDyads[,1:5] <- NULL
DYADS <- data.frame(DYADS, sdAlphaDyads[match(DYADS$dyadID, sdAlphaDyads$dyadID),"sdAlphaSim"])

# import diff in n alpha 
nAlphaDyads <- read.csv2("simNets/diff_n.csv", sep=",", stringsAsFactors = F)[,-1]
nAlphaDyads$dyadIJ <- paste(substr(nAlphaDyads$i, 1, 8), 
                             substr(nAlphaDyads$j, 1, 8), sep = "_")
nAlphaDyads$dyadJI <- paste(substr(nAlphaDyads$i, 1, 8), 
                             substr(nAlphaDyads$j, 1, 8), sep = "_")
nAlphaDyads$dyadID <- ifelse(nAlphaDyads$dyadIJ %in% DYAD_ID_order, nAlphaDyads$dyadIJ, nAlphaDyads$dyadJI)

nAlphaDyads$nAlphaSim <- -as.numeric(nAlphaDyads$diff_n)
nAlphaDyads[,1:5] <- NULL
DYADS <- data.frame(DYADS, nAlphaDyads[match(DYADS$dyadID, nAlphaDyads$dyadID),"nAlphaSim"])

# import similarity in wording as explanation
wordingDyads <- read.csv2("simNets/SimilarWording.csv", sep=",", stringsAsFactors = F)[,-1]
head(wordingDyads)
wordingDyads$dyadIJ <- paste(substr(wordingDyads$i, 1, 8), 
                             substr(wordingDyads$j, 1, 8), sep = "_")
wordingDyads$dyadJI <- paste(substr(wordingDyads$i, 1, 8), 
                             substr(wordingDyads$j, 1, 8), sep = "_")
wordingDyads$dyadID <- ifelse(wordingDyads$dyadIJ %in% DYAD_ID_order, wordingDyads$dyadIJ, wordingDyads$dyadJI)

wordingDyads$wordingCosSim <- as.numeric(wordingDyads$cosSimTerm)
wordingDyads[,1:5] <- NULL
DYADS <- data.frame(DYADS, wordingDyads[match(DYADS$dyadID, wordingDyads$dyadID),"wordingCosSim"])

# import similarity in citation as explanation
citationDyads <- read.csv2("simNets/SimilarCitations.csv", sep=",", stringsAsFactors = F)[,-1]
citationDyads$dyadIJ <- paste(citationDyads$i, citationDyads$j, sep = "_")
citationDyads$dyadJI <- paste(citationDyads$j, citationDyads$i, sep = "_")
citationDyads$dyadID <- ifelse(citationDyads$dyadIJ %in% DYAD_ID_order, citationDyads$dyadIJ, citationDyads$dyadJI)
citationDyads$citationCosSim <- as.numeric(citationDyads$cosSim)
citationDyads[,1:5] <- NULL
DYADS <- data.frame(DYADS, citationDyads[match(DYADS$dyadID, citationDyads$dyadID),"citationCosSim"])

# import similarity in discipline as explanation
disciplineDyads <- read.csv2("simNets/SimilarDisciplinesCited_rel.csv", sep=",", stringsAsFactors = F)[,-1]
disciplineDyads$dyadIJ <- paste(disciplineDyads$i, disciplineDyads$j, sep = "_")
disciplineDyads$dyadJI <- paste(disciplineDyads$j, disciplineDyads$i, sep = "_")
disciplineDyads$dyadID <- ifelse(disciplineDyads$dyadIJ %in% DYAD_ID_order, disciplineDyads$dyadIJ, disciplineDyads$dyadJI)
disciplineDyads$disciplineCosSim <- as.numeric(disciplineDyads$cosSim)
disciplineDyads[,1:5] <- NULL
DYADS <- data.frame(DYADS, disciplineDyads[match(DYADS$dyadID, disciplineDyads$dyadID),"disciplineCosSim"])

# import similarity in countries studied as control
countryDyads <- read.csv2("simNets/SimilarCountries.csv", sep=",", stringsAsFactors = F)[,-1]
countryDyads$dyadIJ <- paste(countryDyads$i, countryDyads$j, sep = "_")
countryDyads$dyadJI <- paste(countryDyads$j, countryDyads$i, sep = "_")
countryDyads$dyadID <- ifelse(countryDyads$dyadIJ %in% DYAD_ID_order, countryDyads$dyadIJ, countryDyads$dyadJI)
countryDyads$countryCosSim <- as.numeric(countryDyads$cosSimC)
countryDyads[,1:5] <- NULL
DYADS <- data.frame(DYADS, countryDyads[match(DYADS$dyadID, countryDyads$dyadID),"countryCosSim"])

# import similarity in decades studied as control
decadeDyads <- read.csv2("simNets/SimilarDecades.csv", sep=",", stringsAsFactors = F)[,-1]
decadeDyads$dyadIJ <- paste(decadeDyads$i, decadeDyads$j, sep = "_")
decadeDyads$dyadJI <- paste(decadeDyads$j, decadeDyads$i, sep = "_")
decadeDyads$dyadID <- ifelse(decadeDyads$dyadIJ %in% DYAD_ID_order, decadeDyads$dyadIJ, decadeDyads$dyadJI)
decadeDyads$decadeCosSim <- as.numeric(decadeDyads$cosSimCt)
decadeDyads[,1:5] <- NULL
DYADS <- data.frame(DYADS, decadeDyads[match(DYADS$dyadID, decadeDyads$dyadID),"decadeCosSim"])

# import similarity in city definition used as control
cityDefDyads <- read.csv2("simNets/SimilarCities.csv", sep=",", stringsAsFactors = F)[,-1]
cityDefDyads$dyadIJ <- paste(cityDefDyads$i, cityDefDyads$j, sep = "_")
cityDefDyads$dyadJI <- paste(cityDefDyads$j, cityDefDyads$i, sep = "_")
cityDefDyads$dyadID <- ifelse(cityDefDyads$dyadIJ %in% DYAD_ID_order, cityDefDyads$dyadIJ, cityDefDyads$dyadJI)
cityDefDyads$cityDefCosSim <- as.numeric(cityDefDyads$cosSimCt)
cityDefDyads[,1:5] <- NULL
DYADS <- data.frame(DYADS, cityDefDyads[match(DYADS$dyadID, cityDefDyads$dyadID),"cityDefCosSim"])

colnames(DYADS) <- c("dyadID","meanAlphaSim", "sdAlphaSim", "nAlphaSim", 
                     "fullTextCosSim", "externalCitationCosSim", "disciplineCitationCosSim",
                     "countriesCosSim", "decadesCosSim", "cityDefCosSim")

# scale all variables
DYADS$meanAlphaSim_scaled <- scale(DYADS$meanAlphaSim)
DYADS$sdAlphaSim_scaled <- scale(DYADS$sdAlphaSim)
DYADS$nAlphaSim_scaled <- scale(DYADS$nAlphaSim)
DYADS$fullTextCosSim_scaled <- scale(DYADS$fullTextCosSim)
DYADS$externalCitationCosSim_scaled <- scale(DYADS$externalCitationCosSim)
DYADS$disciplineCitationCosSim_scaled <- scale(DYADS$disciplineCitationCosSim)
DYADS$countriesCosSim_scaled <- scale(DYADS$countriesCosSim)
DYADS$decadesCosSim_scaled <- scale(DYADS$decadesCosSim)
DYADS$cityDefCosSim_scaled <- scale(DYADS$cityDefCosSim)

DYADSwithNA <- DYADS
DYADS <- DYADS[complete.cases(DYADS),]

modeldata<- DYADS[,c("dyadID","meanAlphaSim_scaled", "sdAlphaSim_scaled", "nAlphaSim_scaled",
                     "fullTextCosSim_scaled", "externalCitationCosSim_scaled", "disciplineCitationCosSim_scaled",
                     "countriesCosSim_scaled", "decadesCosSim_scaled", "cityDefCosSim_scaled")]
colnames(modeldata) <- c("ID", "alpha", "sd_alpha", "n_estim","text", "citation", "discipline", "country", "decade", "city")

```


## Models of mean alpha
```{r meanAlphaModels}

ftm_a <- lm(data = modeldata,
            formula = alpha ~ text,
            na.action = na.omit)
ecm_a <- lm(data = modeldata,
            formula = alpha ~ citation,
            na.action = na.omit)
dcm_a <- lm(data = modeldata,
            formula = alpha ~ discipline,
            na.action = na.omit)
ctrm_a <- lm(data = modeldata,
            formula = alpha ~ n_estim +
             country * decade * city,
            na.action = na.omit)
fullmodel_a <- lm(data = modeldata,
            formula = alpha ~ 
              text + citation + discipline + n_estim +
              country * decade * city,
            na.action = na.omit)

sgm<- stargazer(ftm_a, ecm_a, dcm_a,
          ctrm_a,
          fullmodel_a,
          type="text", out="results/model_comp_results_meanAlpha.html", font.size = "small", column.sep.width = "1pt")

```

## Models of sd alpha
```{r sdAlphaModels}

ftm_sd <- lm(data = modeldata,
            formula = sd_alpha ~ text,
            na.action = na.omit)
ecm_sd <- lm(data = modeldata,
            formula = sd_alpha ~ citation,
            na.action = na.omit)
dcm_sd <- lm(data = modeldata,
            formula = sd_alpha ~ discipline,
            na.action = na.omit)
 n_sd <- lm(data = modeldata,
             formula = sd_alpha ~ n_estim,
             na.action = na.omit)
ctrm_sd <- lm(data = modeldata,
             formula = sd_alpha ~ 
               country * decade * city,
             na.action = na.omit)
fullmodel_sd <- lm(data = modeldata,
                  formula = sd_alpha ~ 
                    text + citation + discipline + n_estim +
                    country * decade * city,
                  na.action = na.omit)

sgm<- stargazer(ftm_sd, ecm_sd, dcm_sd,n_sd,
                ctrm_sd,
                fullmodel_sd,
                type="text", out="results/model_comp_results_sdAlpha.html", font.size = "small", column.sep.width = "1pt")

```

## Residuals analysis
```{r residuals, warning=FALSE}

modeldata$res <- fullmodel_a$residuals
head(modeldata[order(-modeldata$res),])
modeldata$i <- substr(modeldata$ID, 1,8)
modeldata$j <- substr(modeldata$ID, 10, 17)

cs.residuals_resPos <- modeldata[modeldata$res > 1.1,c("i", "j", "res")]
g.residuals <- graph_from_data_frame(cs.residuals_resPos, directed=F)
resPos <- as.data.frame(apply(alphamat[rownames(alphamat) %in% V(g.residuals)$name,],1, FUN = norm_vec))
colnames(resPos) <- "residuals"
resPos$ref <- rownames(resPos)
orderedName <- data.frame(V(g.residuals)$name)
orderedresidualsRes <- data.frame(orderedName, resPos[match(orderedName$V.g.residuals..name, resPos$ref),])[,"residuals"]
orderedresidualsRes <- orderedresidualsRes[!is.na(orderedresidualsRes)]

layout <- layout_nicely(g.residuals,2)
g.residuals$layout <- layout
par(mar = c(0,0,1,0))
plot(g.residuals, edge.width = sqrt(cs.residuals_resPos$res) * 2, vertex.size = 2,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = "green", edge.color = "coral3", 
     main = "Most positive residuals in mean alpha similarity")



cs.residuals_resNeg <- modeldata[modeldata$res < -2,c("i", "j", "res")]
g.residuals <- graph_from_data_frame(cs.residuals_resNeg, directed=F)
resNeg <- as.data.frame(apply(alphamat[rownames(alphamat) %in% V(g.residuals)$name,],1, FUN = norm_vec))
colnames(resNeg) <- "residuals"
resNeg$ref <- rownames(resNeg)
orderedName <- data.frame(V(g.residuals)$name)
orderedresidualsRes <- data.frame(orderedName, resNeg[match(orderedName$V.g.residuals..name, resNeg$ref),])[,"residuals"]
orderedresidualsRes <- orderedresidualsRes[!is.na(orderedresidualsRes)]

layout <- layout_nicely(g.residuals,2)
g.residuals$layout <- layout
par(mar = c(0,0,1,0))
plot(g.residuals, edge.width = sqrt(cs.residuals_resNeg$res) * 2, vertex.size = 2,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = "orange", edge.color = "dodgerblue3", 
     main = "Most negative residuals in mean alpha similarity")



modeldata$res <- fullmodel_sd$residuals
modeldata$i <- substr(modeldata$ID, 1,8)
modeldata$j <- substr(modeldata$ID, 10, 17)

cs.residuals_resPos <- modeldata[modeldata$res > 1,c("i", "j", "res")]
g.residuals <- graph_from_data_frame(cs.residuals_resPos, directed=F)
resPos <- as.data.frame(apply(alphamat[rownames(alphamat) %in% V(g.residuals)$name,],1, FUN = norm_vec))
colnames(resPos) <- "residuals"
resPos$ref <- rownames(resPos)
orderedName <- data.frame(V(g.residuals)$name)
orderedresidualsRes <- data.frame(orderedName, resPos[match(orderedName$V.g.residuals..name, resPos$ref),])[,"residuals"]
orderedresidualsRes <- orderedresidualsRes[!is.na(orderedresidualsRes)]

#clln.residuals <- cluster_louvain(g.residuals)
layout <- layout_nicely(g.residuals,2)
g.residuals$layout <- layout
par(mar = c(0,0,1,0))
plot(g.residuals, edge.width = sqrt(cs.residuals_resPos$res) * 2, vertex.size = 2,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = "green", edge.color = "coral3", 
     main = "Most positive residuals in sd alpha similarity")



cs.residuals_resNeg <- modeldata[modeldata$res < -1.5,c("i", "j", "res")]
g.residuals <- graph_from_data_frame(cs.residuals_resNeg, directed=F)
resNeg <- as.data.frame(apply(alphamat[rownames(alphamat) %in% V(g.residuals)$name,],1, FUN = norm_vec))
colnames(resNeg) <- "residuals"
resNeg$ref <- rownames(resNeg)
orderedName <- data.frame(V(g.residuals)$name)
orderedresidualsRes <- data.frame(orderedName, resNeg[match(orderedName$V.g.residuals..name, resNeg$ref),])[,"residuals"]
orderedresidualsRes <- orderedresidualsRes[!is.na(orderedresidualsRes)]

layout <- layout_nicely(g.residuals,2)
g.residuals$layout <- layout
par(mar = c(0,0,1,0))
plot(g.residuals, edge.width = sqrt(cs.residuals_resNeg$res) * 2, vertex.size = 2,
     vertex.label.cex = 0.7,  edge.curved=.2, vertex.color = "orange", edge.color = "dodgerblue3", 
     main = "Most negative residuals in sd alpha similarity")


```