---
title: "**Supplementary Material**"
subtitle: |
  \vspace{30pt}
  | Online Media Use and COVID-19 Vaccination in 
  | Real-World Personal Networks: Quantitative Study
  \vspace{10pt}
  | (Multimedia Appendix 1)



author: | 
   | \fcolorbox{green}{white}{\textbf{Iulian Oană}$^1$$^-$$^2$$^*$} \vspace{7pt} \hfill \phantom{phantomphantommm}
   | \fcolorbox{green}{white}{\textbf{Marian-Gabriel Hâncean}$^1$$^-$$^3$} \vspace{7pt} \hfill \phantom{phanto }
   | \fcolorbox{green}{white}{\textbf{Matjaž Perc}$^4$$^-$$^7$} \vspace{7pt} \hfill \phantom{phantomphantommm}
   | \fcolorbox{green}{white}{\textbf{Jürgen Lerner}$^8$} \vspace{7pt} \hfill \phantom{phantomphantommm}
   | \fcolorbox{green}{white}{\textbf{Bianca-Elena Mihăilă}$^1$$^-$$^2$} \vspace{7pt} \hfill \phantom{phantomph}
   | \fcolorbox{green}{white}{\textbf{Marius Geantă}$^2$} \vspace{7pt} \hfill \phantom{phantomphantommm}
   | \fcolorbox{green}{white}{\textbf{José Luis Molina}$^9$} \vspace{7pt} \hfill \phantom{phantomphantom}
   | \fcolorbox{green}{white}{\textbf{Isabela Tincă}$^1$$^-$$^2$} \vspace{7pt} \hfill \phantom{phantomphantomm}
   | \fcolorbox{green}{white}{\textbf{Carolina Espina}$^1$$^0$} \vspace{7pt} \hfill \phantom{phantomphantom}

# use the "date" for affiliations     
date: 
    \vspace{30pt}
    \begin{flushleft}
        \scriptsize 
        \textit{$^1$Department of Sociology, University of Bucharest, Panduri 90-92, 050663, Bucharest, Romania} \newline 
        \textit{$^2$Center for Innovation in Medicine, Th. Pallady 42J, 032266, Bucharest, Romania} \newline
        \textit{$^3$The Research Institute of the University of Bucharest, University of Bucharest, Panduri 90-92, 050663, Bucharest, Romania} \newline
        \textit{$^4$Faculty of Natural Sciences and Mathematics, University of Maribor, Koroška Cesta 160, 2000, Maribor, Slovenia} \newline
        \textit{$^5$Community Healthcare Center Dr. Adolf Drolc Maribor, Vošnjakova Ulica 2, 2000, Maribor, Slovenia} \newline
        \textit{$^6$Complexity Science Hub Vienna, Josefstädterstraße 39, 1080, Vienna, Austria} \newline
        \textit{$^7$Department of Physics, Kyung Hee University, 26 Kyungheedae-Ro, Dongdaemun-Gu, Seoul, Republic of Korea} \newline
        \textit{$^8$Department of Computer and Information Science, University of Konstanz, 78457, Konstanz, Germany} \newline
        \textit{$^9$GRAFO - Department of Social and Cultural Anthtropology, Universitat Autònoma de Barcelona, 08193, Bellaterra (Cerdanyola del Vallès), Barcelona, Spain} \newline
        \textit{$^1$$^0$International Agency for Research on Cancer(IARC/WHO), Environment and Lifestyle Epidemiology Branch, 25 avenue Tony Garnier CS 90627, 69366, Lyon CEDEX 07, France} \newline
      \newline  
      \newline
      \normalsize
      \text{$^*$Corresponding author; email:\fcolorbox{red}{white}{\textbf{\texttt{iulian.oana@sas.unibuc.ro}}}} \newline
    \end{flushleft}

    
header-includes:
   - \usepackage{fvextra}
   - \DefineVerbatimEnvironment{Highlighting}{Verbatim}{breaklines,commandchars=\\\{\}}
   - \renewcommand{\thetable}{S\arabic{table}}
   - \renewcommand{\thefigure}{S\arabic{figure}}

output: pdf_document
---

\newpage

<!-- Add table of contents -->

<!-- Add "\listoffigures" and "\listoftables" in the tex chunk if you want separate lists of figures and tables -->

\tableofcontents

\newpage

This is the Supplementary Material file corresponding to the manuscript titled: **Online Media Use and COVID-19 Vaccination in Real-World Personal Networks: Quantitative Study**.

This document contains the `R code` for replicating the results mentioned in the study. We divide this text into several sections corresponding to descriptive statistics, model estimates, model diagnostics, and miscellaneous.

### Libraries needed to run the code

```{r libraries, echo=TRUE, warning=FALSE, message=FALSE}
library(tidyverse)
library(knitr)
# for kableExtra a version lower than 1.4 is needed! Otherwise, the <scale_down> function might not work
library(kableExtra)
library(broom)
library(broom.mixed)
library(lme4)
library(summarytools)
library(report)
library(rempsyc)
library(Hmisc)
library(scales)
library(MASS)
library(jtools)
library(car)
library(sjPlot)
library(janitor)
library(miceadds)
library(pROC)
library(spsUtil)
library(eurostat)
library(stringi)
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}
# remove scientific notation
options(scipen = 999)
```

# **1. Load data**

\scriptsize

```{r load data}
# read ego data
ego_data <- readRDS("ego_data.rds") 

# read alter data
alter_data <- readRDS("alter_data.rds") 
```

\normalsize

## **1.1 About data frames:**

-   the data frames used in this script are divided into **ego data** and **alter data**.

-   ego data frames contain data about the egos and personal networks' structure and composition.

-   alter data frames contain data about the egos, alters, and personal networks' structure.

## **1.2 About variables:**

-   the suffix "**.f**" indicates columns where variables are stored as factors.

-   the suffix "**.mc**" indicates columns where variables are stored as scaled (mean centered) variables.

-   columns not ending in "**.f**" or "**.mc**" are numeric, storing the variables' original score.

-   the prefix "**prop\_**" refers to data describing proportions of certain features inside personal networks.

## **1.3 Variables' names and labels:**

### table S1

\scriptsize

```{r ego data variable labels, results='asis', warning=FALSE, message=FALSE}
#
variables_labels <-  data.frame(Hmisc::label(ego_data))
#
variables_labels <- variables_labels %>% rownames_to_column()
#
colnames(variables_labels) <- c("Variable", "Variable label")
#
variables_labels %>%
  kbl(caption = "Variable names and labels in ego data",
      booktabs = T, linesep = "") %>%
  kable_styling(full_width = T, latex_options = c("hold_position", "scale_down")) %>%
  column_spec(1, width = "5cm", italic = T) %>%
  row_spec(0, bold = T)

```

\normalsize

\newpage

### table S2

\scriptsize

```{r alter data variable labels, results='asis', warning=FALSE, message=FALSE}
#
variables_labels <-  data.frame(Hmisc::label(alter_data))
#
variables_labels <- variables_labels %>% rownames_to_column()
#
colnames(variables_labels) <- c("Variable", "Variable label")
#
variables_labels %>%
  kbl(caption = "Variable names and labels in alter data",
      booktabs = T, linesep = "") %>%
  kable_styling(full_width = F, latex_options = c("hold_position", "scale_down")) %>%
  column_spec(1, width = "5cm", italic = T) %>%
  row_spec(0, bold = T)

```

\normalsize

\newpage

# **2. Descriptive statistics**

## **2.1 Egos: numeric variables of interest and network level features**

### table S3

\scriptsize

```{r egos numeric descriptives, results='markup', warning=FALSE, message=FALSE}
#
ego_desc_num <- ego_data %>%
  # select variables of interest
  dplyr::select(ego.age, ego.education, 
                network_size, network_density, 
                network_components) %>%
  # compute summary statistics
  summarytools::descr() %>%
  round(2) %>%
  as.data.frame() %>%
  # reorder columns
  dplyr::select(1, 2, 5, 4, 3)

colnames(ego_desc_num) <- c("Egos' age", "Egos' education", "Network size",
                            "Network density", "Network components")

ego_desc_num %>%
  kbl(caption = "Descriptive statistics for numeric variables: 
      Ego characteristics and network level measures (level 2)", 
      booktabs = T, linesep = "") %>%
  kable_classic_2(full_width = T) %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

## **2.2 Alters: numeric variables of interest and ego-alter relation intensity**

### table S4

\scriptsize

```{r alter numeric descriptives, results='asis', warning=FALSE, message=FALSE}
#
alter_desc_num <- alter_data %>%
  dplyr::select(alter.age, alter.education, 
                ego.alter.intensity, alter.betweenness, alter.degree,
                vaccination.assortativity) %>%
  summarytools::descr() %>%
  round(2) %>%
  as.data.frame() %>%
  # reorder columns
  dplyr::select(1, 4, 5, 3, 2, 6)

