Descriptive analyses

Last compiled on May, 2025

---

This is the code with which we render our descriptive analyses

---

1 Initatiating R environment

Start out with a custom function to load a set of required packages.

# packages and read data
rm(list = ls())

# (c) Jochem Tolsma
fpackage.check <- function(packages) {
    lapply(packages, FUN = function(x) {
        if (!require(x, character.only = TRUE)) {
            install.packages(x, dependencies = TRUE, repos = "http://cran.us.r-project.org")
            library(x, character.only = TRUE)
        }
    })
}
packages = c("haven", "coda", "matrixStats", "parallel", "MASS", "doParallel", "dplyr", "cowplot", "tidyverse",
    "naniar", "dotwhisker", "gt", "reshape2", "VGAM", "expss", "Hmisc", "poweRlaw", "fitdistrplus", "grid",
    "ggplotify", "smplot2")
fpackage.check(packages)

#> [[1]]
#> NULL
#> 
#> [[2]]
#> NULL
#> 
#> [[3]]
#> NULL
#> 
#> [[4]]
#> NULL
#> 
#> [[5]]
#> NULL
#> 
#> [[6]]
#> NULL
#> 
#> [[7]]
#> NULL
#> 
#> [[8]]
#> NULL
#> 
#> [[9]]
#> NULL
#> 
#> [[10]]
#> NULL
#> 
#> [[11]]
#> NULL
#> 
#> [[12]]
#> NULL
#> 
#> [[13]]
#> NULL
#> 
#> [[14]]
#> NULL
#> 
#> [[15]]
#> NULL
#> 
#> [[16]]
#> NULL
#> 
#> [[17]]
#> NULL
#> 
#> [[18]]
#> NULL
#> 
#> [[19]]
#> NULL
#> 
#> [[20]]
#> NULL
#> 
#> [[21]]
#> NULL

rm(packages)
load("data/dutch_netsize_analyses_revision_2.rda")

---

2 Descriptives independent variables

The code below generates the descriptive values which are used for Tables 1 and 2.

# descriptive table not automated: TABLE 1
table(df$work)

#> 
#>   0   1 
#> 583 666

psych::describe(df$work)

#>    vars    n mean  sd median trimmed mad min max range  skew kurtosis   se
#> X1    1 1249 0.53 0.5      1    0.54   0   0   1     1 -0.13    -1.98 0.01

psych::describe(df$hhsize)

#>    vars    n mean   sd median trimmed  mad min max range skew kurtosis   se
#> X1    1 1249 2.13 1.12      2    1.98 1.48   1   8     7 1.34     2.26 0.03

table(df$agecat)

#> 
#>   4   1   2   3 
#> 299 208 230 512

# 1 Maj 2 west 3 nonwest


table(df$income)

#> 
#>   2   1   3 
#> 440 627 182

# 1 < modal 2 > modal 3 unknown

psych::describe(df$worthhouse)

#>    vars    n mean   sd median trimmed  mad  min  max range skew kurtosis   se
#> X1    1 1249 2.67 0.89   2.54    2.58 0.73 0.95 9.96  9.01 1.75     7.12 0.03

table(df$woman)

#> 
#>   0   1 
#> 609 640

psych::describe(df$woman)

#>    vars    n mean  sd median trimmed mad min max range  skew kurtosis   se
#> X1    1 1249 0.51 0.5      1    0.52   0   0   1     1 -0.05       -2 0.01

table(df$opl2)

#> 
#>   0   1 
#> 783 466

summary(df$opl2)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.0000  0.0000  0.0000  0.3731  1.0000  1.0000

# 1 prim/sec 2 lower tert 3 higher tert

# psych::describe(df$neighdens)


# this is for the correlation table df$educ3 <- df$opl df$educ3 <-
# as.numeric(as.character(df$educ3))

table(df$migr3)

#> 
#>    1    2    3 
#> 1131   62   56

# cor(as.matrix(df[df$income!=3,c('income','worthhouse')]))

# rcorr(as.matrix(df[,c('work', 'hhsize', 'leeftijd10', 'woman', 'educ3', 'neighdens', 'worthhouse'
# )]), type = 'pearson')