colnames(alter_desc_num) <- c("Alters' age", "Alters' education", "Ego-Alter intensity",
                              "Alters' degree", "Alters' betweenness", "Vaccination assortativity")

alter_desc_num %>%
  kbl(caption = "Descriptive statistics for numeric variables: 
      Alter characteristics and ego-alter intensity (level 1)", 
      booktabs = T, linesep = "") %>%
  kable_classic_2(full_width = T) %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

## **2.3 Egos: categorical variables of interest**

### table S5

\scriptsize

```{r ego freq vacc status, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(ego_data$ego.vaccination.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Egos' vaccination status", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

\newpage

### table S6

\scriptsize

```{r ego freq sex, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(ego_data$ego.sex.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Egos' sex", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

### table S7

\scriptsize

```{r ego freq education, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(ego_data$ego.education.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Egos' education - last finished level", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

\newpage

### table S8

\scriptsize

```{r ego freq single, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(ego_data$ego.single.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Egos' being single", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

### table S9

\scriptsize

```{r ego freq employed, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(ego_data$ego.employed.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Egos' being employed", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

### table S10

\scriptsize

```{r ego freq media, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(ego_data$ego.media.use.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Egos' media use for information about health and prevention", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

## **2.4 Alters: categorical variables of interest**

### table S11

\scriptsize

```{r alter freq vacc status, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(alter_data$alter.vaccination.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Alters' vaccination status", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

### table S12

\scriptsize

```{r alter freq sex, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(alter_data$alter.sex.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Alters' sex", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

### table S13

\scriptsize

```{r alter freq education, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(alter_data$alter.education.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Alters' education - last finished level", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

### table S14

\scriptsize

```{r alter freq single, results='asis', warning=FALSE, message=FALSE}
#
results <- freq(alter_data$alter.single.f) %>% as.data.frame()
rownames(results)[rownames(results) == "<NA>"] <- "Missing"
results <- results %>% replace(is.na(.), 0)
#
results %>% dplyr::select(-c(2, 3)) %>%
  kbl(digits = 2,
      caption = "Alters' being single", 
      booktabs = T, linesep = "") %>%
  kable_classic_2() %>%
  column_spec(1, bold = T) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

# **3. Bivariate analyses**

## **3.1 Contingency table for egos' and alters' vaccination status**

### table S15

\scriptsize

```{r contingecy alter vacc by ego vacc, results='asis', warning=FALSE, message=FALSE}
#
results <- sjt.xtab(alter_data$ego.vaccination.f,
         alter_data$alter.vaccination.f,
         show.row.prc = T,
         show.col.prc = F,
         show.na = F,
         show.summary = F)
# convert results to data frame
results = sjtable2df::xtab2df(results)
# rename columns
colnames(results) <- c("Egos' vaccination status", "No", "Yes", "Total")
# print results
results %>% kbl(caption = "Contingency table for egos' and alters' vaccination status",
                booktabs = T, linesep = "") %>%
  add_header_above(c(" " = 1, "Alters' vaccination status" = 2, " " = 1)) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

#### chi-squared test for table S15

\scriptsize

```{r  chi-square alter vacc by ego vacc, results='asis', warning=FALSE, message=FALSE}
# print chi-squared results
#
chisqtable = (table(alter_data$ego.vaccination.f,
                    alter_data$alter.vaccination.f))
#
chisq.test(chisqtable)
```

\normalsize

## **3.2 Contingency table for egos' media use and alters' vaccination status**

### table S16

\scriptsize

```{r contingecy alter vacc by ego media, results='asis', warning=FALSE, message=FALSE}
#
results <- sjt.xtab(alter_data$ego.media.use.f,
         alter_data$alter.vaccination.f,
         show.row.prc = T,
         show.col.prc = F,
         show.na = F,
         show.summary = F)
# convert results to data frame
results = sjtable2df::xtab2df(results)
# rename columns
colnames(results) <- c("Egos' media use", "No", "Yes", "Total")
# print results
results %>% kbl(caption = "Contingency table for egos' media use and alters' vaccination status",
                booktabs = T, linesep = "") %>%
  add_header_above(c(" " = 1, "Alters' vaccination status" = 2, " " = 1)) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

#### chi-squared test for table S16

\scriptsize

```{r  chi-square alter vacc by ego media, results='asis', warning=FALSE, message=FALSE}
# print chi-squared results
#
chisqtable = (table(alter_data$ego.media.use.f,
                    alter_data$alter.vaccination.f))
#
chisq.test(chisqtable)
```

\normalsize

## **3.3 Contingency table for egos' media use and their vaccination status**

### table S17

\scriptsize

```{r contingecy ego vacc by ego media, results='asis', warning=FALSE, message=FALSE}
#
results <- sjt.xtab(ego_data$ego.media.use.f,
         ego_data$ego.vaccination.f,
         show.row.prc = T,
         show.col.prc = F,
         show.na = F,
         show.summary = F)
# convert results to data frame
results = sjtable2df::xtab2df(results)
# rename columns
colnames(results) <- c("Egos' media use", "No", "Yes", "Total")
# print results
results %>% kbl(caption = "Contingency table for egos' media use by their vaccination status",
                booktabs = T, linesep = "") %>%
  add_header_above(c(" " = 1, "Egos' vaccination status" = 2, " " = 1)) %>%
  kable_styling(latex_options = "hold_position")
```

\normalsize

#### chi-squared test for table S17

\scriptsize

```{r  chi-square ego vacc by ego media, results='asis', warning=FALSE, message=FALSE}
# print chi-squared results
#
chisqtable = (table(ego_data$ego.media.use.f,
                    ego_data$ego.vaccination.f))
#
chisq.test(chisqtable)
fisher.test(chisqtable)
```

\normalsize

## **3.3 T-test of independent samples: vaccination assortativity by alters' vaccination status**

### table S18

\scriptsize

```{r levene assortativity ttest1, results='asis', warning=FALSE, message=FALSE}
# test for variance by group
lt <- car::leveneTest(vaccination.assortativity ~ alter.vaccination.f, data = alter_data)
lt <- data.frame(lt) %>% mutate(df = lt[1,1]) %>% mutate(df2 = lt[2,1]) %>% filter(!is.na(F.value))
lt <- lt %>% dplyr::select(-c(1)) %>% rename("F value" = F.value, "P value" = Pr..F.)
row.names(lt) <- NULL

#
lt %>% 
  kbl(caption = "Levene test of variance equality", booktabs = T, digits = 5) %>% 
  kable_classic() %>%
  kable_styling(latex_options = "hold_position", full_width = F)
# variance equal: P>.05; variances unequal: P<0.05

```

\normalsize

### table S19

\scriptsize

```{r assortativity ttest1, results='asis', warning=FALSE, message=FALSE}
results <- t.test(vaccination.assortativity ~ alter.vaccination.f, 
       data = alter_data,
       var.equal = T)
# select necessary columns
t.table <- as.data.frame(report(results)) %>% 
  dplyr::select(3:11)
# rename columns
colnames(t.table) <- c("Mean - unvaccinated group", "Mean - vaccinated group",
                      "Mean difference", "CI", "CI low", "CI high", "t", "df", "p")
# print results
t.table %>% 
  kbl(digits = 4, 
      caption = "Vaccination assortativity by alters' vaccination status - independent samples t-test (two-sided)",
      booktabs = T, linesep = "") %>%
  kable_styling(latex_options = "hold_position", full_width = F) %>%
  # row_spec(0, align = "l") %>%
  column_spec(1, width = "3cm") %>% 
  column_spec(2, width = "2.7cm") %>% 
  row_spec(0, align = "l")
```

\normalsize

\newpage

## **3.4 T-test of independent samples: vaccination assortativity by egos' vaccination status**

### table S20

\scriptsize

```{r levene assortativity ttest2, results='asis', warning=FALSE, message=FALSE}
# test is variances are equal by group
lt <- car::leveneTest(vaccination.assortativity ~ ego.vaccination.f, data = alter_data)
lt <- data.frame(lt) %>% mutate(df = lt[1,1]) %>% mutate(df2 = lt[2,1]) %>% filter(!is.na(F.value))
lt <- lt %>% dplyr::select(-c(1)) %>% rename("F value" = F.value, "P value" = Pr..F.)
row.names(lt) <- NULL

#
lt %>% 
  kbl(caption = "Levene test of variance equality", booktabs = T, digits = 5) %>% 
  kable_classic() %>%
  kable_styling(latex_options = "hold_position", full_width = F)
# variance equal: P>.05; variances unequal: P<0.05
```

### table S21

```{r assortativity ttest2, results='asis', warning=FALSE, message=FALSE}
results <- t.test(vaccination.assortativity ~ ego.vaccination.f, 
       data = alter_data,
       var.equal = T)
# select necessary columns
t.table <- as.data.frame(report(results)) %>% 
  dplyr::select(3:11)
# rename columns
colnames(t.table) <- c("Mean - unvaccinated group", "Mean - vaccinated group",
                      "Mean difference", "CI", "CI low", "CI high", "t", "df", "p")
# print results
t.table %>% 
  kbl(digits = 4, 
      caption = "Vaccination assortativity by egos' vaccination status - independent samples t-test (two-sided)",
      booktabs = T, linesep = "") %>%
  kable_styling(latex_options = "hold_position", full_width = F) %>%
  # row_spec(0, align = "l") %>%
  column_spec(1, width = "3cm") %>% 
  column_spec(2, width = "2.7cm") %>% 
  row_spec(0, align = "l")
```

\normalsize

\newpage

## **3.5 Contingency table for alters' assortativity quartile by egos' media use x egos' vaccination status**

### table S22

\scriptsize

```{r chi-sq assort data, warning=FALSE, message=FALSE}
# create data frame with needed variables
my_tab <- alter_data %>%
  dplyr::select(networkCanvasEgoUUID, networkCanvasUUID, 
                ego.media.use.f, ego.vaccination.f,
                alter.vaccination.f, vaccination.assortativity) %>%
  # remove egos with missing data on the media use variable
  filter(!is.na(ego.media.use.f)) %>%
  # remove alters with missing data on the vaccination variable
  filter(!is.na(alter.vaccination.f)) %>%
  # remove alters with missing data on the vaccination assortativity variable
  filter(!is.na(vaccination.assortativity)) %>%
  # create variable combining egos' vaccination status and media use
  mutate(ego.media.use.f = as.character(ego.media.use.f)) %>%
  mutate(ego.vaccination.f = as.character(ego.vaccination.f)) %>%
  mutate(ego.media.vaccinated = paste(ego.media.use.f, ego.vaccination.f, sep = "-")) %>%
  # rename categories
  mutate(ego.media.vaccinated = str_replace_all(ego.media.vaccinated, "Both media-No",
                                       "Unvaccinated & uses both media")) %>%
    mutate(ego.media.vaccinated = str_replace_all(ego.media.vaccinated, "Both media-Yes",
                                       "Vaccinated & uses both media")) %>%
    mutate(ego.media.vaccinated = str_replace_all(ego.media.vaccinated, "Online media-No",
                                       "Unvaccinated & uses online media")) %>%
    mutate(ego.media.vaccinated = str_replace_all(ego.media.vaccinated, "Online media-Yes",
                                       "Vaccinated & uses online media")) %>%
    mutate(ego.media.vaccinated = str_replace_all(ego.media.vaccinated, "Traditional media-No",
                                       "Unvaccinated & uses traditional media")) %>%
    mutate(ego.media.vaccinated = str_replace_all(ego.media.vaccinated, "Traditional media-Yes",
                                       "Vaccinated & uses traditional media")) %>%
  mutate(ego.media.vaccinated = factor(ego.media.vaccinated)) %>%
  # create quartiles for alters' vaccination assortativity
  mutate(assortativity.quartile = ntile(vaccination.assortativity, 4))
```

```{r chi-sq assort table, warning=FALSE, message=FALSE}
# print results
results <- sjt.xtab(my_tab$ego.media.vaccinated,
         my_tab$assortativity.quartile,
         show.row.prc = T,
         show.col.prc = F,
         show.na = F,
         show.summary = F)
# convert results to data frame
results = sjtable2df::xtab2df(results)
# rename columns
colnames(results) <- c("Egos' vaccination status & media use", 
                       "Quartile 1", "Quartile 2", "Quartile 3", "Quartile 4", "Total")
# add needed spaces
results <- results %>% mutate(`Egos' vaccination status & media use` = str_replace_all(
  `Egos' vaccination status & media use`, "uses", "uses "
))
# print results
results %>% kbl(caption = "Alters' assortativity quartile by egos' vaccination status x media use",
                booktabs = T, linesep = " ") %>%
  add_header_above(c(" " = 1, "Alters' assortativity quartile" = 4, " " = 1)) %>%
  kable_styling(latex_options = c("hold_position", "scale down"), full_width = F, font_size = 7)
```

\normalsize

#### chi-squared test for table S22

\scriptsize

```{r  chi-square alter q by ego vacc-media, results='asis', warning=FALSE, message=FALSE}
# print chi-squared results
#
chisqtable = (table(my_tab$ego.media.vaccinated,
                    my_tab$assortativity.quartile))
#
chisq.test(chisqtable)
```

\normalsize

# **4. Multilevel logistic regression models estimation: alters-alter ties between friends**

### **4.1 Null model (Model 1)**

\scriptsize

```{r multilevel null M1, warning=FALSE, message=FALSE, results='asis'}
m1 <- lme4::glmer(
  alter.vaccination.f # Dependent variable (alter's vaccination status)
   ~ 1 +
    (1 | networkCanvasEgoUUID), # Intercept (1) varies in level-2 units (egos)
  family = binomial("logit"), # Model class (logistic)
  data = na.omit(alter_data),
  set.seed(1234),
  nAGQ = 100,
  control = glmerControl(optimizer = "bobyqa",
                         optCtrl = list(maxfun = 2e5))
)
```

\normalsize

### **4.2 Attributes model (Model 2)**

\scriptsize

```{r multilevel attributes M2, warning=FALSE, message=FALSE, results='asis'}
#
m2 <- lme4::glmer(
  alter.vaccination.f # Dependent variable (alter's vaccination status)
    
    ### alter attributes
    ~ alter.sex.f
    + alter.education.mc
    + alter.single.f
    + alter.age.mc
    
    ### ego attributes
    + ego.sex.f
    + ego.education.mc
    + ego.age.mc
    + ego.employed.f
    + ego.vaccination.f
    
    ### ego media use 
    + ego.media.use.f
    +
    (1 | networkCanvasEgoUUID), # Intercept (1) varies in level-2 units (egos)
  family = binomial("logit"), # Model class (logistic)
  data = na.omit(alter_data),
  set.seed(1234),
  nAGQ = 100,
  control = glmerControl(optimizer = "bobyqa",
                         optCtrl = list(maxfun = 2e5)))
```

\normalsize

### **4.3 Network model (Model 3)**

\scriptsize

```{r multilevel network M3, warning=FALSE, message=FALSE, results='asis'}
#
m3 <- lme4::glmer(
  alter.vaccination.f # Dependent variable (alter's vaccination status)
    
    ### network data: ego-alter relation
    ~ ego.alter.intensity
    
    ### network data: node level features
    + alter.betweenness.mc
    + vaccination.assortativity.mc
    
    ### network data: network level features
    + network_size.mc
    + network_density.mc
    + network_components.mc
    +
    (1 | networkCanvasEgoUUID), # Intercept (1) varies in level-2 units (egos)
  family = binomial("logit"), # Model class (logistic)
  data = na.omit(alter_data),
  set.seed(1234),
  nAGQ = 100,
  control = glmerControl(optimizer = "bobyqa",
                         optCtrl = list(maxfun = 2e5)))
```

\normalsize

### **4.4 Full model (Model 4)**

\scriptsize

```{r multilevel full M4, warning=FALSE, message=FALSE, results='asis'}
m4 <- lme4::glmer(
  alter.vaccination.f # Dependent variable (alter's vaccination status)
    
    ### alter attributes
    ~ alter.sex.f
    + alter.education.mc
    + alter.single.f
    + alter.age.mc
  
    ### ego attributes
    + ego.sex.f
    + ego.education.mc
    + ego.age.mc
    + ego.employed.f
    + ego.vaccination.f
  
    ### ego media use
    + ego.media.use.f
    
    ### network data: ego-alter relation
    + ego.alter.intensity
  
    ### network data: node level features
    + alter.betweenness.mc
    + vaccination.assortativity.mc
  
    ### network data: network level features
    + network_size.mc
    + network_density.mc
    + network_components.mc
    +
    
    (1 | networkCanvasEgoUUID), # Intercept (1) varies in level-2 units (egos)
  family = binomial("logit"), # Model class (logistic)
  data = na.omit(alter_data),
  set.seed(1234),
  nAGQ = 100,
  control = glmerControl(optimizer = "bobyqa",
                         optCtrl = list(maxfun = 2e5))
)
```

\normalsize

### **4.5 Print results**

\scriptsize

```{r multilevel tabel and labels, warning=FALSE, message=FALSE, results='asis'}
# Create data frame with all models
# -- create an  object with all the models
table.models = sjPlot::tab_model(m1, m2, m3, m4, 
                                     show.r2 = T, show.aic = T, show.loglik = T)
# -- transform the object into a data frame
table.models_df = sjtable2df::mtab2df(table.models, n_models = 4) %>% as.data.frame()
# Run code from "multilevel_labels.R" script to apply transformations to the data frame in order
# to recode labels and attach missing model fit estimates
source("multilevel_labels.R")
```

\normalsize

### table S23

\scriptsize

```{r multilevel models results, warning=FALSE, message=FALSE, results='asis'}

table.models_df %>% 
  kable(col.names = c("Predictors", "Odds Ratios", "95% CIs", "P value", "Odds Ratios", "95% CIs", "P value",
                      "Odds Ratios", "95% CIs", "P value", "Odds Ratios", "95% CIs", "P value"),
        caption = "Mixed multilevel regression models predicting alters' vaccination status",
        booktabs = T, linesep = "") %>%
  kable_classic() %>%
  add_header_above(c(" " = 1, "Model 1\n(null model)" = 3, 
                     "Model 2\n(attributes only)" = 3, 
                     "Model 3\n(network only)" = 3, 
                     "Model 4\n(full model)" = 3)) %>%
  kable_styling(latex_options = c("hold_position", "scale_down", "striped"), full_width = F) %>%
  row_spec(18,  hline_after = T)
```

\normalsize

\newpage

### figure S1

\scriptsize

```{r multilevel friends summary plot, fig.cap="Summary plot of multilevel models", fig.height=6.5, out.width="77%", fig.align='center', results='asis'}
#
jtools::plot_summs(m2, m3, m4,
                   coefs = c(
                    "Alter sex (female])" = "alter.sex.fFemale", "Alter education" = "alter.education.mc",
                    "Alter being single (yes)" = "alter.single.fYes", "Alter age" = "alter.age.mc",
                    "Ego sex (female)" = "ego.sex.fFemale", "Ego education" = "ego.education.mc",
                    "Ego age" = "ego.age.mc", "Ego being employed (yes)" = "ego.employed.fEmployed",
                    "Ego being vaccinated (yes)" = "ego.vaccination.fYes",
                    "Ego media use (online)" = "ego.media.use.fOnline media",
                    "Ego media use (both)" = "ego.media.use.fBoth media",
                    "Ego-alter intensity" = "ego.alter.intensity", "Alter betweenness" = "alter.betweenness.mc",
                    "Vaccination assortativity" = "vaccination.assortativity.mc",
                    "Network size" = "network_size.mc", "Network density" = "network_density.mc",
                    "Network components (strong)" = "network_components.mc"),
                   exp = TRUE,
                   model.names = c("Model 2\n(attributes only)", 
                                   "Model 3\n(network only)", 
                                   "Model 4\n(full model)")) +
  theme(plot.margin = margin(t = 0.75,  # Top margin
                             r = 0,  # Right margin
                             b = 0,  # Bottom margin
                             l = 0,  # Left margin
                             unit = "cm"),
        axis.title.x = element_text(margin = margin(t = 10), 
                                    hjust = 0.43, family = "Times"),
        # text=element_text(family="Times"),
        axis.text.x = element_text(family="Times", face = "bold"),
        axis.text.y = element_text(family="Times", hjust = 0),
        legend.text = element_text(family="Times", face = "bold", hjust = 0.5),
        legend.title = element_text(family="Times")) +
  labs(x = "Odds ratios & 95% CIs")
```

\normalsize

\newpage

# **5. Model diagnostics**

## **5.1 Predicting alters' vaccination status using GLM with robust clustered errors**

### **5.1.1 GLM model 1**

\scriptsize

```{r GLM M1, warning=FALSE, message=FALSE, results='asis'}
#
m1.glm <- miceadds::glm.cluster(alter.vaccination.f # Dependent variable (alter's vaccination status)
              ~ 1
              ### control for proportion of unvaccinated alters (minus the alter of interest)
              + prop.vacc.alters_ex.alter.mc, 
              
              ### cluster errors by ego's ID & set random number seed for reproducibility of results
              cluster="networkCanvasEgoUUID",
              data = na.omit(alter_data),
              set.seed(1234),
              family = "binomial"
)
```

\normalsize

### **5.1.2 GLM model 2**

\scriptsize

```{r GLM M2, warning=FALSE, message=FALSE, results='asis'}
#
m2.glm <- miceadds::glm.cluster(alter.vaccination.f # Dependent variable (alter's vaccination status)
    ### alter attributes
    ~ alter.sex.f
    + alter.education.mc
    + alter.single.f
    + alter.age.mc
  
    ### ego attributes
    + ego.sex.f
    + ego.education.mc
    + ego.age.mc
    + ego.employed.f
    + ego.vaccination.f
  
    ### ego media use 
    + ego.media.use.f
    
    ### control for proportion of unvaccinated alters (minus the alter of interest)
    + prop.vacc.alters_ex.alter.mc,
    
    ### cluster errors by ego's ID & set random number seed for reproducibility of results
    cluster="networkCanvasEgoUUID",
    data = na.omit(alter_data),
    set.seed(1234),
    family = "binomial"
)
```

\normalsize

### **5.1.3 GLM model 3**

\scriptsize

```{r GLM M3, warning=FALSE, message=FALSE, results='asis'}
#
m3.glm <- miceadds::glm.cluster(alter.vaccination.f # Dependent variable (alter's vaccination status)
    ### network data: ego-alter relation
    ~ ego.alter.intensity
  
    ### network data: node level features
    + alter.betweenness.mc
    + vaccination.assortativity.mc
  
    ### network data: network level features
    + network_size.mc
    + network_density.mc
    + network_components.mc

    ### control for proportion of unvaccinated alters (minus the alter of interest)
    + prop.vacc.alters_ex.alter.mc,
    
    # cluster errors by ego's ID & set random number seed for reproducibility of results
    cluster="networkCanvasEgoUUID",
    data = na.omit(alter_data),
    set.seed(1234),
    family = "binomial"
)
```

\normalsize

### **5.1.4 GLM model 4**

\scriptsize

```{r GLM M4, warning=FALSE, message=FALSE, results='asis'}
#
m4.glm <- miceadds::glm.cluster(alter.vaccination.f # Dependent variable (alter's vaccination status)
    ### alter attributes
    ~ alter.sex.f
    + alter.education.mc
    + alter.single.f
    + alter.age.mc
  
    ### ego attributes
    + ego.sex.f
    + ego.education.mc
    + ego.age.mc
    + ego.employed.f
    + ego.vaccination.f

    ### ego media use 
    + ego.media.use.f
              
    ### network data: ego-alter relation
    + ego.alter.intensity
  
    ### network data: node level features
    + alter.betweenness.mc
    + vaccination.assortativity.mc
  
    ### network data: network level features
    + network_size.mc
    + network_density.mc
    + network_components.mc

    ### control for proportion of unvaccinated alters (minus the alter of interest)
    + prop.vacc.alters_ex.alter.mc,
    
    ### cluster errors by ego's ID & set random number seed for reproducibility of results
    cluster="networkCanvasEgoUUID",
    data = na.omit(alter_data),
    set.seed(1234),
    family = "binomial"
)
```

\normalsize

### **5.1.5 Print results**

\scriptsize

```{r GLM tabel and labels, warning=FALSE, message=FALSE, results='asis'}
# Run code from "glm_labels.R" script to apply transformations to the data frames in order to create a single table, recode labels, and attach missing model fit estimates
source("glm_labels.R")
```

\normalsize

### table S24

\scriptsize

```{r glm models results, warning=FALSE, message=FALSE, results='asis'}
all.glm %>% 
  kable(col.names = c("Predictors",
                      "Odds Ratios",
                      "95% CIs",
                      "P value",
                      "Odds Ratios",
                      "95% CIs",
                      "P value",
                      "Odds Ratios",
                      "95% CIs",
                      "P value",
                      "Odds Ratios",
                      "95% CIs",
                      "P value"),
        caption = "GLM regression models predicting alters' vaccination status",
        booktabs = T, linesep = "") %>%
  kable_classic() %>%
  add_header_above(c(" " = 1, "GLM Model 1\n(null model)" = 3, 
                     "GLM Model 2\n(attributes only)" = 3, 
                     "GLM Model 3\n(network only)" = 3, 
                     "GLM Model 4\n(full model)" = 3)) %>%
  row_spec(19,  hline_after = T) %>%
  kable_styling(latex_options = c("hold_position", "scale_down", "striped"), full_width = F)
```

\normalsize

\newpage

## **5.2 Multicolliniarity**

### **5.2.1 VIF scores for multilevel models**

### table S25

\scriptsize

```{r vif m2, warning=FALSE, message=FALSE, results='asis'}
# vif scores for multilevel model 2
vif_values_m2 <- vif(m2)
vif_values_m2 %>% kbl(caption = "VIF scores for multilevel regression model 2 (attributes only)",
                       booktabs = T, linesep = "") %>% 
  kable_classic() %>% kable_styling(latex_options = "hold_position", full_width = F)
```

\normalsize

### table S26

\scriptsize

```{r vif m3, warning=FALSE, message=FALSE, results='asis'}
# vif scores for multilevel model 3
vif_values_m3 <- vif(m3)
vif_values_m3 <- data.frame(vif_values_m3) %>%
  # rename first column and create second column with degrees of freedom (all 1)
  rename(GVIF = vif_values_m3) %>% mutate(Df = 1)
vif_values_m3 %>% kbl(caption = "VIF scores for multilevel regression model 3 (network only)",
                      booktabs = T, linesep = "") %>% 
  kable_classic() %>% kable_styling(latex_options = "hold_position", full_width = F)
```

\normalsize

### table S27

\scriptsize

```{r vif m4, warning=FALSE, message=FALSE, results='asis'}
# vif scores for multilevel model 4
vif_values_m4 <- vif(m4)
vif_values_m4 %>% kbl(caption = "VIF scores for multilevel regression model 4 (full model)",
                      booktabs = T, linesep = "") %>% 
  kable_classic() %>% kable_styling(latex_options = "hold_position", full_width = F)
```

\normalsize

\newpage

### **5.2.2 VIF scores for GLM models with clustered errors**

### table S28

\scriptsize

```{r vif m2.glm, warning=FALSE, message=FALSE, results='asis'}
# vif scores for GLM model 2
vif_values_m2.glm <- vif(m2.glm$glm_res)
vif_values_m2.glm %>% kbl(caption = "VIF scores for GLM regression model 2 (attributes only)",
                          booktabs = T, linesep = "") %>% 
  kable_classic() %>% kable_styling(latex_options = "hold_position", full_width = F)
```

\normalsize

### table S29

\scriptsize

```{r vif m3.glm, warning=FALSE, message=FALSE, results='asis'}
# vif scores for GLM model 3
vif_values_m3.glm <- vif(m3.glm$glm_res)
vif_values_m3.glm <- data.frame(vif_values_m3.glm) %>%
  # rename first column and create second column with degrees of freedom (all 1)
  rename(GVIF = vif_values_m3.glm) %>% mutate(Df = 1)
vif_values_m3.glm %>% kbl(caption = "VIF scores for GLM regression model 3 (network only)",
                          booktabs = T, linesep = "") %>% 
  kable_classic() %>% kable_styling(latex_options = "hold_position", full_width = F)
```

\normalsize

### table S30

\scriptsize

```{r vif m4.glm, warning=FALSE, message=FALSE, results='asis'}
# vif scores for GLM model 4
vif_values_m4.glm <- vif(m4.glm$glm_res)
vif_values_m4.glm %>% kbl(caption = "VIF scores for GLM regression model 4 (full model)",
                          booktabs = T, linesep = "") %>% 
  kable_classic() %>% kable_styling(latex_options = "hold_position", full_width = F)
```

\normalsize

## **5.3 Predictive performance**

### **5.3.1 ROC curve for multilevel models**

### figure S2

\scriptsize

```{r roc m1 plot, fig.cap="ROC curve for multilevel model 1 (null model)", fig.height=5, out.width="55%", fig.align='center', results='asis', warning=FALSE, message=FALSE}
# roc_curve multilevel model 1 (m1)
data_m1 <- getData(m1)  
predicted_probs_m1 <- predict(m1, type = "response")
roc_curve <- roc(data_m1$alter.vaccination.f, predicted_probs_m1)
auc_value <- auc(roc_curve)
plot(roc_curve, main = sprintf("ROC Curve (AUC = %.2f)", auc_value))
```

\normalsize

### figure S3

\scriptsize

```{r roc m2 plot, fig.cap="ROC curve for multilevel model 2 (attributes only)", fig.height=5, out.width="55%", fig.align='center', results='asis', warning=FALSE, message=FALSE}
# roc_curve multilevel model 2 (m2)
data_m2 <- getData(m2)  
predicted_probs_m2 <- predict(m2, type = "response")
roc_curve <- roc(data_m2$alter.vaccination.f, predicted_probs_m2)
auc_value <- auc(roc_curve)
plot(roc_curve, main = sprintf("ROC Curve (AUC = %.2f)", auc_value))
```

\normalsize

### figure S4

\scriptsize

```{r roc m3 plot, fig.cap="ROC curve for multilevel model 3 (network only)", fig.height=5, out.width="55%", fig.align='center', results='asis', warning=FALSE, message=FALSE}
# roc_curve multilevel model 3 (m3)
data_m3 <- getData(m3)  
predicted_probs_m3 <- predict(m3, type = "response")
roc_curve <- roc(data_m3$alter.vaccination.f, predicted_probs_m3)
auc_value <- auc(roc_curve)
plot(roc_curve, main = sprintf("ROC Curve (AUC = %.2f)", auc_value))
```

\normalsize

### figure S5

\scriptsize

```{r roc m4 plot, fig.cap="ROC curve for multilevel model 4 (full model)", fig.height=5, out.width="55%", fig.align='center', results='asis', warning=FALSE, message=FALSE}
# roc_curve multilevel model 4 (m4)
data_m4 <- getData(m4)  
predicted_probs_m4 <- predict(m4, type = "response")
roc_curve <- roc(data_m4$alter.vaccination.f, predicted_probs_m4)
auc_value <- auc(roc_curve)
plot(roc_curve, main = sprintf("ROC Curve (AUC = %.2f)", auc_value))
```

\normalsize

### 5.3.1 ROC curve for GLM models with clustered errors

### figure S6

\scriptsize

```{r roc m1.glm plot, fig.cap="ROC curve for GLM model 1 (null model)", fig.height=5, out.width="50%", fig.align='center', results='asis', warning=FALSE, message=FALSE}
# roc_curve GLM model 1 (m1.glm)
data_m1.glm <- m1.glm$glm_res$data  
predicted_probs_m1.glm <- predict(m1.glm$glm_res, type = "response")
roc_curve <- roc(data_m1.glm$alter.vaccination.f, predicted_probs_m1.glm)
auc_value <- auc(roc_curve)
plot(roc_curve, main = sprintf("ROC Curve (AUC = %.2f)", auc_value))
```

\normalsize

### figure S7

\scriptsize

```{r roc m2.glm plot, fig.cap="ROC curve for GLM model 2 (attributes only)", fig.height=5, out.width="55%", fig.align='center', results='asis', warning=FALSE, message=FALSE}

# roc_curve GLM model 2 (m2.glm)
data_m2.glm <- m2.glm$glm_res$data  
predicted_probs_m2.glm <- predict(m2.glm$glm_res, type = "response")
roc_curve <- roc(data_m2.glm$alter.vaccination.f, predicted_probs_m2.glm)
auc_value <- auc(roc_curve)
plot(roc_curve, main = sprintf("ROC Curve (AUC = %.2f)", auc_value))
```

\normalsize

### figure S8

\scriptsize

```{r roc m3.glm plot, fig.cap="ROC curve for GLM model 3 (network only)", fig.height=5, out.width="55%", fig.align='center', results='asis', warning=FALSE, message=FALSE}
# roc_curve GLM model 3 (m3.glm)
data_m3.glm <- m3.glm$glm_res$data  
predicted_probs_m3.glm <- predict(m3.glm$glm_res, type = "response")
roc_curve <- roc(data_m3.glm$alter.vaccination.f, predicted_probs_m3.glm)
auc_value <- auc(roc_curve)
plot(roc_curve, main = sprintf("ROC Curve (AUC = %.2f)", auc_value))
```

\normalsize

### figure S9

\scriptsize

```{r roc m4.glm plot, fig.cap="ROC curve for GLM model 4 (full model)", fig.height=5, out.width="55%", fig.align='center', results='asis', warning=FALSE, message=FALSE}
# roc_curve GLM model 4 (m4.glm)
data_m4.glm <- m4.glm$glm_res$data  
predicted_probs_m4.glm <- predict(m4.glm$glm_res, type = "response")
roc_curve <- roc(data_m4.glm$alter.vaccination.f, predicted_probs_m4.glm)
auc_value <- auc(roc_curve)
plot(roc_curve, main = sprintf("ROC Curve (AUC = %.2f)", auc_value))
```

\normalsize

\newpage

## **5.4 Residual plots**

### **5.4.1 Residual plots for multilevel models**

### figure S10

\scriptsize

```{r res plot m1, fig.cap="Residuals plot for multilevel model 1 (null model)", fig.height=5, out.width="51%", fig.align='center', results='asis'}

# multilevel model 1 (m1)
fitted_values_m1 <- fitted(m1)
residuals_m1 <- resid(m1)
data_frame <- data.frame(Fitted = fitted_values_m1, Residuals = residuals_m1)
ggplot(data_frame, aes(x = Fitted, y = Residuals)) + geom_point() +geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() + xlab("Fitted Values") + ylab("Residuals")
```

\normalsize

### figure S11

\scriptsize

```{r res plot m2, fig.cap="Residuals plot for multilevel model 2 (attributes only)", fig.height=5, out.width="51%", fig.align='center', results='asis'}
# multilevel model 2 (m2)
fitted_values_m2 <- fitted(m2)
residuals_m2 <- resid(m2)
data_frame <- data.frame(Fitted = fitted_values_m2, Residuals = residuals_m2)
ggplot(data_frame, aes(x = Fitted, y = Residuals)) + geom_point() + geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() + xlab("Fitted Values") + ylab("Residuals")
```

\normalsize

### figure S12

\scriptsize

```{r res plot m3, fig.cap="Residuals plot for multilevel model 3 (network only)", fig.height=5, out.width="51%", fig.align='center', results='asis'}
# multilevel model 3 (m3)
fitted_values_m3 <- fitted(m3)
residuals_m3 <- resid(m3)
data_frame <- data.frame(Fitted = fitted_values_m3, Residuals = residuals_m3)
ggplot(data_frame, aes(x = Fitted, y = Residuals)) + geom_point() + geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() + xlab("Fitted Values") + ylab("Residuals")
```

\normalsize

### figure S13

\scriptsize

```{r res plot m4, fig.cap="Residuals plot for multilevel model 4 (full model)", fig.height=5, out.width="51%", fig.align='center', results='asis'}
# multilevel model 4 (m4)
fitted_values_m4 <- fitted(m4)
residuals_m4 <- resid(m4)
data_frame <- data.frame(Fitted = fitted_values_m4, Residuals = residuals_m4)
ggplot(data_frame, aes(x = Fitted, y = Residuals)) + geom_point() + geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() + xlab("Fitted Values") + ylab("Residuals")
```

\normalsize

### **5.4.2 Residual plots for GLM models with clustered errors**

### figure S14

\scriptsize

```{r res plot m1.glm, fig.cap="Residuals plot for GLM model 1 (null model)", fig.height=5, out.width="51%", fig.align='center', results='asis'}
# GLM model 1 (m1.glm)
fitted_values_m1.glm <- fitted(m1.glm$glm_res)
residuals_m1.glm <- resid(m1.glm$glm_res)
data_frame <- data.frame(Fitted = fitted_values_m1.glm, Residuals = residuals_m1.glm)
ggplot(data_frame, aes(x = Fitted, y = Residuals)) + geom_point() + geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() + xlab("Fitted Values") + ylab("Residuals")
```

\normalsize

### figure S15

\scriptsize

```{r res plot m2.glm, fig.cap="Residuals plot for GLM model 2 (attributes only)", fig.height=5, out.width="51%", fig.align='center', results='asis'}

# GLM model 2 (m2.glm)
fitted_values_m2.glm <- fitted(m2.glm$glm_res)
residuals_m2.glm <- resid(m2.glm$glm_res)
data_frame <- data.frame(Fitted = fitted_values_m2.glm, Residuals = residuals_m2.glm)
ggplot(data_frame, aes(x = Fitted, y = Residuals)) + geom_point() + geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() + xlab("Fitted Values") + ylab("Residuals")
```

\normalsize

### figure S16

\scriptsize

```{r res plot m3.glm, fig.cap="Residuals plot for GLM model 3 (network only)", fig.height=5, out.width="51%", fig.align='center', results='asis'}

# GLM model 3 (m3.glm)
fitted_values_m3.glm <- fitted(m3.glm$glm_res)
residuals_m3.glm <- resid(m3.glm$glm_res)
data_frame <- data.frame(Fitted = fitted_values_m3.glm, Residuals = residuals_m3.glm)
ggplot(data_frame, aes(x = Fitted, y = Residuals)) + geom_point() + geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() + xlab("Fitted Values") + ylab("Residuals")
```

\normalsize

### figure S17

\scriptsize

```{r res plot m4.glm, fig.cap="Residuals plot for GLM model 4 (full model)", fig.height=5, out.width="51%", fig.align='center', results='asis'}
# GLM model 4 (m4.glm)
fitted_values_m4.glm <- fitted(m4.glm$glm_res)
residuals_m4.glm <- resid(m4.glm$glm_res)
data_frame <- data.frame(Fitted = fitted_values_m4.glm, Residuals = residuals_m4.glm)
ggplot(data_frame, aes(x = Fitted, y = Residuals)) + geom_point() + geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() + xlab("Fitted Values") + ylab("Residuals")
```

\normalsize

## **5.5 Overdispersion estimation**

### **5.5.1 Residual plots for multilevel models**

### table S31

\scriptsize

```{r overdispersion test m1, warning=FALSE, message=FALSE, results='asis'}
# multilevel model 1 (m1)
print(DHARMa::testDispersion(m1, plot = F))
```

\normalsize

### figure S18

\scriptsize

```{r overdispersion plot m1, fig.cap="Overdispersion plots for multilevel model 1 (null model)", fig.width=9, fig.height=6.5, out.width="100%", fig.align='center', results='asis'}
# multilevel model 1 (m1)
plot(DHARMa::simulateResiduals(m1))
```

\normalsize

\newpage

### table S32

\scriptsize

```{r overdispersion test m2, warning=FALSE, message=FALSE, results='asis'}
# multilevel model 2 (m2)
print(DHARMa::testDispersion(m2, plot = F))
```

\normalsize

### figure S19

\scriptsize

```{r overdispersion plot m2, fig.cap="Overdispersion plots for multilevel model 2 (attributes only)", fig.width=9, fig.height=6.5, out.width="100%", fig.align='center', results='asis'}
# multilevel model 2 (m2)
plot(DHARMa::simulateResiduals(m2))
```

\normalsize

\newpage

### table S33

\scriptsize

```{r overdispersion test m3, warning=FALSE, message=FALSE, results='asis'}
# multilevel model 3 (m3)
print(DHARMa::testDispersion(m3, plot = F))
```

\normalsize

### figure S20

\scriptsize

```{r overdispersion plot m3, fig.cap="Overdispersion plots for multilevel model 3 (network only)", fig.width=9, fig.height=6.5, out.width="100%", fig.align='center', results='asis'}
# multilevel model 3 (m3)
plot(DHARMa::simulateResiduals(m3))
```

\normalsize

\newpage

### table S34

\scriptsize

```{r overdispersion test m4, warning=FALSE, message=FALSE, results='asis'}
# multilevel model 4 (m4)
print(DHARMa::testDispersion(m4, plot = F))
```

\normalsize

### figure S21

\scriptsize

```{r overdispersion plot m4, fig.cap="Overdispersion plots for multilevel model 4 (full model)",fig.width=9, fig.height=6.5, out.width="100%", fig.align='center', results='asis'}
# multilevel model 4 (m4)
plot(DHARMa::simulateResiduals(m4))
```

\normalsize

\newpage

### **5.5.2 Residual plots for GLM models with clustered errors**

### table S35

\scriptsize

```{r overdispersion test m1.glm, warning=FALSE, message=FALSE, results='asis'}
# GLM model 1 (m1.glm)
print(DHARMa::testDispersion(m1.glm$glm_res, plot = F))
```

\normalsize

### figure S22

\scriptsize

```{r overdispersion plot m1.glm, fig.cap="Overdispersion plots for GLM model 1 (null model)",fig.width=9, fig.height=6.5, out.width="100%", fig.align='center', results='asis'}
# GLM model 1 (m1.glm)
plot(DHARMa::simulateResiduals(m1.glm$glm_res))
```

\normalsize

\newpage

### table S36

\scriptsize

```{r overdispersion test m2.glm, warning=FALSE, message=FALSE, results='asis'}
# GLM model 2 (m2.glm)
print(DHARMa::testDispersion(m2.glm$glm_res, plot = F))
```

\normalsize

### figure S23

\scriptsize

```{r overdispersion plot m2.glm, fig.cap="Overdispersion plots for GLM model 2 (attributes only)", fig.width=9, fig.height=6.5, out.width="100%", fig.align='center', results='asis'}
# GLM model 2 (m2.glm)
plot(DHARMa::simulateResiduals(m2.glm$glm_res))
```

\normalsize

\newpage

### table S37

\scriptsize

```{r overdispersion test m3.glm, warning=FALSE, message=FALSE, results='asis'}
# GLM model 3 (m3.glm)
print(DHARMa::testDispersion(m3.glm$glm_res, plot = F))
```

\normalsize

### figure S24

\scriptsize

```{r overdispersion plot m3.glm, fig.cap="Overdispersion plots for GLM model 3 (network only)", fig.width=9, fig.height=6.5, out.width="100%", fig.align='center', results='asis'}
# GLM model 3 (m3.glm)
plot(DHARMa::simulateResiduals(m3.glm$glm_res))
```

\normalsize

\newpage

### table S38

\scriptsize

```{r overdispersion test m4.glm, warning=FALSE, message=FALSE, results='asis'}
# GLM model 4 (m4.glm)
print(DHARMa::testDispersion(m4.glm$glm_res, plot = F))
```

\normalsize

### figure S25

\scriptsize

```{r overdispersion plot m4.glm, fig.cap="Overdispersion plots for GLM model 4 (full model)", fig.width=9, fig.height=6.5, out.width="100%", fig.align='center', results='asis'}
# GLM model 4 (m4.glm)
plot(DHARMa::simulateResiduals(m4.glm$glm_res))
```

\normalsize

\newpage

## **5.6 Plots of the residuals for random effects in multilevel models**

### figure S26

\scriptsize

```{r res random eff plot m1, fig.cap="Residuals for random effects plots for multilevel model 1 (null model)", fig.height=6.5, out.width="77%", fig.align='center', results='asis'}
# model 1 (m1)
ranef_resid <- ranef(m1, condVar = TRUE)
plot(ranef_resid, which = "random_effect")
```

\normalsize

\newpage

### figure S27

\scriptsize

```{r res random eff plot m2, fig.cap="Residuals for random effects plots for multilevel model 2 (attributes only)", fig.height=6.5, out.width="77%", fig.align='center', results='asis'}
# model 2 (m2)
ranef_resid <- ranef(m2, condVar = TRUE)
plot(ranef_resid, which = "random_effect")
```

\normalsize

\newpage

### figure S28

\scriptsize

```{r res random eff plot m3, fig.cap="Residuals for random effects plots for multilevel model 3 (network only)", fig.height=6.5, out.width="77%", fig.align='center', results='asis'}
# model 3 (m3)
ranef_resid <- ranef(m3, condVar = TRUE)
plot(ranef_resid, which = "random_effect")
```

\normalsize

\newpage

### figure S29

\scriptsize

```{r res random eff plot m4, fig.cap="Residuals for random effects plots for multilevel model 4 (full model)", fig.height=6.5, out.width="77%", fig.align='center', results='asis'}
# model 4 (m4)
ranef_resid <- ranef(m4, condVar = TRUE)
plot(ranef_resid, which = "random_effect")
```

\normalsize

\newpage

# **6. Miscellaneous (Eurostat data)**

We mention that, as these data are dynamically downloaded using the `get_eurostat()` function from the `eurostat` package, slight changes in percentages might appear, when comparing the values from the article with those from the subsequent figures.

\scriptsize

```{r, eu c labels, warning=FALSE, message=FALSE}
# Get labels for EU countries

eu_c.labels <- eurostat::eu_countries %>% rename(geo = code, geo.name = name) %>% dplyr::select(geo, geo.name)
```

\normalsize

## **6.1 Eurostat data from dataset "Individuals' level of digital skills (from 2021 onwards)"**

**`Online data code`**`: isoc_sk_dskl_i21`

**`link`**`: https://ec.europa.eu/eurostat/databrowser/view/isoc_sk_dskl_i21/default/table?lang=en`

\vspace{15pt}

\scriptsize

```{r, eu data1, warning=FALSE, message=FALSE}
# Download data
eurostat_data1 <- get_eurostat("isoc_sk_dskl_i21")

# Create labels for table indicators
eu_i.labels <- data.frame(
  indic_is = c(
    "I_DSK2_BAB",
    "I_DSK2_X"
  ),
  indic_is.label = c(
    "Individuals with basic or above basic overall digital skills (all five component indicators are at basic or above basic level)",
    "Individuals with no overall digital skills"
  )
)

# Subset data
eurostat_data1 <- eurostat_data1 %>%
  # filter data for 2023
  filter(str_sub(TIME_PERIOD, 1, 4) == "2023") %>%
  # select all types of respondents
  filter(ind_type == "IND_TOTAL") %>%
  # select respondents who used the internet in the last 3 months
  filter(unit == "PC_IND_IU3") %>%
  # filter data for countries in the EU (also keep the EU mean)
  filter(geo %in% eu_c.labels$geo | geo == "EU27_2020") %>%
  # rename EU27 label
  mutate(geo = ifelse(geo == "EU27_2020", "EU27 2020", geo)) %>%
  # select only needed columns
  dplyr::select(geo, indic_is, values) %>%
  # select only needed indicators
  filter(indic_is == "I_DSK2_BAB" | indic_is == "I_DSK2_X") %>%
  # attach indicator labels
  left_join(eu_i.labels) %>% relocate(indic_is.label, .after = indic_is) %>%
  # attach country names
  left_join(eu_c.labels) %>% relocate(geo.name, .after = geo)
```

\normalsize

### figure S30

\scriptsize

```{r, I_DSK2_BAB, fig.cap="Individuals with basic or above basic overall digital skills (all five component indicators are at basic or above basic level)", out.width="95%", fig.width=10, fig.height=9, fig.align='center'}
#
I_DSK2_BAB <- eurostat_data1 %>%
  filter(indic_is == "I_DSK2_BAB")
#
I_DSK2_BAB %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_DSK2_BAB$geo == "RO", "red", 
                  ifelse(I_DSK2_BAB$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals with basic or above basic overall digital skills \n(all five component indicators are at basic or above basic level)",
    subtitle = "Eurostat: Individuals' level of digital skills (from 2021 onwards) (code: isoc_sk_dskl_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) +
  scale_y_continuous(limits=c(0, 100)) +
  ggthemes::theme_pander() +
  theme(
    plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", 
                                margin = margin(t = 0, r = 20, b = 0, l = 10),
                                size = 15),
    axis.title.x = element_text(face = "bold", 
                                margin = margin(t = 20, r = 20, b = 0, l = 0),
                                size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), 
                size = 3.8, 
                vjust = 0.5, nudge_x = 0, hjust = -0.5,
                color = "black",
                angle = 90,
                fontface = "bold")
```

\normalsize

\newpage

### figure S31

\scriptsize

```{r, I_DSK2_X, fig.cap="Individuals with no overall digital skills", out.width="90%", fig.width=10, fig.height=9, fig.align='center'}
#
I_DSK2_X <- eurostat_data1 %>%
  filter(indic_is == "I_DSK2_X")
#
I_DSK2_X %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_DSK2_X$geo == "RO", "red", ifelse(I_DSK2_X$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals with no overall digital skills",
    subtitle = "Eurostat: Individuals' level of digital skills (from 2021 onwards) (code: isoc_sk_dskl_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 20)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10),size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0),size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

\newpage

## **6.2 Eurostat data from dataset "Individuals - internet activities"**

**`Online data code`**`: isoc_ci_ac_i`

**`link`**`: https://ec.europa.eu/eurostat/databrowser/view/isoc_ci_ac_i/default/table?lang=en`

\vspace{15pt}

\scriptsize

```{r, eu data2, warning=FALSE, message=FALSE}
# Download data
eurostat_data2 <- get_eurostat("isoc_ci_ac_i",
                               filters = list(ind_type = "IND_TOTAL", unit = "PC_IND_IU3"))

# Create labels for table indicators
eu_i.labels <- data.frame(
  indic_is = c(
    "I_IUSNET",
    "I_IHIF"
  ),
  indic_is.label = c(
    "Internet use: participating in social networks (creating user profile, posting messages or other contributions to facebook, twitter, etc.)",
    "Internet use: seeking health information"
  )
)

# Subset data
eurostat_data2 <- eurostat_data2 %>%
  # filter data for 2023
  filter(str_sub(time, 1, 4) == "2023") %>%
  # select all types of respondents
  filter(ind_type == "IND_TOTAL") %>%
  # select respondents who used the internet in the last 3 months
  filter(unit == "PC_IND_IU3") %>%
  # filter data for countries in the EU (also keep the EU mean)
  filter(geo %in% eu_c.labels$geo | geo == "EU27_2020") %>%
  # rename EU27 label
  mutate(geo = ifelse(geo == "EU27_2020", "EU27 2020", geo)) %>%
  # select only needed columns
  dplyr::select(geo, indic_is, values) %>%
  # select only needed indicators
  filter(indic_is == "I_IUSNET" | indic_is == "I_IHIF") %>%
  # attach indicator labels
  left_join(eu_i.labels) %>% relocate(indic_is.label, .after = indic_is) %>%
  # attach country names
  left_join(eu_c.labels) %>% relocate(geo.name, .after = geo)
```

\normalsize

### figure S32

\scriptsize

```{r, I_IUSNET, fig.cap="Internet use: participating in social networks (creating user profile, posting messages \nor other contributions to facebook, twitter, etc.)", out.width="95%", fig.width=10, fig.height=9, fig.align='center'}
#
I_IUSNET <- eurostat_data2 %>%
  filter(indic_is == "I_IUSNET")
#
I_IUSNET %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_IUSNET$geo == "RO", "red", ifelse(I_IUSNET$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Internet use: participating in social networks (creating user profile, posting messages \nor other contributions to facebook, twitter, etc.)",
    subtitle = "Eurostat: Individuals - internet activities (code: isoc_ci_ac_i)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 110)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

\newpage

### figure S33

\scriptsize

```{r, I_IHIF, fig.cap="Internet use: seeking health information", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_IHIF <- eurostat_data2 %>%
  filter(indic_is == "I_IHIF")
#
I_IHIF %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_IHIF$geo == "RO", "red", ifelse(I_IHIF$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Internet use: seeking health information",
    subtitle = "Eurostat: Individuals - internet activities (code: isoc_ci_ac_i)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 100)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

\newpage

## **6.3 Eurostat data from dataset "Evaluating data, information and digital content (2021 onwards)"**

**`Online data code`**`: isoc_sk_edic_i21`

**`link`**`: https://ec.europa.eu/eurostat/databrowser/view/isoc_sk_edic_i21/default/table?lang=en`

\vspace{15pt}

\scriptsize

```{r, eu data3, warning=FALSE, message=FALSE}
# Download data
eurostat_data3 <- get_eurostat("isoc_sk_edic_i21")

# Create labels for table indicators
eu_i.labels <- data.frame(
  indic_is = c(
    "I_UDI",
    "I_TIC",
    "I_TICCSFOI",
    "I_TICIDIS",
    "I_TICNIDIS",
    "I_TIC2",
    "I_TICXND",
    "I_TICXSKL",
    "I_TICXOTH"
  ),
  indic_is.label = c(
  "Individuals have seen untrue or doubtful information or content on the internet news sites or social media (3 months)",
  "Individuals have checked the truthfulness of the information or content they found on the internet news sites or social media (3 months)",
  "Individuals have checked the truthfulness of the information or content found on the internet by checking the sources or finding other information on the internet (3 months) (Method 1)",
  "Individuals have checked the truthfulness of the information or content found on the internet by following or taking part in discussion on internet  regarding the information (3 months) (Method 2)",
  "Individuals have checked the truthfulness of the information or content found on the internet by discussing the information offline with other persons or using sources not on internet (3 months) (Method 3)",
  "Individuals have checked the truthfulness of the information or content found on the internet (3 months) by using methods 1, 2, or 3",
  "Individuals have not checked the truthfulness of the information or content found on the internet because the individuals already knew that information, content or source was not reliable (3 months)",
  "Individuals have not checked the truthfulness of the information or content found on the internet because the individuals lacked skills or knowledge (3 months)",
  "Individuals have not checked the truthfulness of the information or content found on the internet because of other reasons (3 months)"
  )
)

# Subset data
eurostat_data3 <- eurostat_data3 %>%
  # filter data for 2023
  filter(str_sub(TIME_PERIOD, 1, 4) == "2023") %>%
  # select all types of respondents
  filter(ind_type == "IND_TOTAL") %>%
  # select respondents who used the internet in the last 3 months
  filter(unit == "PC_IND_IU3") %>%
  # filter data for countries in the EU (also keep the EU mean)
  filter(geo %in% eu_c.labels$geo | geo == "EU27_2020") %>%
  # rename EU27 label
  mutate(geo = ifelse(geo == "EU27_2020", "EU27 2020", geo)) %>%
  # select only needed columns
  dplyr::select(geo, indic_is, values) %>%
  # attach indicator labels
  left_join(eu_i.labels) %>% relocate(indic_is.label, .after = indic_is) %>%
  # attach country names
  left_join(eu_c.labels) %>% relocate(geo.name, .after = geo)
```

\normalsize

### figure S34

\scriptsize

```{r, I_UDI, fig.cap="Individuals have seen untrue or doubtful information or content on the internet news sites or social media", out.width="95%", fig.width=10, fig.height=9, fig.align='center'}
#
I_UDI <- eurostat_data3 %>%
  filter(indic_is == "I_UDI")
#
I_UDI %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_UDI$geo == "RO", "red", ifelse(I_UDI$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have seen untrue or doubtful information \nor content on the internet news sites or social media",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 100)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8,  vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

\newpage

### figure S35

\scriptsize

```{r, I_TIC, fig.cap="Individuals have checked the truthfulness of the information \nor content they found on the internet news sites or social media", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_TIC <- eurostat_data3 %>%
  filter(indic_is == "I_TIC")
#
I_TIC %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_TIC$geo == "RO", "red", ifelse(I_TIC$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have checked the truthfulness of the information \nor content they found on the internet news sites or social media",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 80)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

### figure S36

\scriptsize

```{r, I_TICCSFOI, fig.cap="Individuals have checked the truthfulness of the information or content found on the internet by checking the sources or finding other information on the internet (Method 1)", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_TICCSFOI <- eurostat_data3 %>%
  filter(indic_is == "I_TICCSFOI")
#
I_TICCSFOI %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_TICCSFOI$geo == "RO", "red", ifelse(I_TICCSFOI$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have checked the truthfulness of the information \nor content found on the internet by checking the sources or \nfinding other information on the internet (Method 1)",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 80)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

### figure S37

\scriptsize

```{r, I_TICIDIS, fig.cap="Individuals have checked the truthfulness of the information or content found on the internet by following or taking part in discussion on internet regarding the information (Method 2)", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_TICIDIS <- eurostat_data3 %>%
  filter(indic_is == "I_TICIDIS")
#
I_TICIDIS %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_TICIDIS$geo == "RO", "red", ifelse(I_TICIDIS$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have checked the truthfulness of the information \nor content found on the internet by following or taking part in \ndiscussion on internet regarding the information (Method 2)",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 40)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

### figure S38

\scriptsize

```{r, I_TICNIDIS, fig.cap="Individuals have checked the truthfulness of the information or content found on the internet by discussing the information offline with other persons or using sources not on internet (Method 3)", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_TICNIDIS <- eurostat_data3 %>%
  filter(indic_is == "I_TICNIDIS")
#
I_TICNIDIS %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_TICNIDIS$geo == "RO", "red", ifelse(I_TICNIDIS$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have checked the truthfulness of the information or \ncontent found on the internet by discussing the information offline \nwith other persons or using sources not on internet (Method 3)",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 40)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

### figure S39

\scriptsize

```{r, I_TIC2, fig.cap="Individuals have checked the truthfulness of the information or content found on the internet by using methods 1, 2, or 3", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_TIC2 <- eurostat_data3 %>%
  filter(indic_is == "I_TIC2")
#
I_TIC2 %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_TIC2$geo == "RO", "red", ifelse(I_TIC2$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have checked the truthfulness of the information \nor content found on the internet by using methods 1, 2, or 3",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 80)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

### figure S40

\scriptsize

```{r, I_TICXND, fig.cap="Individuals have not checked the truthfulness of the information or content found on the internet because the individuals already knew that information, content or source was not reliable", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_TICXND <- eurostat_data3 %>%
  filter(indic_is == "I_TICXND")
#
I_TICXND %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_TICXND$geo == "RO", "red", ifelse(I_TICXND$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have not checked the truthfulness of the information \nor content found on the internet because the individuals already \nknew that information, content or source was not reliable",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 50)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

### figure S41

\scriptsize

```{r, I_TICXSKL, fig.cap="Individuals have not checked the truthfulness of the information or content found on the internet because the individuals lacked skills or knowledge", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_TICXSKL <- eurostat_data3 %>%
  filter(indic_is == "I_TICXSKL")
#
I_TICXSKL %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_TICXSKL$geo == "RO", "red", ifelse(I_TICXSKL$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have not checked the truthfulness of the information \nor content found on the internet because the individuals lacked \nskills or knowledge",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 30)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize

### figure S42

\scriptsize

```{r, I_TICXOTH, fig.cap="Individuals have not checked the truthfulness of the information or content found on the internet because of other reasons", out.width="85%", fig.width=10, fig.height=9, fig.align='center'}
#
I_TICXOTH <- eurostat_data3 %>%
  filter(indic_is == "I_TICXOTH")
#
I_TICXOTH %>%
  ggplot() +
  aes(x = reorder(geo, values), y = values) + 
  geom_col(fill = ifelse(I_TICXOTH$geo == "RO", "red", ifelse(I_TICXOTH$geo ==  "EU27 2020", "orange", "#0C4C8A"))) +
  labs(
    # title = "Individuals have not checked the truthfulness of the information \nor content found on the internet because of other reasons",
    subtitle = "Eurostat: Evaluating data, information and digital content (2021 onwards) (code: isoc_sk_edic_i21)",
    x = "Country",
    y = "% of individuals who used the internet \nin the last 3 months (2023 data)"
  ) + scale_y_continuous(limits=c(0, 40)) +
  ggthemes::theme_pander() +
  theme(plot.title = element_text(size = 20),
    plot.subtitle = element_text(size = 9, colour = "gray40", face = "bold"),
    axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 10), size = 15),
    axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 20, b = 0, l = 0), size = 15),
    axis.text.x = element_text(face = "bold", size = 10, angle = 90, colour = "gray40"),
    axis.text.y = element_text(face = "bold", size = 10)
  )  +
  geom_text(aes(label = values), size = 3.8, vjust = 0.5, nudge_x = 0, hjust = -0.5, color = "black",
                angle = 90, fontface = "bold")
```

\normalsize