# table 2
summary(df$perwomen)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.4800  0.6571  0.6108  0.7861  1.0000      11

summary(df[df$woman == 1, c("perwomen")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.5498  0.7148  0.6564  0.8158  1.0000       4

summary(df[df$woman == 0, c("perwomen")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.4223  0.5994  0.5626  0.7346  1.0000       7

summary(df$permen)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.2139  0.3429  0.3892  0.5200  1.0000      11

summary(df[df$woman == 1, c("permen")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.1842  0.2852  0.3436  0.4502  1.0000       4

summary(df[df$woman == 0, c("permen")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.2654  0.4006  0.4374  0.5777  1.0000       7

summary(df$samegender)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.3503  0.5709  0.5499  0.7590  1.0000      11

summary(df[df$woman == 1, c("samegender")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.5498  0.7148  0.6564  0.8158  1.0000       4

summary(df[df$woman == 0, c("samegender")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.2654  0.4006  0.4374  0.5777  1.0000       7

summary(df$pereduchigh)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>   0.000   0.400   0.750   0.652   1.000   1.000     640

summary(df[df$opl2 == 1, c("pereduchigh")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.6667  1.0000  0.8085  1.0000  1.0000     191

summary(df[df$opl2 == 0, c("pereduchigh")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.1952  0.5000  0.5231  1.0000  1.0000     449

summary(df$pereduclow)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>   0.000   0.000   0.250   0.348   0.600   1.000     640

summary(df[df$opl2 == 1, c("pereduclow")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.0000  0.0000  0.1915  0.3333  1.0000     191

summary(df[df$opl2 == 0, c("pereduclow")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.0000  0.5000  0.4769  0.8048  1.0000     449

summary(df$sameeduc)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.3333  0.6667  0.6266  1.0000  1.0000     640

summary(df[df$opl2 == 1, c("sameeduc")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.6667  1.0000  0.8085  1.0000  1.0000     191

summary(df[df$opl2 == 0, c("sameeduc")])

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>  0.0000  0.0000  0.5000  0.4769  0.8048  1.0000     449

summary(df$netsover6)

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    98.3   289.6   487.6   568.0   718.3  3590.9

---

3 Comparing naive estimations

The code below generates Figure 1 from the paper and runs a number of correlations found on page X.

# comps <- gather(df[, c('nsize_naive', 'nsize_b1b2', 'nsize_b3', 'nsize_n1', 'nsize_n2',
# 'nsize_n3', 'nsize_n123')]) comparisons <- ggplot(comps[comps$value < 2000, ], aes(x=value, color
# = key)) + geom_density(size = 0.5) + theme_minimal() + labs(x = 'Network size estimate', y =
# 'Density') + scale_color_manual(name = 'Estimation scenario', breaks = c('nsize_naive',
# 'nsize_b1b2', 'nsize_b3', 'nsize_n1', 'nsize_n2', 'nsize_n3', 'nsize_n123'), labels = c('Naive
# estimator: all X's', 'Battery 1: Tert. educ., elec vehicle, etc.', 'Battery 2: all names', 'Names
# 1: Sophie, Anna, Thomas, Willem', 'Names 2: Julia, Elisabeth, Max, Ali', 'Names 3: Sanne,
# Cornelia, Kevin, Mohammed', 'Names1+2+3'), values = c('#000000', '#E69F00', '#56B4E9', '#009E73',
# '#0072B2', '#D55E00', '#CC79A7', 'darkgrey')) # save ggsave('output/comparisons.pdf', plot =
# comparisons, device = 'pdf', scale = 1, width = 10, height = 5.5, units = c('in'), dpi =
# 'retina') # # correlation battery 1 and 2 # cor(df$nsize_b1, df$nsize_b2) # # # correlation
# battery 1 and all # cor(df$nsize_b1, df$nsize_naive) # correlation battery 1+2 and all
# cor(df$nsize_b1b2, df$nsize_naive) # correlation battery 3 (names) and all cor(df$nsize_b3,
# df$nsize_naive) # correlation names 1 and all cor(df$nsize_n1, df$nsize_naive) # correlation
# names 2 and all cor(df$nsize_n2, df$nsize_naive) # correlation names 3 and all cor(df$nsize_n3,
# df$nsize_naive) # correlation names 123 and all cor(df$nsize_n123, df$nsize_naive) # sample 5
# columns of scenarios and correlate with naive estimand set.seed(1987) sample(13:183, 5) #
# correlation scenario 85 and naive cor(df$nsize_naive, df[, 85]) # correlation scenario 36 and
# naive cor(df$nsize_naive, df[, 36]) # correlation scenario 150 and naive cor(df$nsize_naive, df[,
# 150]) # correlation scenario 14 and naive cor(df$nsize_naive, df[, 14]) # correlation scenario
# 125 and naive cor(df$nsize_naive, df[, 125]) comparisons

---

4 Comparing Bayesian estimations

The code below generates Figure 2 from the paper associated with the text on page X.

# netsizes <- read.table(file = 'data/dutch_netsize_desc.txt') # VIZ of netsize netsize_l <-
# gather(netsizes) dens1 <- ggplot(netsize_l[netsize_l$value < 2000,], aes(x=value, color = key)) +
# geom_density(alpha = .2, size = 0.1) + theme_minimal() + theme(legend.position = 'none') + labs(x
# = 'Network size estimate', y = 'Density') +
# geom_vline(xintercept=as.numeric(psych::describe(netsize_l[, 2])[3]), color = 'darkgrey',
# linetype = 2) + geom_vline(xintercept=as.numeric(psych::describe(netsize_l[, 2])[5]), color =
# 'darkgrey', linetype = 2) + annotate('text', x = 490, y = 0.0004, color = 'darkgrey', angle = 90,
# label = paste0('Mean = ', round(as.numeric(psych::describe(netsize_l[, 2])[3]), digits = 0))) +
# annotate('text', x = 340, y = 0.0004, color = 'darkgrey', angle = 90, label = paste0('Median = ',
# round(as.numeric(psych::describe(netsize_l[, 2])[5]), digits = 0))) + ggtitle('B) Distribution of
# network sizes, all scenarios') # Get lower triangle of the correlation matrix
# get_lower_tri<-function(cormat){ cormat[upper.tri(cormat)] <- NA return(cormat) } mat <-
# cor(netsizes) lower_tri <- get_lower_tri(mat) melted_cormat <- reshape2::melt(lower_tri) #
# data.table also has melt funciton that won't work on matrices melted_cormat$Var2 <-
# as.character(melted_cormat$Var2) melted_cormat$Var1 <- as.character(melted_cormat$Var1)
# melted_cormat <- melted_cormat[!melted_cormat$Var2 == melted_cormat$Var1,] melted_cormat <-
# melted_cormat[!is.na(melted_cormat$value), ] # Viz of correlations between netsize estimates
# dens2 <- ggplot(melted_cormat, aes(x=value)) + geom_density(size = 0.3) + xlim(0.915, 1) +
# theme_minimal() + labs(x='Pearson correlation', y = 'Density') + geom_vline(xintercept = .95,
# color = 'darkgrey', linetype = 2) + ggtitle('A) Distributions of correlations between network
# size scenarios') + geom_segment(x = 0.95, y = 125, xend = .99, yend = 125, linetype = 2, color =
# 'darkgrey', arrow = arrow(length = unit(0.25, 'cm'))) + annotate('text', x = .97, y = 135, color
# = 'darkgrey', angle = 0, label = paste0(round(nrow(melted_cormat[melted_cormat$value > .95,
# ])/nrow(melted_cormat)*100, digits = 1),'% of correlations are > .95')) # lay 'm out on the grid
# denses <- plot_grid(dens2, dens1, nrow = 1) # save ggsave('output/densities.pdf', plot = denses,
# device = 'pdf', scale = 1, width = 12, height = 4, units = c('in'), dpi = 'retina') # denses

---

5 Figure 1

The code below generates Figure 1 from the paper associated with the text on page X.

netsize_l <- data.frame(df[, c("netsover6")])
names(netsize_l) <- c("value")
round(summary(df$netsover6)[3], digits = 0)

#> Median 
#>    488

round(summary(df$netsover6)[4], digits = 0)

#> Mean 
#>  568

IQR(df$netsover6)

#> [1] 428.6466

dens1 <- ggplot(netsize_l, aes(x = value)) + geom_density(alpha = 0.2, size = 0.5) + theme_minimal() +
    theme(legend.position = "none") + labs(x = "Network size estimate", y = "Density") + scale_color_manual(values = c("#56B4E9")) +
    geom_vline(xintercept = 488, color = "darkgrey", linetype = 2) + geom_vline(xintercept = 568, color = "darkgrey",
    linetype = 2) + annotate("text", x = 630, y = 4e-04, color = "black", angle = 90, label = "Mean = 568") +
    annotate("text", x = 360, y = 4e-04, color = "black", angle = 90, label = "Median = 488") + xlim(c(0,
    2500))

ggsave("output/densities1_revision.pdf", plot = dens1, device = "pdf", scale = 1, width = 5, height = 4,
    units = c("in"), dpi = "retina")
ggsave("output/densities1_revision.jpeg", plot = dens1, device = "jpeg", scale = 1, width = 5, height = 4,
    units = c("in"), dpi = "retina")


cor(df[, c("netsover1", "netsover2", "netsover4", "netsover5", "netsover6")])

#>           netsover1 netsover2 netsover4 netsover5 netsover6
#> netsover1 1.0000000 0.9999047 0.9425976 0.9990875 0.8936700
#> netsover2 0.9999047 1.0000000 0.9428511 0.9990386 0.8939040
#> netsover4 0.9425976 0.9428511 1.0000000 0.9400682 0.9585224
#> netsover5 0.9990875 0.9990386 0.9400682 1.0000000 0.8913483
#> netsover6 0.8936700 0.8939040 0.9585224 0.8913483 1.0000000

# # VIZ of netsize scenarios netsize_l <- gather(df[, c('netsover1', 'netsover2', 'netsover4',
# 'netsover5', 'netsover6')]) round(summary(df$netsover6)[3], digits = 0)
# round(summary(df$netsover6)[4], digits = 0) dens2 <- ggplot(netsize_l, aes(x=value, color = key))
# + geom_density(alpha = .2, size = 0.5) + theme_minimal() + theme(legend.position = 'none') +
# labs(x = 'Network size estimate', y = 'Density') + scale_color_manual(values = c('grey', 'grey',
# 'grey', 'grey', '#56B4E9')) + geom_vline(xintercept=488, color = 'darkgrey', linetype = 2) +
# geom_vline(xintercept=568, color = 'darkgrey', linetype = 2) + annotate('text', x = 560, y =
# 0.0008, color = 'black', angle = 90, label = 'Mean = 568') + annotate('text', x = 400, y =
# 0.0008, color = 'black', angle = 90, label = 'Median = 488') + xlim(c(0, 2500)) # save
# ggsave('output/densities2_revision.pdf', plot = dens1, device = 'pdf', scale = 1, width = 5,
# height = 4, units = c('in'), dpi = 'retina') dens2

6 Distributional fit

We now fit a “complementary cumulative degree distribution” and check how it fits.

#df$netsize <- round(rowSums(df[,c(14:185)]) / length(14:185), 0)
df <- df[!df$netsover6 > 2500, ]
degrees <- round(df$netsover6, 0)

degr_pl <- displ$new(degrees) # create a discrete powerlaw distribution object 

# with estimating lower threshold:
degr_ln_xmin <- dislnorm$new(degrees)
est = estimate_xmin(degr_ln_xmin)
degr_ln_xmin$setXmin(est)

# without lower threshold:
degr_ln_noxmin <-fitdist(degrees, "lnorm")

# Goodness of fit for lognorm, no xmin:
gof_degr_ln_noxmin <- gofstat(degr_ln_noxmin, fitnames = "Log-normal", discrete = TRUE)
gof_degr_ln_noxmin

#> Chi-squared statistic:  397.1896 
#> Degree of freedom of the Chi-squared distribution:  26 
#> Chi-squared p-value:  4.716077e-68 
#> Chi-squared table:
#>         obscounts theocounts
#> <= 114   42.00000   19.04267
#> <= 155   42.00000   38.18254
#> <= 189   43.00000   47.59966
#> <= 220   41.00000   53.06053
#> <= 235   45.00000   28.07115
#> <= 277   42.00000   83.65007
#> <= 286   47.00000   18.48203
#> <= 329   42.00000   88.69643
#> <= 337   42.00000   16.36047
#> <= 380   42.00000   85.81643
#> <= 389   42.00000   17.37445
#> <= 427   41.00000   70.57309
#> <= 438   42.00000   19.51967
#> <= 448   42.00000   17.37170
#> <= 490   44.00000   68.97153
#> <= 505   42.00000   23.05357
#> <= 543   42.00000   54.71768
#> <= 554   44.00000   14.87218
#> <= 599   44.00000   56.50844
#> <= 634   41.00000   39.38943
#> <= 663   42.00000   29.83273
#> <= 718   43.00000   50.22114
#> <= 770   41.00000   40.66522
#> <= 820   41.00000   33.68975
#> <= 888   41.00000   38.58097
#> <= 981   41.00000   41.78179
#> <= 1111  41.00000   42.47278
#> <= 1297  41.00000   39.19962
#> > 1297   61.00000   66.24228
#> 
#> Goodness-of-fit criteria
#>                                Log-normal
#> Akaike's Information Criterion   17686.70
#> Bayesian Information Criterion   17696.95

x_seq <- seq(0, max(degrees), length = 500)

# For plotting the no x_min version of LN
ln_fit_ccdf <- 1 - plnorm(x_seq, meanlog = degr_ln_noxmin$estimate[1], sdlog =  degr_ln_noxmin$estimate[2])

# capture data that make the plot
pd1 <- plot(degr_pl, draw = F) # plot data 1, empirical distr
pd2 <- lines(degr_ln_xmin, draw = F) # plot data 2 no xmin
pd3 <- data.frame(cbind(x_seq, ln_fit_ccdf)) # plot data 3 with xmin

# now render it a ggplot
ccdfplot <- ggplot() + 
              geom_point(data=pd1, aes(x = x, y = y), shape=1) + 
              #geom_line(data=pd2, aes(x = x, y = y, color = "Log-normal fit, w/ x_min")) +
              geom_line(data=pd3, aes(x = x_seq, y = ln_fit_ccdf, color = "Log-normal fit, no x_min")) +
              labs(x="Extended network size", y="CCDF") + 
              scale_color_manual(name = "Fit", values = c(#"Log-normal fit, w/ x_min" = "green", 
                                                          "Log-normal fit, no x_min" = "purple")) +
              theme_bw() +
              theme(legend.position = c(0.22, 0.2)) +
              scale_x_log10(breaks = c(50,100,200,500,1000,2000), limits = c(min(pd1$x), max(pd1$x))) + scale_y_log10() +
              ggtitle("CCDF plot of average degree")

# save
ggsave("output/ccdfplot.pdf", plot = ccdfplot, device = "pdf",
       scale = 1, width = 6, height = 4, units = c("in"),
       dpi = "retina")
ggsave("output/ccdfplot.jpeg", plot = ccdfplot, device = "jpeg",
       scale = 1, width = 6, height = 4, units = c("in"),
       dpi = "retina")

# tryout <- plot_grid(dens2, dens1, ccdfplot, nrow = 1)
# 
# ggsave("output/comps_combined.pdf", plot = tryout, device = "pdf",
#        scale = 1, width = 16, height = 5, units = c("in"),
#        dpi = "retina")

### original plot in Base r
# # Plotting two versions of LN
# plot(degr_pl, xlab="Number of contacts", ylab="CCDF")
# lines(degr_ln_xmin, col = "green", lwd = 2)
# lines(x_seq, ln_fit_ccdf, col = "purple", lwd = 2)
# legend(x = "bottomleft",
#        legend = c("Log-normal fit, w/ x_min",
#                   "Log-normal fit, no x_min" 
#                   ),
#        col = c("green","purple"),
#        lwd = 2)
ccdfplot

dfans <- df[, c("uni", "hbo", "mbo", "dochterzoon", "tweeling", "corona", "elauto", "scooter", "vegan",
    "Sophie", "Julia", "Sanne", "Lisa", "Laura", "Maria", "Linda", "Johanna", "Monique", "Ester", "Anna",
    "Elisabeth", "Cornelia", "Wilhelmina", "Amira", "Samira", "Sara", "Daan", "Sem", "Thomas", "Max",
    "Kevin", "Johannes", "Dennis", "Jeroen", "Jan", "Marcel", "Cornelis", "Hendrik", "Petrus", "Willem",
    "Ali", "Mohammed", "Noor")]  # take mean of each of the cats
means <- colMeans(dfans)





ref1 <- c(84957, 75214, 145600, 168066, 2500, 1558549)
# ORDER: uni, hbo, mbo, dochter/zoon, tweeling, corona


ref2 <- c(273259, 460618, 261000)
# ORDER: elecauto, scooter, vegan,

ref3 <- c(15276, 16350, 27394, 21200, 25681, 334502, 29955, 266522, 39481, 2692, 136296, 110231, 112807,
    98208, 1386, 2186, 11640, 22704, 13276, 40543, 17024, 23167, 307032, 36411, 49182, 186746, 35973,
    134956, 118610, 86500, 102296, 4213, 5003, 4517)
# ORDER: Sophie, Julia,Sanne,Lisa,Laura,Maria,Linda,Johanna,Monique,Ester,
# Anna,Elisabeth,Cornelia,Wilhelmina,Amira,Samira,Sara,Daan,Sem,
# Thomas,Max,Kevin,Johannes,Dennis,Jeroen,Jan,Marcel,Cornelis,Hendrik,Petrus,
# Willem,Ali,Mohammed,Noor

pops <- c(ref1, ref2, ref3)

fig <- data.frame(cbind(means, pops))



cor(fig[!fig$pops > 5e+05, ])

#>           means      pops
#> means 1.0000000 0.3410926
#> pops  0.3410926 1.0000000

cor(fig[!fig$pops > 2e+05, ])

#>           means      pops
#> means 1.0000000 0.4616083
#> pops  0.4616083 1.0000000

cor(fig[fig$pops < 50000, ])

#>           means      pops
#> means 1.0000000 0.8918255
#> pops  0.8918255 1.0000000

cor(fig[fig$pops < 1e+05, ])

#>           means      pops
#> means 1.0000000 0.3374189
#> pops  0.3374189 1.0000000

ex3 <- ggplot(fig, aes(x = pops, y = means)) + geom_point() + xlim(x = c(0, 50000)) + ylim(c(0, 2.5)) +
    geom_smooth(method = "lm") + geom_text(label = rownames(fig), nudge_y = 0.1) + theme_minimal() +
    xlab("Population size") + ylab("Average number mentioned") + sm_statCorr() + ggtitle("C) Population size <50k")
ex3

ex2 <- ggplot(fig, aes(x = pops, y = means)) + geom_point() + xlim(x = c(0, 5e+05)) + ylim(c(0, 2.5)) +
    geom_smooth(method = "lm") + geom_text(label = rownames(fig), nudge_y = 0.1) + theme_minimal() +
    xlab("Population size") + ylab("Average number mentioned") + sm_statCorr() + ggtitle("B) Population size <500k")
ex2

ex1 <- ggplot(fig, aes(x = pops, y = means)) + geom_point() + geom_smooth(method = "lm") + geom_text(label = rownames(fig),
    nudge_y = 0.1) + theme_minimal() + xlab("Population size") + ylab("Average number mentioned") + sm_statCorr() +
    ggtitle("A) All population sizes")
ex1

ex <- cowplot::plot_grid(ex1, ex2, ex3, nrow = 1)

ggsave("output/popsizes.pdf", plot = ex, device = "pdf", scale = 1, width = 16, height = 5, units = c("in"),
    dpi = "retina")
ggsave("output/popsizes.jpeg", plot = ex, device = "jpeg", scale = 1, width = 16, height = 5, units = c("in"),
    dpi = "retina")
ex