ms.Rmd

---
title: Global well-being and mental health in the internet age
shorttitle: Global well-being
leftheader: Global well-being
author: 
  - name: Matti Vuorre
    affiliation: "1,2"
    address: Tilburg University
    email: m.j.vuorre@tilburguniversity.edu
  - name: Andrew K. Przybylski
    affiliation: "2"
    corresponding: yes
    address: Oxford Internet Institute
    email: andy.przybylski@oii.ox.ac.uk
affiliation:
  - id: 1
    institution: Tilburg School of Social and Behavioral Sciences, Tilburg University
  - id: 2
    institution: Oxford Internet Institute, University of Oxford
authornote: |
  \noindent \textbf{This pre-print is not yet peer-reviewed.}
abstract: | 
  In the last two decades the widespread adoption of internet technologies has inspired concern that they have negatively impacted mental health and psychological well-being. However, research on the topic is contested and hampered by methodological shortcomings leaving the broader consequences of internet adoption unknown. We show that the past two decades have seen only small and inconsistent changes in global well-being and mental health that are not suggestive of the idea that the adoption of internet and mobile broadband is consistently linked to negative psychological outcomes. Further investigation of this topic requires transparent study of online behaviours where they occur, namely on online platforms. We call for increased collaborative efforts between independent scientists and the internet technology sector.
  
keywords: well-being, mental health, internet technology, technology effects
wordcount: '`r scales::number(wordcountaddin:::word_count("ms.Rmd"), big.mark = ",")`'
bibliography: "references.bib"
floatsintext: yes
linenumbers: no
draft: no
mask: no
figurelist: no
tablelist: no
footnotelist: no
csl: "`r system.file('rmd', 'apa7.csl', package = 'papaja')`"
documentclass: apa7
classoption: jou
output: 
  papaja::apa6_pdf:
    number_sections: false
    keep_tex: true
  papaja::apa6_docx: 
    number_sections: false
editor_options: 
  chunk_output_type: console
---

```{r setup, include = FALSE}
# Packages
library(papaja)
library(scales)
library(tidybayes)
library(ggdist)
library(lubridate)
library(brms)
library(patchwork)
library(knitr)
library(posterior)
library(bayestestR)
library(tidyverse)

opts_chunk$set(
  eval = TRUE,
  cache = TRUE,
  include = FALSE,
  message = FALSE,
  warning = FALSE,
  error = TRUE,
  dpi = 300,
  fig.align = "center"
)

# Plotting options
theme_set(
  theme_linedraw(
    base_size = if_else(knitr::is_latex_output(), 8, 7)
  ) +
    theme(
      panel.grid = element_blank(),
      strip.text = element_text(margin = margin(4, 4, 4, 4, "pt")),
      axis.text = element_text(size = rel(0.7)),
      axis.ticks = element_line(size = rel(0.2))
    )
)
dir.create("cache", FALSE)
```

```{r data}
dat <- read_rds("data/data-all.rds")
gwp <- read_rds("data/gwp.rds")
itu <- read_rds("data/itu.rds")
gbd <- read_rds("data/gbd.rds")
z <- read_rds("data/scaling-factors.rds")
coefficients <- read_rds("cache/coefficients.rds")
epred <- read_rds("cache/epred.rds")
```

```{r variables}
# For ordering outcome factor
oo <- c(
  "Life_satisfaction", 
  "Negative_experiences", 
  "Positive_experiences",
  "Anxiety", "Depression", "Selfharm"
)
```

```{r functions}
# Function to replace extreme percentages
percent2 <- function(x, accuracy = .1) {
  x <- percent(x, accuracy = accuracy)
  x <- if_else(x == "100.0%", ">99.9%", x)
  x <- if_else(x == "0.0%", "<0.1%", x)
  x
}
```

In 2005, an estimated 17% of the global population used the internet, but by 2020 this number was already 59% [@ituITUICTStatistics2021]. Accompanying the rapid spread of the internet, worries have proliferated that its broad adoption—and technologies enabled by it such as online games, smartphones, and social media—is actively harming its users, particularly adolescents [@carrShallowsHowInternet2010; @turkleAloneTogetherWhy2011]. As a response, users and a growing number of national governments have acted in limiting access to online technologies [@buckleyChinaTightensLimits2021; @departmentfordigitalculturemediasportDraftOnlineSafety2021; @kattulaGamePoliciesContexts2021]. 

However, evidence for widespread harms of online technologies is limited. Initial reports of internet-facilitated harms have been challenged by later work informed by many methodologies including longitudinal [@jensenYoungAdolescentsDigital2019], specification curve analyses [@orbenAssociationAdolescentWellbeing2019], meta-analyses [@appelAreSocialMedia2020], and systematic reviews [@bestOnlineCommunicationSocial2014; @dickson2018screen; @odgersAnnualResearchReview2020; @ophirNewMediaScreenTime2020]. Reviews indicate early research has been hampered by inaccurate measurements of engagement with internet and related technologies [@davidsonFuzzyConstructsTechnology2022; @parrySystematicReviewMetaanalysis2021; @scharkowAccuracySelfReportedInternet2016], biased convenience samples drawn predominantly from countries in the Global North [@ghaiLackSampleDiversity2023; @ghaiSocialMediaAdolescent2022], studying a limited range of well-being outcomes [@orbenSocialMediaEnduring2019], and reliance on self-reported evaluations in place of clinical estimates of important mental health outcomes [@campbellInternetUseSocially2006]. A comprehensive test of the overall association between internet adoption and well-being and mental health, broadly defined, has therefore not been conducted. As a consequence, much of the evidence purporting to show that internet technology adoption is associated with broad changes in well-being and mental health, or even causes negative outcomes, remains equivocal [@dickson2018screen; @hawkesCMOReportUnable2019].

In this research we present two studies of global well-being and mental health in the internet age. In the first, we focus on three aspects of psychological well-being, and contrast them with yearly per-capita internet users and mobile broadband subscriptions, across 168 countries and 16 years. In the second we examine three mental health outcomes across 202 countries and 19 years. Our aim is to better understand: 1) How well-being and mental health have changed, on a global scale, during the past two decades of dramatic proliferation of internet technologies and connectivity; 2) How per-capita internet users and mobile broadband subscriptions predict country-level well-being and mental health within a given country and across countries; and 3) Assess the extent to which associations between internet technology adoption and well-being and mental health differ across age and sex and if they are specific to previously suggested vulnerable populations, such as young women.

# Study 1: Well-being

```{r data-summary-wb}
gwp_sum <- gwp %>% 
  group_by(outcome) %>% 
  summarise(n = sum(n), cells = n()) %>% 
  mutate(across(n:cells, ~number(., big.mark = ",")))
gwp_sum
```

In Study 1, we focused on psychological well-being, as reflected in self-reports of life satisfaction, positive experiences, and negative experiences from `r gwp_sum$n[1]` individuals aged 15 to 89 across `r length(unique(gwp$country))` countries from `r paste0(range(gwp$year), collapse = " to ")`. We contrasted those observations with time series data of the countries' per capita internet users and mobile broadband subscriptions [@ituITUICTStatistics2021] to examine if and how internet and mobile broadband adoption predicted psychological well-being over the past two decades.

## Transparency and Openness

This study was not preregistered. All code used in this project, along with the GBD data, link to the ITU data, and a synthetic version of the GWP data are available at <https://doi.org/10.5281/zenodo.7004053>. Our sample size was determined by the number of available years, countries, and demographics in each respective dataset.

## Methods

### Internet adoption

```{r data-summary-itu}
itu %>% 
  summarise(
    i = sum(!is.na(internet)),
    m = sum(!is.na(mobile))
  )
```

We identified the International Telecommunication Union's (ITU) database of information and communications technology (ICT) as the most comprehensive country-level source of time-series data on internet adoption [@ituITUICTStatistics2021]. The ITU has collated, from 222 countries' statistical and telecommunications agencies, the yearly percentages of population using the internet from 2000 to 2021, and yearly per capita mobile broadband subscriptions from 2007 to 2021. Our analyses used 4,200 yearly internet percentages and 2,466 broadband subscription rates.

### Subjective well-being

We examined subjective well-being indicators from the Gallup World Poll (GWP), a nationally representative annual survey of 1,000 civilian, non-institutionalised individuals aged 15 years or older from `r length(unique(gwp$country))` countries from 2005 to 2022. The surveys are conducted face-to-face or via telephone, in the respondents' native language and by local interviewers. For details on the GWP sampling and survey methodology, see [@gallupGallupWorldPoll2020; @gallupHowDoesGallup2014]. 

GWP measures subjective well-being with the positive and negative experience indices, which measure respondents' experienced well-being on the day before the survey with five items each. For positive experiences, these items are "Did you feel well-rested yesterday?"; "Were you treated with respect all day yesterday?"; "Did you smile or laugh a lot yesterday?"; "Did you learn or do something interesting yesterday?"; and "(Did you experience the following feelings during a lot of the day yesterday?) How about enjoyment?" And for negative experiences, the items are responses to "Did you experience the following feelings during a lot of the day yesterday?" with prompts "How about physical pain?"; "How about worry?"; "How about sadness?"; "How about stress?"; and "How about anger?" We aggregated both scales for each respondent by taking a mean of the five items. 

Life satisfaction in the moment was measured with one 11-step Likert item, "Please imagine a ladder, with steps numbered from 0 at the bottom to 10 at the top. The top of the ladder represents the best possible life for you and the bottom of the ladder represents the worst possible life for you. On which step of the ladder would you say you personally feel you stand at this time?", similar to the Cantril self-anchoring scale [@cantril1965pattern; @kapteynDimensionsSubjectiveWellbeing2015]. For analyses, we converted these variables to percentages, and aggregated the `r gwp_sum[1, 2]` individuals' data to means and standard errors for each outcome, country, year, sex, and age (5-year age groups from 15 to 89 years) combination (`r gwp_sum[1, 3]` cells).

### Data analysis

We analysed the data with meta-analytic bayesian hierarchical regression models [@burknerBrmsPackageBayesian2017; @gelmanDataAnalysisUsing2007; @standevelopmentteamStanModelingLanguage2022]. Bayesian methods are especially suitable to our reserch questions, because they involve many potential contrasts between e.g. countries and ages. Using bayesian hierarchical models allowed us to accurately quantify uncertainty at different levels of analyses while avoiding problems with multiple comparisons (see below).

We first asked how the well-being outcomes had changed over time. To answer, we modelled each outcome (e.g. life satisfaction) on population level intercepts (we centred time on 2010), contrasts of time (in continuous years), sex (Male: -1, Female: 1), and their interaction. We allowed all coefficients to vary across countries, age groups, and the age by country interactions, and treated all outcomes as normally distributed.

Before modelling the data, we aggregated the individual-level well-being values to means, standard deviations, and counts for each country-year-age-sex combination, in order to facilitate the Hamiltonian Monte Carlo computations on large data, and to harmonise the analyses and inferences between Studies 1 and 2. Importantly, we note that because we modelled unweighted demographic- and country-level data, our unit of analyses are demographic groups within countries, with no adjustment to e.g. countries' population sizes.

Our second research question asked how changes in internet technology adoption predicted well-being. To answer, we expanded the above model to include within-country centred per capita internet users (or mobile broadband subscriptions in a separate model) and its interaction with sex as predictors of the well-being outcome, and all their associated varying effects. We also included the between-country centered variable as a predictor to examine between-country associations between internet technology adoption and well-being. 

Importantly, we did not make the necessary strong assumptions required for identifying causal effects, and therefore highlight the descriptive nature of the resulting associations. For instance, we did not adjust for either time-invariant or -varying confounders, which could bias observed between-country and within-country associations, respectively. Nevertheless, we included time as a predictor to adjust for simple temporal trends. In addition, we within-country centred the internet technology adoption variables in order to isolate within-country associations from between-country associations.

More formally, we specified the model as

\begin{align*} 
y_i &\sim \text{Normal}(\mu_i, \sqrt{\sigma^{2} + \text{se}^2_i}), \\
\mu_i &= \alpha_{0} + \beta_{0\text{country}[i]} + \gamma_{0\text{age}[i]} + \delta_{0\text{age:country}[i]} + \\
&\quad\ (\alpha_{1} + \beta_{1\text{country}[i]} + \gamma_{1\text{age}[i]} + \delta_{1\text{age:country}[i]})\text{Time}_i + \\
&\quad\ (\alpha_{2} + \beta_{2\text{country}[i]} + \gamma_{2\text{age}[i]} + \delta_{2\text{age:country}[i]})\text{Sex}_i + \\
&\quad\ (\alpha_{3} + \beta_{3\text{country}[i]} + \gamma_{3\text{age}[i]} + \delta_{3\text{age:country}[i]})\text{Sex}_i \times \text{Time}_i + \\
&\quad\ (\alpha_{4} + \beta_{4\text{country}[i]} + \gamma_{4\text{age}[i]} + \delta_{4\text{age:country}[i]})\text{X}^{\text{CW}}_i + \\
&\quad\ (\alpha_{5} + \beta_{5\text{country}[i]} + \gamma_{5\text{age}[i]} + \delta_{5\text{age:country}[i]})\text{Sex}_i \times \text{X}^{\text{CW}}_i + \\
&\quad\ \alpha_6 \text{X}^{\text{CB}}_i,
\\
\pmb{\beta} &\sim \text{MVN}(\pmb{0}, \Sigma^\text{country}), \\
\pmb{\gamma} &\sim \text{MVN}(\pmb{0}, \Sigma^\text{age}), \\
\pmb{\delta} &\sim \text{MVN}(\pmb{0}, \Sigma^\text{age:country}),
\end{align*}

where $i$ indexes rows in the data and $\text{se}_i$ are the known standard errors. $\text{X}^{\text{CW}}_i$ and $\text{X}^{\text{CB}}_i$ are the within- and between-country centred internet technology adoption predictor (per capita internet users or mobile broadband subscriptions, in separate models), respectively. We used default noninformative priors [@burknerBrmsPackageBayesian2017] on all parameters except the following population-level coefficients, that we modelled as

\begin{align*}
\alpha_{1,3-6} &\sim \text{Normal}(0, 1), \\
\alpha_{2} &\sim \text{Normal}(0, 3).
\end{align*}

We used these priors to help model convergence, and on the assumption that differences in well-being between males and females are likely greater than differences associated with years, one-percent increases in technology adoption predictors, or the interaction terms.

Our third research question asked whether any trends in well-being or its associations with internet or mobile broadband adoption were specific to adolescents or other demographic groups. To answer, the model allowed all parameters to vary randomly across the age groups and the age by country interaction. Therefore, each age group, on average and within each country, received their own partially pooled estimates [@gelmanDataAnalysisUsing2007]. This bayesian approach to estimating age-specific associations is beneficial, given that the large number of age groups, especially within countries, would otherwise present difficulties with uncertain estimates and multiple comparisons. In addition, bayesian methods allowed quantifying uncertainties at different levels of analysis, rather than providing e.g. only point estimates of country-specific quantities. Therefore, we could compare younger age groups to older age groups, and investigate other contrasts, with confidence and without additional post-hoc adjustment procedures [@gelmanWhyWeUsually2012]. 

We conducted all data analyses with R [@rcoreteamLanguageEnvironmentStatistical2022] and estimated the models using Stan's Hamiltonian Monte Carlo sampling via the brms R package [@burknerBrmsPackageBayesian2017; @standevelopmentteamStanModelingLanguage2022]. 

Our main descriptive focus was on raw regression coefficients, which indicate contrasts in the outcomes on the percentage scale [@cohenProblemUnitsCircumstance1999] as a function of one unit (year / percent) increase in the predictor. We report these parameters' 95% Credibility Intervals (CI) and posterior probabilities of direction ($p^d$). The latter indicate the degree of certainty in the observed direction of a coefficient, and are sometimes numerically similar to one minus a one-sided *p*-value. 

In addition, on a reviewer's request we report standardized coefficients to describe the magnitudes in an alternative context. These standardized coefficients indicate contrasts in the z-scored outcomes (life satisfaction SD: `r number(z$val_sd[1], .1)`; negative experiences: `r number(z$val_sd[4], .1)`; positive experiences: `r number(z$val_sd[7], .1)`) as a function of the average year-on-year changes in internet (`r number(z$xdev[2], .1)`%) or mobile broadband (`r number(z$xdev[3], .1)`%) adoption. To determine the practical significance of the standardized associations, we report posterior probabilities that their magnitudes exceed the region of practical equivalence to zero ($p^{\text{ROPE}}$) [@kruschkeRejectingAcceptingParameter2018] as defined by [-0.1, 0.1] [@fergusonProvidingLowerboundEstimate2021]. 

## Results

Figure \@ref(fig:fig-means-wb)A describes the levels of internet and mobile broadband adoption, and the three well-being variables over time averaged over sex and age. Figure \@ref(fig:fig-means-wb)A shows the universal global penetration of internet technologies in the past two decades, and suggests that in contrast, changes in well-being are likely to be small.

```{r tbl-avg-wb}
#| include: true

tbl_average <- coefficients %>% 
  filter(
    level == "average",
    sex == "average"
  ) %>% 
  select(
    Predictor = predictor, 
    Outcome = outcome, 
    `Raw estimate` = Raw, 
    `Standardized estimate` = Scaled
  )

tbl_average %>% 
  filter(Outcome %in% oo[1:3]) %>% 
  mutate(
    Predictor = if_else(Outcome %in% oo[c(1, 4)], Predictor, ''),
    Outcome = str_replace(Outcome, "_", " ")
  ) %>% 
  # Justify considering negative estimates
  mutate(
    across(
      3:4, 
      ~if_else(
        str_starts(., "-"), 
        ., 
        paste0("\\phantom{-}", .)
      )
    )
  ) %>% 
  apa_table(
    caption = "Well-being: Population-level estimates",
    span_text_columns = TRUE, 
    font_size = "small",
    escape = FALSE,
    note = "Raw estimates are contrasts in outcome percentages associated with a one year or percentage increase in the predictor. Numbers indicate posterior means, [95\\%CIs], and ($p^d$; posterior probabilities of direction). Standardized estimates refer to contrasts in the z-scored outcome as function of the average year-on-year change in the predictor. The latter percentages indicate posterior probabilities inside a region of practical equivalence to zero ([-0.1, 0.1]; $p^{\\text{ROPE}}$)."
  )
```

### Trends

We found that, for the average country, life satisfaction had remained relatively stable, whereas both negative and positive experiences had increased (Table \@ref(tab:tbl-avg-wb); blue lines and ribbons in Figure \@ref(fig:fig-means-wb)A). Interestingly, the increase in negative experiences (second row in Table \@ref(tab:tbl-avg-wb)) was approximately fivefold to that in positive experiences. However, simple sign tests do not inform about the practical significance of magnitudes of changes. On the standardized scale, the coefficients were small ("Standardized estimate" in Table \@ref(tab:tbl-avg-wb)), and in all cases practically equivalent to zero with very high confidence. Because these trends were mixed across outcomes, and practically equivalent to zero in all cases, these results do not offer convincing evidence for consistent or meaningfully large global changes in well-being during this period of global internet technology adoption, as would be expected---all else being equal---if the latter had broad negative consequences. 

### Average associations

We then turned to our main question; understanding to what extent internet technology adoption predicted changes in well-being. We expanded the model of temporal changes to also predict well-being from within- and between-country centred per capita internet users or mobile broadband subscriptions (Table \@ref(tab:tbl-avg-wb), "Internet" and "Mobile"). The raw coefficients describe the extent to which a one-percent increase in within-country centered per capita internet users (or mobile broadband subscriptions) predict that country's well-being, adjusting for linear temporal trends. The standardized estimates describe the same quantities, but in terms of contrasts in the z-scored outcomes and average year-on-year percentage changes in the outcomes.

```{r fig-wb-ce}
#| cache: false
# Conditional effects plots

draw_double_axes <- function(pdata, predictor, ys = 1, ym = 0, xs = 1, outcome) {
  ylab <- str_replace(outcome, "_", " ")

  # pct_to_rate <- if_else(outcome %in% oo[1:3], 1, 1000)
  
  p <- pdata %>% 
    ggplot(aes(effect1__, estimate__)) +
    geom_ribbon(
      aes(ymin = lower__, ymax = upper__),
      alpha = .2, fill = "dodgerblue1"
    ) +
    geom_line(
      color = "dodgerblue3",
      linewidth = 0.75
    ) + 
    scale_y_continuous(
      ylab,
      expand = expansion(mult = 1),
      breaks = extended_breaks(),
      sec.axis = sec_axis(
        trans = ~ (. - ym) / ys, 
        breaks = extended_breaks(7),
        name = str_glue("{ylab}\n(standardized)")
        )
    ) +
    scale_x_continuous(
      predictor,
      breaks = extended_breaks(7),
      minor_breaks = seq(-50, 50, by = 5),
      sec.axis = sec_axis(
        trans = ~ . / xs, 
        breaks = extended_breaks(7),
        name = str_glue("{predictor} (scaled)")
      )
    )
  if (predictor == "Year") {
    p <- p + 
      scale_x_continuous(
        "Year",
        minor_breaks = seq(-5, 12, by = 1),
        labels = ~ . + 2010
      )
  }
  if (predictor %in% c("Internet", "Mobile")) {
    p <- p %+% filter(p$data, between(effect1__, -30, 30))
  }
  
  p
}
epred <- epred %>% 
  mutate(
    p = pmap(
      list(pdata, predictor, val_sd, val_mean, xdev, outcome),
      ~draw_double_axes(..1, ..2, ..3, ..4, ..5, ..6))
  )
```

```{r fig-means-wb}
#| include: true
#| fig.height: 6
#| fig.width: 8
#| fig.env: "figure*"
#| fig.cap: \emph{A}. Time courses of per capita internet users, mobile broadband subscriptions, and three psychological well-being outcomes. Thin dark lines indicate countries' yearly means, aggregated across sex and age. Blue lines and ribbons indicate model-implied regression lines (exploratory generalized additive model fits for Internet and Mobile broadband adoption.) \emph{B-D}. Estimated conditional means of three well-being metrics on within-country centered per-capita internet users. Primary x- (bottom) and y-axes (left) display raw predictor and outcome values, secondary axes (top, right) display z-scored outcomes and average year-on-year scaled within-country internet adoption. \emph{E-G}. Same as B-D but for mobile broadband adoption.

p1_data <- itu %>%
  select(
    country, year,
    Internet = internet,
    Mobile = mobile,
  ) %>%
  pivot_longer(
    c(Internet, Mobile),
    names_to = "outcome", values_to = "val"
  ) %>% 
  bind_rows(gwp, gbd) %>%
  group_by(country, year, outcome) %>%
  summarise(mean = mean(val, na.rm = TRUE)) %>%
  ungroup() %>% 
  mutate(outcome = factor(outcome, levels = c("Internet", "Mobile", oo))) %>% 
  drop_na(mean)

p1_data_model <- epred %>% 
  filter(predictor == "Year") %>% 
  select(predictor, outcome, pdata) %>% 
  unnest(pdata)

fig_ce_wb <- wrap_plots(
  wrap_plots(epred$p[c(2, 5, 8)]),
  wrap_plots(epred$p[c(3, 6, 9)]),
  ncol = 1
) & theme(
  panel.grid.major = element_line(linewidth = .1, color = "grey50")
)

p1 <- p1_data %>% 
  filter(!(outcome %in% oo[4:6])) %>% 
  ggplot(aes(year, mean)) +
  scale_y_continuous(
    "Value (%)",
    # Display percentages
    breaks = pretty_breaks(),
    expand = expansion(.01)
  ) +
  scale_x_continuous(
    "Year",
    breaks = c(2000, 2005, 2010, 2015, 2020),
    labels = c("'00", "'05", "'10", "'15", "'20"),
    expand = expansion(.02)
  ) +
  coord_cartesian(ylim = c(0, 100)) +
  geom_smooth(
    data = filter(p1_data, outcome %in% c("Internet", "Mobile")),
    method = "gam",
    col = "dodgerblue3",
    fill = alpha("dodgerblue1", .3),
    linewidth = 0.7
  ) +
  geom_lineribbon(
    data = p1_data_model %>% filter(outcome %in% oo[1:3]),
    aes(y = estimate__, ymin = lower__, ymax = upper__, x = year + 2010),
    col = "dodgerblue3",
    fill = alpha("dodgerblue1", .3),
    linewidth = 0.7
  ) +
  geom_line(
    aes(
      group = interaction(country, outcome)
    ),
    linewidth = .02, alpha = .2
  ) +
  facet_wrap(
    "outcome",
    scales = "free_x",
    nrow = 1,
    labeller = as_labeller(function(x) str_replace(x, "_", "\n"))
  ) +
  theme(
    legend.position = "none"
  )
(p1 / fig_ce_wb) +
  plot_layout(heights = c(22, 78)) +
  plot_annotation(tag_levels = "A")
```

Per capita internet users did not credibly predict any of the three well-being metrics: A one percent increase in per capita internet users predicted small increases in life satisfaction (Figure \@ref(fig:fig-means-wb)B), negative experiences (Figure \@ref(fig:fig-means-wb)C), and positive experiences (Figure \@ref(fig:fig-means-wb)D) for the average country, but the probabilities of direction ($p^d$) did not exceed the 95% threshold. Thus, based on simple sign tests, per-capita internet adoption was not a credible predictor of well-being. Moreover, in standardized terms, all three associations were squarely within the [-0.1, 0.1] ROPE: An average year-on-year increase in internet adoption predicted smaller than [-0.1, 0.1] changes in life satisfaction, negative and positive experiences with greater than 99.9% posterior probability. Thus, there was evidence against a practically meaningfully large association between per capita internet adoption and all three well-being metrics. 

In addition, (within country centered) per capita mobile broadband subscriptions were at best weak predictors of well-being. While per capita mobile broadband subscriptions predicted life satisfaction positively (0.06%, $p^d$ > 99.9%; Figure \@ref(fig:fig-means-wb)E), this association was negative and not credibly different from zero for negative (Figure \@ref(fig:fig-means-wb)F) or positive experiences (Figure \@ref(fig:fig-means-wb)G). In standardized terms, all three associations between per capita mobile broadband adoption and well-being were credibly equivalent to zero (Table \@ref(tab:tbl-avg-wb), "Standardized estimate").

In sum, we found that internet technology adoption did not predict life satisfaction, or negative or positive experiences to a meaningfully large degree.

### Countries

```{r tbl-cb-wb}
coefficients %>% 
  filter(
    outcome %in% oo[1:3],
    level == "cb",
    predictor != "Year"
  ) %>% 
  mutate(
    predictor = if_else(outcome == "Life_satisfaction", predictor, ''),
    outcome = str_replace(outcome, "_", " ")
  ) %>% 
  select(
    Predictor = predictor, 
    Outcome = outcome, 
    `Raw estimate` = Raw, 
    `Standardized estimate` = Scaled
  ) %>% 
  mutate(
    across(
      3:4, 
      ~if_else(
        str_starts(., "-"), 
        ., 
        paste0("\\phantom{-}", .)
      )
    )
  ) %>% 
  apa_table(
    caption = "Well-being: Between-country associations",
    span_text_columns = TRUE,
    font_size = "small",
    escape = FALSE,
    placement = "htb",
    note = "Raw estimates are contrasts in outcome percentages associated with a one percentage increase in the predictor. Numbers indicate posterior means, [95\\%CIs], and (posterior probabilities of direction). Standardized estimates refer to contrasts in the z-scored outcome as function of the average year-on-year change in the predictor. The latter percentages indicate posterior probabilities inside a region of practical equivalence to zero ([-0.1, 0.1])."
  )
```

Our model also assessed the extent to which internet technology adoption predicted well-being between countries. These between-country associations indicated, with high confidence, that countries with greater average levels of internet and mobile broadband adoption tended to report greater average levels of life satisfaction and positive experiences, and lower average levels of negative experiences (Table \@ref(tab:tbl-cb-wb)). However, based on standardized magnitudes and ROPE tests, only the positive association between mobile broadband adoption and life satisfaction was credibly different from zero with >90% confidence. Other between-country associations were equivalent to zero on the standardized scale. To be clear, between-country associations are likely to reflect multiple other causes, such as socioeconomic factors and levels of inequality, impinging on both well-being and internet technology adoption. For example while wealthier nations are likely to have greater internet coverage, they also tend to have greater healthcare coverage, which in turn is likely to affect levels of well-being. 

```{r fig-country-wb}
#| cache: true
#| include: true
#| fig.height: 6.4
#| fig.width: 8
#| fig.env: "figure*"
#| fig.cap: "Country-specific yearly changes in three well-being outcomes (top row), and their associations with per capita internet users (second row) and mobile broadband subscriptions (third row). Points and lines indicate individual countries' posterior means and 95\\% CIs. Countries are sorted from the most negative to most positive estimate. Filled points indicate that the parameter's 95\\% credibility interval excludes zero (in both panels). Raw associations (bottom x-axes) indicate changes in the outcome percentage as function of percentage (or one year) change in the predictor. Scaled associations (top x-axes) indicate changes in standardized outcome as function of the average year-on-year change in the predictor. "

figure_country <- function(y = "Life_satisfaction", x = "Year") {
  
  d <- coefficients %>% 
    filter(outcome == y, predictor == x, level == "country", sex == "average") %>% 
    left_join(z)
  
  p <- d %>% 
    # Arrange per panel
    group_by(predictor, outcome) %>% 
    arrange(predictor, outcome, mean) %>% 
    mutate(y = 1:n()) %>% 
    ggplot(aes(mean, y, col = pd > .975)) +
    geom_vline(xintercept = 0, linewidth = .33, linetype = "dashed", col = "grey70") +
    scale_color_manual(values = c("dodgerblue1", "dodgerblue3")) +
    scale_y_discrete(
      "Country",
      expand = expansion(.015)
    ) +
    scale_x_continuous(
      "Raw association",
      breaks = extended_breaks(6), 
      expand = expansion(c(0.01)),
      sec.axis = sec_axis(
        trans = ~ . * unique(d$zfactor), 
        breaks = extended_breaks(7),
        name = "Scaled association"
      )
    ) +
    scale_shape_manual(values = c(21, 19)) +
    scale_linetype_manual(values = c(1, 1, 1)) +
    geom_pointinterval(
      aes(xmin = q2.5, xmax = q97.5, shape = pd > 0.975),
      fill = "white",
      interval_size_range = c(.1, .1),
      fatten_point = 8,
      stroke = .2
    ) +
    guides(
      linetype = "none", 
      shape = "none",
      color = guide_legend(override.aes = list(size = .5, stroke = 1))
    ) +
    theme(
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      panel.grid.major.x = element_line(
        color = "black", linetype = "dotted", size = .1
      ),
      legend.position = "none"
    )
  
  if (y %in% oo[c(3, 6)]) {
    p <- p +
      facet_grid(
        rows = vars(predictor), scales = "free_y",
        labeller = as_labeller(~str_replace(., "_", " "))
      ) +
      theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
        )
  }
  
  if (y %in% oo[c(2, 5)]) {
    p <- p +
      theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
        )
  }
  if (x == "Year") {
    p <- p +
      facet_grid(
        cols = vars(outcome), scales = "free_y",
        labeller = as_labeller(~str_replace(., "_", " "))
      )
  }
  if (y %in% oo[c(3, 6)] & x == "Year") {
    p <- p +
      facet_grid(
        rows = vars(predictor),
        cols = vars(outcome),
        scales = "free_y",
        labeller = as_labeller(~str_replace(., "_", " "))
      )
  }
  
  p
}
tmp <- distinct(coefficients, outcome, predictor) %>% 
  distinct(outcome, predictor) %>% 
  filter(outcome %in% oo[1:3]) %>% 
  arrange(predictor, outcome) %>% 
  mutate(p = map2(outcome, predictor, ~figure_country(.x, .y))) %>% 
  pull(p)
wrap_plots(tmp)
```

```{r}
tmp <- coefficients %>% 
  filter(
    predictor == "Year", 
    outcome == "Life_satisfaction", 
    level == "country", 
    sex == "average"
  ) %>% 
  arrange(mean)
tmp$Raw[[1]]
tmp$Raw[[nrow(tmp)]]
```

While these trends and associations are meaningful summaries of the average country among the 168 examined, or associations between them, they do not represent specific countries which might have unique patterns of associations. Figure \@ref(fig:fig-country-wb) describes the country-level trends over time (top row) and associations linking per capita internet users (middle row) and mobile broadband subscriptions (bottom row) to each of the three outcomes (columns). We observed that all three were highly variable across countries.

First, there appeared to be no universal temporal trend that applied to all countries. Instead, the country-specific estimates varied from substantial decreases (greatest decrease in e.g. life satisfaction: `r tmp$Raw[[1]]`) to large increases in well-being (greatest increase in e.g. life satisfaction: `r tmp$Raw[[nrow(tmp)]]`). We observed similar variability in country-specific associations linking well-being to internet adoption and to mobile broadband adoption (Figure \@ref(fig:fig-country-wb)).

```{r tbl-country-prop-wb}
tbl_props <- coefficients %>% 
  filter(level == "country", sex == "average") %>% 
  mutate(
    Sign = case_when(
      pd > .975 & sign(mean)==-1 ~ "Negative",
      pd > .975 & sign(mean)==1 ~ "Positive",
      TRUE ~ "Indecisive"
    ),
    ROPE = case_when(
      rope < .05 & sign(mean)==-1 ~ "Negative",
      rope < .05 & sign(mean)==1 ~ "Positive",
      rope > .95 ~ "Null",
      TRUE ~ "Indecisive"
    )
  ) %>% 
  select(predictor, outcome, group, Sign, ROPE) %>% 
  pivot_longer(c(Sign, ROPE)) %>% 
  count(predictor, outcome, name, value) %>% 
  mutate(
    n = str_glue("{n} ({percent(n/sum(n), 1)})"), 
    .by = c(outcome, predictor, name)
  ) %>% 
  pivot_wider(names_from = value, values_from = n) %>% 
  mutate(across(c(Indecisive, Negative, Null, Positive), ~replace_na(., "0 (0%)"))) %>% 
  rename(Predictor = predictor, Outcome = outcome, Test = name) %>% 
  arrange(Predictor, Outcome, desc(Test)) %>%
  pivot_wider(
    names_from = Test, 
    values_from = c(Indecisive, Negative, Null, Positive), 
    names_vary = "slowest"
  ) %>% 
  select(-Null_Sign)

tbl_props %>% 
  filter(Outcome %in% oo[1:3]) %>% 
    mutate(
    Outcome = fct_inorder(str_replace(Outcome, "_", " "))
  ) %>% 
  as.data.frame() %>% 
  setNames(str_remove(names(.), "_ROPE") %>% str_remove("_Sign")) %>% 
  apa_table(
    placement = "htb",
    span_text_columns = TRUE,
    col_spanners = list(` ` = c(1, 2), `Sign test` = c(3, 5), `ROPE` = c(6, 9)),
    font_size = "small",
    caption = "Well-being: Country-specific test summary",
    note = "Numbers indicate counts and (percentages) of countries with different decisions either based on a sign test or a ROPE of [-0.1, 0.1] on the standardized scale."
  )
```

To summarize these country-specific associations, we calculated proportions of countries with credibly positive, negative, null, and inconclusive associations based on either simple sign tests on the raw scale ($p^d$ > 97.5%), or ROPE tests on the standardized scale ($p^{\text{ROPE}}$ > 95%). Table \@ref(tab:tbl-country-prop-wb) shows that the association between internet adoption and life satisfaction was credibly positive, based on sign tests, for `r tbl_props$Positive_Sign[7]` countries, negative for `r tbl_props$Negative_Sign[7]` countries, and indecisive for `r tbl_props$Indecisive_Sign[7]` countries. Evaluated against a region of practical equivalence (to zero) of [-0.1, 0.1], `r tbl_props$Null_ROPE[7]` countries showed credibly null associations between internet adoption and life satisfaction, whereas the result was indecisive for `r tbl_props$Indecisive_ROPE[7]` countries, positive for `r tbl_props$Positive_ROPE[7]` and negative for `r tbl_props$Negative_ROPE[7]` countries. The country-specific associations between negative and positive experiences and per capita internet adoption were broadly similar to those with life satisfaction: When evaluated against a ROPE on the standardized scale, the majority of associations were either indecisive or credibly null. Moreover, the results regarding country-specific associations between mobile broadband adoption and well-being were broadly similar: The majority of associations were either indecisive or credibly null.

This heterogeneity and lack of consistent associations across countries should qualify any inferences concerning associations for the average country, and is further evidence against the idea that the adoption of internet or mobile broadband has had uniform global negative effects on well-being. 

### Demographics

```{r fig-demographics-wb}
#| cache: true
#| include: true
#| fig.height: 6.4
#| fig.width: 8
#| fig.env: "figure*"
#| fig.cap: "Age- and sex-specific (green: female, red: male) changes in three well-being outcomes (top row), and their associations with per capita internet users (middle row) and mobile broadband subscriptions (bottom row). Bottom x-axes indicate raw associations; top x-axes are scaled associations."

figure_demo <- function(x = "Year", y = "Life_satisfaction") {
  
  # Get coefficients for males and females for average and each age
  d <- coefficients %>% 
    filter(
      outcome == y, 
      predictor == x, 
      level %in% c("average", "age"), 
      sex %in% c("male", "female")
    ) %>%
    mutate(group = if_else(group == "Average", str_to_title(sex), group)) %>% 
    left_join(z) %>% 
    mutate(
      group = fct_rev(fct_relevel(group, "Female", "Male"))
    )
  
  p <- d %>% 
    ggplot(aes(mean, group, col = sex, shape = pd > 0.975)) +
    geom_vline(xintercept = 0, linewidth = .2, linetype = "dashed", col = "grey70") +
    scale_color_brewer(palette = "Set2") +
    scale_shape_manual(values = c(21, 19), breaks = c(FALSE, TRUE)) +
    geom_pointinterval(
      aes(xmin = q2.5, xmax = q97.5),
      fill = "white",
      interval_size_range = c(.4, .4),
      fatten_point = 3,
      stroke = .75,
      position = position_dodgejust()
    ) +
    guides(shape = "none") +
    scale_x_continuous(
      "Raw association",
      breaks = extended_breaks(6), 
      expand = expansion(c(0.05)),
      sec.axis = sec_axis(
        trans = ~ . / unique(d$val_sd) * unique(d$xdev), 
        breaks = extended_breaks(),
        name = "Scaled association"
      )
    ) +
    theme(
      axis.title.y = element_blank(),
      panel.grid.major = element_line(
        color = "black", linetype = "dotted", size = .1
      ),
      legend.position = "none"
    )
  
  if (y %in% oo[c(3, 6)]) {
    p <- p +
      facet_grid(
        rows = vars(predictor),
        labeller = as_labeller(~str_replace(., "_", " "))
      ) +
      theme(
        # axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
        )
  }
  
  if (y %in% oo[c(2, 5)]) {
    p <- p +
      theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
        )
  }
  
  if (x == "Year") {
    p <- p +
      facet_grid(
        cols = vars(outcome),
        labeller = as_labeller(~str_replace(., "_", " "))
      )
  }
  if (y %in% oo[c(3, 6)] & x == "Year") {
    p <- p +
      facet_grid(
        rows = vars(predictor),
        cols = vars(outcome),
        labeller = as_labeller(~str_replace(., "_", " "))
      )
  }
  
  p
}

tmp <- coefficients %>% 
  distinct(predictor, outcome) %>% 
  filter(outcome %in% oo[1:3]) %>% 
  arrange(predictor, outcome) %>% 
  mutate(p = map2(predictor, outcome, ~figure_demo(.x, .y))) %>% 
  pull(p)
wrap_plots(tmp) &
  ylab("Demographic group")
```

These global and country-specific analyses are informative but shed no light on how internet and mobile broadband adoption might be differently associated with well-being across demographic groups. So we then examined variation across age and sex (for the average country) in the relations linking internet and mobile broadband adoption to psychological well-being. Age- and sex-specific estimates, for the average country, are shown in Figure \@ref(fig:fig-demographics-wb). 

```{r}
tmp <- coefficients %>% 
  filter(sex == "female-male", outcome == "Life_satisfaction", level == "average")
```

We found that, for the average country, life satisfaction had increased more for females (`r tmp$Raw[[1]]`), but, on the standardized scale, this interaction term was very small and practically equivalent to zero with near absolute certainty (`r tmp$Scaled[[1]]`). Increases in negative and positive experiences were also greater for females, and while the sign test indicated greater increases in negative experiences for women, both interaction terms were within the ROPE with near absolute certainty. We present these results in detail in Table \@ref(tab:tbl-sex1-wb).

```{r tbl-sex1-wb}
coefficients %>% 
  filter(
    level == "average",
    sex == "female-male",
    outcome %in% oo[1:3]
  ) %>% 
  mutate(
    predictor = if_else(outcome == "Life_satisfaction", predictor, ''),
    outcome = str_replace(outcome, "_", " ")
  ) %>% 
  select(
    Predictor = predictor, 
    Outcome = outcome, 
    `Raw estimate` = Raw, 
    `Standardized estimate` = Scaled
  ) %>% 
  # Justify considering negative estimates
  mutate(
    across(
      3:4, 
      ~if_else(
        str_starts(., "-"), 
        ., 
        paste0("\\phantom{-}", .)
      )
    )
  ) %>% 
  apa_table(
    caption = "Well-being: Predictor-sex interactions",
    span_text_columns = TRUE,
    placement = "htb",
    font_size = "small",
    escape = FALSE,
    note = "Parameter estimates indicate differences between males and females."
  )
```

There were no further differences in associations linking well-being and (mobile) broadband adoption between men and women when evaluated either with sign tests (Table \@ref(tab:tbl-sex1-wb), "Raw estimate") or equivalence tests against [-0.1, 0.1] (Table \@ref(tab:tbl-sex1-wb), "Standardized estimate"). Overall, these demography-specific trends and associations indicated no clear patterns specific to a particular demographic group. Neither did they support the commonly offered narrative that young individuals, particularly young women, have experienced disproportionately large decrements in well-being in association with the adoption of internet technologies.

## Discussion

We examined the extent to which three indicators of psychological well-being had changed over time and in association with per capita internet users and mobile broadband subscriptions. Across 168 countries and 16 years, we found little to no support for changes in well-being over time and in association with internet technology adoption. Although we found that both negative and positive experiences had increased, based on simple sign tests, these changes were small enough to be credibly practically equivalent to zero when evaluated against a smallest association magnitude of 0.1 on the standardized scale.

Moreover, although mobile broadband adoption did positively predict life satisfaction, when evaluated against ROPE on a standardized scale, this association, like all others, was credibly equivalent to zero.

We then examined variation in trends and associations across countries: Variability in sign and magnitude was large, suggesting no consistent connection between internet technology adoption and country-level well-being. Finally, our analysis of differences across demographic groups did not support the idea that young individuals, or young women particularly, were the most at-risk group. Nevertheless, this analysis was necessarily constrained to a limited range of available outcomes reflecting subjective well-being [e.g. @jebbSubjectiveWellBeingWorld2020]. In Study 2, we extended our investigation to mental health outcomes.

# Study 2: Mental health

In the second study we extended our investigation to focus on mental health, using meta-analytic rates of anxiety, depression, and self-harm from 2000 to 2019 among 202 countries, and their associations with internet technology adoption, in place of self-reported measures of psychological well-being. 

## Methods

### Mental health outcomes

We studied meta-analytic prevalence rates (per 100,000 individuals) of anxiety disorders (ICD10 F40-F44.9, F93-F93.2), depressive disorders (ICD10 F32-F33.9, F34.1), and self-harm (ICD10 X60-X64.9, X66-X84.9, Y87.0) in 204 countries from 2000 to 2019 as estimated by the Institute for Health Metrics and Evaluation's (IHME) Global Burden of Disease 2019 (GBD) study [@jamesGlobalRegionalNational2018; @vosGlobalBurden3692020]. The GBD collates heterogeneous data from all WHO member states' censuses, household surveys, civil registration and vital statistics, disease registries, health service use statistics, disease notifications, and other sources. It then aggregates data from these sources with bayesian meta-regression to produce country-specific yearly prevalence estimates. 

The GBD 2019 prevalence rate estimates are based on 19,773 data sources with varying coverage for individual countries; for details of the GBD 2019 methodology, see [@vosGlobalBurden3692020] and especially Appendix 1 therein. The prevalence rates are estimated for females and males in 5-year age groups, and are provided as the IHME meta-regression model's predicted rates and 95% credibility intervals; we converted the latter to approximate standard errors for our meta-analytic modelling strategy (see Study 1, Methods). 

```{r}
gbd %>% 
  count(outcome)
```

We emphasize that the GBD estimates are not observed data, and therefore are accurate only to the extent that the GBD's data collection methods and modelling strategies are valid. We have compared the GBD estimates to the CDC's estimates of self-harm in the United States [@u.s.centersfordiseasecontrolandpreventionWISQARSFatalInjury2022], and found that they are likely to deviate in systematic ways from other authoritative information sources. We nevertheless argue that because the GBD provides the most comprehensive dataset of global mental health, studying these estimates is informative, but emphasize this caveat. The sample size for our analyses (combinations of country, year, sex, and age) was 130,560. 

### Data analysis

Before analyses, we converted the GBD meta-analytic rates (per 100k population) to percentages, for consistency with Study 1. Otherwise, we analysed the data in the same manner as in Study 1, but did not include varying parameters over the age by country interaction because the models did not converge due to invariance in the data. As above, we report coefficients on the raw scale, with 95%CIs and $p^d$, as well as coefficients scaled by the ratio of the predictor's standard deviation to the standard deviation in the outcome. We test the latter for practically large coefficients with a ROPE of [-0.1, 0.1].

## Results

```{r tbl-avg-mh}
#| include: true

tbl_average_mh <- coefficients %>% 
  filter(
    level == "average",
    sex == "average",
    outcome %in% oo[4:6]
  ) %>% 
  left_join(z) %>% 
  mutate(
    Raw = str_glue(
      "{number(mean, .0001)} [{number(q2.5, .0001)},",
      "{number(q97.5, .0001)}] ({percent2(pd, .1)})"
    ),
    Scaled = str_glue(
      "{number(zmean, .0001)} [{number(zq2.5, .0001)},",
      "{number(zq97.5, .0001)}] ({percent2(rope, .1)})"
    )
  ) %>% 
  select(
    Predictor = predictor, 
    Outcome = outcome, 
    `Raw estimate` = Raw, 
    `Standardized estimate` = Scaled
  )
tbl_average_mh %>% 
  mutate(
    Predictor = if_else(Outcome %in% oo[c(1, 4)], Predictor, ''),
    Outcome = str_replace(Outcome, "_", " ")
  ) %>% 
  mutate(
    across(
      3:4, 
      ~if_else(
        str_starts(., "-"), 
        ., 
        paste0("\\phantom{-}", .)
      )
    )
  ) %>% 
  apa_table(
    caption = "Mental health: Population-level estimates",
    span_text_columns = TRUE,
    font_size = "small",
    escape = FALSE,
    note = "Raw estimates are contrasts in outcome percentages associated with a one year or percentage increase in the predictor. Numbers indicate posterior means, [95\\%CIs], and (posterior probabilities of direction). Standardized estimates refer to contrasts in the z-scored outcome as function of the average year-on-year change in the predictor. The latter percentages indicate posterior probabilities inside a region of practical equivalence to zero ([-0.1, 0.1])."
  )
```

### Trends

Because Study 2 concerned the Global Burden of Disease's meta-analytic estimates of mental health rates, we observed less variation in the variables over time than was the case with the well-being outcomes in Study 1 (Figure \@ref(fig:fig-data-mh)A). Our results showed that rates of anxiety had increased (`r tbl_average_mh[["Raw estimate"]][[1]]`; Table \@ref(tab:tbl-avg-mh), "Raw estimate"), whereas those of depression (`r tbl_average_mh[["Raw estimate"]][[2]]`) and self-harm (`r tbl_average_mh[["Raw estimate"]][[3]]`) had decreased, and the model certainty in the directions of these average trends was high (>99.9%). However, evaluated against a ROPE of [-0.1, 0.1] on the standardized scale, all three trends were credibly equivalent to zero (Table \@ref(tab:tbl-avg-mh), "Standardized estimate"). Thus, similarly to well-being trends, although sign tests suggested credibly non-null changes in mental health over time, when evaluated against the [-0.1, 0.1] ROPE on the standardized scale, there was evidence in support of smaller than meaningfully large associations.

```{r fig-data-mh}
#| include: true
#| fig.height: 6
#| fig.width: 8
#| fig.env: "figure*"
#| fig.cap: \emph{A}. Time courses of three mental health outcomes. Thin dark lines indicate countries' yearly means, aggregated across sex and age. Blue lines and ribbons indicate model-implied regression lines. \emph{B-D}. Estimated conditional means of three mental health metrics on within-country centered per-capita internet users. Primary x- (bottom) and y-axes (left) display raw predictor and outcome values, secondary axes (top, right) display z-scored outcomes and average year-on-year scaled within-country internet adoption. \emph{E-G}. Same as B-D but for mobile broadband adoption.

fig_ce_mh <- wrap_plots(
  wrap_plots(epred$p[c(11, 14, 17)]),
  wrap_plots(epred$p[c(12, 15, 18)]),
  ncol = 1
) & theme(
  panel.grid.major = element_line(linewidth = .1, color = "grey50")
)

p1 <- p1_data %>% 
  filter(outcome %in% oo[4:6]) %>% 
  ggplot(aes(year, mean)) +
  scale_y_continuous(
    "Value (%)",
    # Display percentages
    breaks = pretty_breaks(),
    expand = expansion(.01)
  ) +
  scale_x_continuous(
    "Year",
    breaks = c(2000, 2005, 2010, 2015, 2020),
    labels = c("'00", "'05", "'10", "'15", "'20"),
    expand = expansion(.02)
  ) +
  # coord_cartesian(ylim = c(0, 100)) +
  geom_lineribbon(
    data = p1_data_model %>% filter(outcome %in% oo[4:6]),
    aes(y = estimate__, ymin = lower__, ymax = upper__, x = year + 2010),
    col = "dodgerblue3",
    fill = alpha("dodgerblue1", .25),
    linewidth = 0.6
  ) +
  geom_line(
    aes(
      group = interaction(country, outcome)
    ),
    linewidth = .025, alpha = .2
  ) +
  facet_wrap(
    "outcome",
    scales = "free",
    nrow = 1,
    labeller = as_labeller(function(x) str_replace(x, "_", "\n"))
  ) +
  theme(
    legend.position = "none"
  )
(p1 / fig_ce_mh) +
  plot_layout(heights = c(22, 78)) +
  plot_annotation(tag_levels = "A")
```

### Average associations

Then, to answer our primary research question, we examined the within-country associations between internet technology adoption and mental health. We observed no credible (at the 95% level) associations between internet adoption and either anxiety, depression, or self-harm (Figure \@ref(fig:fig-data-mh)B), evaluated either against a point null on the raw scale (anxiety: `r tbl_average_mh[["Raw estimate"]][[4]]`; depression: `r tbl_average_mh[["Raw estimate"]][[5]]`; self-harm: `r tbl_average_mh[["Raw estimate"]][[6]]`), or against a ROPE of [-0.1, 0.1] on the standardized scale (Table \@ref(tab:tbl-avg-mh), "Standardized estimate")). Similarly, per capita mobile broadband subscriptions were not credible predictors of either anxiety, depression, or self-harm (Table \@ref(tab:tbl-avg-mh)). These results suggested that, all else being equal, country-level (mobile) internet adoption does not predict mental health.

We also examined the between-country associations linking internet technology adoption with mental health. Countries with greater levels of per capita internet users, on average, tended to be those with lower levels of depression, and higher levels of self-harm (Figure \@ref(fig:fig-data-mh)B, Between country estimates). Similarly, per capita mobile broadband subscriptions predicted lower rates of both anxiety and depression, but higher rates of self-harm. However, all of these between-country associations were within the region of practical equivalence to zero with very high confidence (>99.9%).

### Countries

```{r tbl-country-prop-mh}
tbl_props_mh <- tbl_props %>% 
  filter(Outcome %in% oo[4:6]) %>% 
  mutate(
    Outcome = fct_inorder(str_replace(Outcome, "_", " "))
  )
tbl_props_mh %>% 
  as.data.frame() %>% 
  setNames(str_remove(names(.), "_ROPE") %>% str_remove("_Sign")) %>% 
  apa_table(
    placement = "htb",
    span_text_columns = TRUE,
    col_spanners = list(` ` = c(1, 2), `Sign test` = c(3, 5), `ROPE` = c(6, 9)),
    font_size = "small",
    caption = "Mental health: Country-specific test summary",
    note = "Numbers indicate counts and (percentages) of countries with different decisions either based on a sign test or a ROPE of [-0.1, 0.1] on the standardized scale."
  )
```

We then examined the heterogeneity in the associations between countries. Although there were considerable differences in changes in anxiety, depression, and self-harm between countries, the associations varied less (Figure \@ref(fig:fig-country-mh)). To summarize these country-level trends and associations, we calculated proportions of countries with negative, positive, null, and indecisive results based on both sign tests and ROPE tests, as in Study 1 (Table \@ref(tab:tbl-country-prop-mh)). First, we find that the vast majority of country-specific trends and associations were indecisive when evaluated against a point null hypothesis (Table \@ref(tab:tbl-country-prop-mh)). Moreover, when evaluated against a ROPE of [-0.1, 0.1] all country-level coefficients were credibly equivalent to zero.

Thus, analyses of country-specific associations provided stronger evidence that (mobile) internet penetration does not meaningfully predict fluctuations in countries' mental health.

```{r tbl-cb-mh}
coefficients %>% 
  filter(
    outcome %in% oo[4:6],
    level == "cb",
    predictor != "Year"
  ) %>% 
  mutate(
    predictor = if_else(outcome %in% oo[c(1, 4)], predictor, ''),
    outcome = str_replace(outcome, "_", " ")
  ) %>% 
  select(
    Predictor = predictor, 
    Outcome = outcome, 
    `Raw estimate` = Raw, 
    `Standardized estimate` = Scaled
  ) %>% 
  # Justify considering negative estimates
  mutate(
    across(
      3:4, 
      ~if_else(
        str_starts(., "-"), 
        ., 
        paste0("\\phantom{-}", .)
      )
    )
  ) %>% 
  apa_table(
    caption = "Mental health: Between-country associations",
    span_text_columns = TRUE,
    font_size = "small",
    escape = FALSE,
    placement = "htb",
    note = "Raw estimates are contrasts in outcome percentages associated with a one percentage increase in the predictor. Numbers indicate posterior means, [95\\%CIs], and (posterior probabilities of direction). Standardized estimates refer to contrasts in the z-scored outcome as function of the average year-on-year change in the predictor. The latter percentages indicate posterior probabilities inside a region of practical equivalence to zero ([-0.1, 0.1])."
  )
```

```{r fig-country-mh}
#| cache: true
#| include: true
#| fig.height: 6.4
#| fig.width: 8
#| fig.env: "figure*"
#| fig.cap: "Country-specific yearly changes in three mental health outcomes (top row), and their associations with per capita internet users (second row) and mobile broadband subscriptions (third row). Points and lines indicate individual countries' posterior means and 95\\% CIs. Countries are sorted from the most negative to most positive estimate. Filled points indicate that the parameter's 95\\% credibility interval excludes zero (in both panels). Raw associations (bottom x-axes) indicate changes in the outcome percentage as function of percentage (or one year) change in the predictor. Scaled associations (top x-axes) indicate changes in standardized outcome as function of the average year-on-year change in the predictor. "

tmp <- distinct(coefficients, outcome, predictor) %>% 
  distinct(outcome, predictor) %>% 
  filter(outcome %in% oo[4:6]) %>% 
  arrange(predictor, outcome) %>% 
  mutate(p = map2(outcome, predictor, ~figure_country(.x, .y))) %>% 
  pull(p)
wrap_plots(tmp)
```

### Demographics

```{r fig-demographics-mh}
#| cache: true
#| include: true
#| fig.height: 6.6
#| fig.width: 8
#| fig.env: "figure*"
#| fig.cap: "Age- and sex-specific (green: female, red: male) changes in three mental health outcomes (top row), and their associations with per capita internet users (middle row) and mobile broadband subscriptions (bottom row). Bottom x-axes indicate raw associations; top x-axes are scaled associations."

tmp <- coefficients %>% 
  distinct(predictor, outcome) %>% 
  filter(outcome %in% oo[4:6]) %>% 
  arrange(predictor, outcome) %>% 
  mutate(p = map2(predictor, outcome, ~figure_demo(.x, .y))) %>% 
  pull(p)
wrap_plots(tmp) &
  ylab("Demographic group")
```

We then turned to our third research question, and studied how trends in mental health, or its associations with internet technology adoption varied across demographic groups. Changes in rates of anxiety showed a complex pattern, whereby older males' rates of anxiety increased more than younger males', except the very youngest age group (Figure \@ref(fig:fig-demographics-mh)). For females, on the other hand, rates of anxiety increased least for the oldest and youngest groups. Overall, the difference in changes in anxiety between males and females was credibly null (Table \@ref(tab:tbl-sex1-mh)). Depression had decreased more for females than males (Table \@ref(tab:tbl-sex1-mh)), and for the middle age groups vs. the youngest and oldest groups. There were no average sex differences in rates of decline of self-harm (Table \@ref(tab:tbl-sex1-mh)), but the declines were greater for older than younger age groups (Figure \@ref(fig:fig-demographics-mh)).

Overall, the demography-based differences in the mental health-internet adoption associations were smaller than within the mental health trends (Figure \@ref(fig:fig-demographics-mh)). There were no sex differences in associations linking either per capita internet users or mobile broadband subscriptions to either mental health outcome that were credibly different from zero. Tested against a ROPE of [-0.1, 0.1], all sex differences were credibly null. Age-based differences in links with anxiety were small. However, depression was negatively and credibly associated with internet adoption for older age groups, but positively and credibly associated with per capita mobile broadband subscriptions for both sexes in the 10 to 14 and 15 to 19 year-old age brackets, based on sign tests. The magnitudes of these associations were very small, and supported null associations on the standardized scale. Similarly, per capita internet users predicted self-harm negatively for older age groups in the 65 to 89 year-old brackets, but positively in the 10 to 34 year-old brackets. These results indicated some support for the idea that younger individuals' mental health may be more negatively associated with the adoption of internet technologies, although this finding was not consistent across the two different predictors or three different mental health outcomes, and in fact null when considered in the context of a ROPE of [-0.1, 0.1] on the standardized scale.

```{r tbl-sex1-mh}
coefficients %>% 
  filter(
    level == "average",
    sex == "female-male",
    outcome %in% oo[4:6]
  ) %>% 
  mutate(
    predictor = if_else(outcome %in% oo[c(1, 4)], predictor, ''),
    outcome = str_replace(outcome, "_", " ")
  ) %>% 
  select(
    Predictor = predictor, 
    Outcome = outcome, 
    `Raw estimate` = Raw, 
    `Standardized estimate` = Scaled
  ) %>% 
  # Justify considering negative estimates
  mutate(
    across(
      3:4, 
      ~if_else(
        str_starts(., "-"), 
        ., 
        paste0("\\phantom{-}", .)
      )
    )
  ) %>% 
  apa_table(
    caption = "Mental health: Predictor-sex interactions",
    span_text_columns = TRUE,
    placement = "htb",
    font_size = "small",
    escape = FALSE,
    note = "Parameter estimates indicate differences between males and females."
  )
```

## Discussion

In Study 2, we replicated the analyses regarding well-being in Study 1 using three mental health outcomes from a different dataset. We found mixed results regarding changes in the rates of mental health problems, whereby rates of anxiety had increased, but rates of depression and self-harm decreased, for the average country and when tested against a point null hypothesis. However, when tested against a ROPE of [-0.1, 0.1] on a standardized scales, all associations were credibly equivalent to that null interval. Similar to the average contrasts, country-, age-, and sex-specific coefficients were all credibly equivalent to zero. With the important caveat that the Global Burden of Disease mental health rates analysed here are noisy meta-analytic estimates, our results suggested only small if any average changes in mental health globally, and moreover no associations between mental health and internet technology adoption.

# General discussion

The idea that the rapid and global penetration of the internet, and technologies enabled by it, is affecting psychological well-being and mental health is compelling but not adequately tested. In two studies, we focused on the degree to which well-being and mental health had changed in the past two decades across most world countries. Based on sign tests, we found mostly mixed and inconclusive results; based on ROPE tests and a smallest coefficient size of interest of 0.1, we found evidence against meaningfully large changes in well-being and mental health. These results align with other results that tend to be mixed or find few changes over time [@borAreChildAdolescent2014; @collishawAnnualResearchReview2015; @cosmaCrossNationalTimeTrends2020; @keyesRecentIncreasesDepressive2019]. In addition, although previous reports have suggested increasing mental health problems specific to young females [e.g. @borAreChildAdolescent2014], we did not find support for such changes for the average country.

We then measured the extent to which per capita internet users and mobile broadband subscriptions predicted levels of well-being and mental health (within countries). Although confidence for direction of these coefficients exceeded 95% in one case out of twelve, on a standardized scale all 12 associations were practically equivalent to zero with very high confidence. The small associations were also qualified by variability across countries and demographics. Furthermore, countries where internet technologies were more widely adopted overall were not those with lowest levels of well-being and mental health. 

We also examined the moderating roles of sex and age: Both depression and self-harm were more positively associated with per capita internet users and mobile broadband subscriptions in the younger age brackets---including the youngest studied here at 10 to 14 years---than they were for older age groups. However, evaluated against a smallest association of interest, all sex- and age-based associations provided evidence for null associations. Studies on the association between internet technology use and mental health across wide age ranges are scarce. Nevertheless, while we found few gender differences, our results are in line with one study reporting that cross-sectional associations between social media use and life satisfaction are more negative at specific time windows in adolescence [@orbenWindowsDevelopmentalSensitivity2022].

Our results do not provide evidence supporting the view that the internet, and technologies enabled by it such as smartphones with internet access, are actively promoting or harming either well-being or mental health globally. Our results also suggest that age- and sex-based differences in associations linking internet technology adoption and well-being are small. However, we did find some evidence suggesting that internet technology adoption is more negatively associated with young individuals' mental health than older groups' (whose associations were, in some cases, positive). Nevertheless, these associations, although sometimes credibly different from zero, appeared small in magnitude and were credibly in the null region when evaluated against a SESOI of 0.1 on the standardized scale. Because of their apparently small magnitude and the fact that the data are not at the individual level, we highlight that these contrasts do not speak to whether the associations might be clinically relevant.

It is important to note that the data and theory required to address this question at the causal level are absent. Consequently, our analyses cannot account for potential confounders in the associations linking mental health and well-being to internet technology adoption. Our descriptions, therefore, are suggestive but are not intended to provide evidence for, or against, causal relations. Although the first to study this topic within a global perspective, future studies should investigate other potential differences based on demographic, economic, and cultural factors, which we could not address here. In addition, our results regarding mental health, while comprehensive, are necessarily and considerably less certain than those on well-being, due to the lack of diverse global data and our own reliance on GBD estimates in place of objective measurements of mental health. In addition, the data aggregation necessary to link internet adoption and well-being at the national level preclude investigating other important and relevant questions, such as comorbidities.

Research on the effects of internet technologies is stalled because the data most urgently needed are collected and held behind closed doors by technology companies and online platforms. It is crucial to study, in more detail and with more transparency from all stakeholders, data on individual adoption of and engagement with internet-based technologies. These data exist and are continuously analysed by global technology firms for marketing and product improvement, but unfortunately are not accessible for independent research. It remains a fundamental challenge to this field of inquiry to ensure that this information is accessible to independent scholars. Until these data can be transparently analysed for the public good, the potential harmful effects of the internet and other digital environments will remain unknown.

# Material availability and ethics

All the data-analytic code and a synthetic GWP dataset to facilitate verification, extension, and critique of our procedures are available at <https://doi.org/10.5281/zenodo.7004053>. The ITU data, along with a complete list of years, countries, and sources of data are freely available at <https://www.itu.int/en/ITU-D/Statistics/Pages/stat/default.aspx>. The GWP data are available to subscribing institutions through Gallup. The GBD data are freely available at <http://ghdx.healthdata.org/gbd-results-tool>.

This study and methods therein were approved by the University of Oxford Central University Research Ethics Committee (SSH_OII_CIA_21_084).

# Author contributions

MV: Conceptualization, data curation, formal analysis, investigation, methodology, project administration, software, validation, visualization, writing - original draft preparation, writing - review & editing. AKP: Conceptualization, funding acquisition, investigation, methodology, project administration, resources, supervision, visualization, writing - original draft preparation, writing - review & editing.

# Conflicts of interest

The authors declare no competing interests. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

# Acknowledgements

We thank Niklas Johannes for feedback on this work. The authors would like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work http://dx.doi.org/10.5281/zenodo.22558.

# Funding

This research was supported by the Huo Family Foundation and Economic and Social Research Council Grant ES/T008709/1. 

# References

::: {#refs}
:::

\clearpage
\onecolumn

<!-- Prepend table & figure numbers with A and reset counters -->
\setcounter{page}{1}
\setcounter{table}{0}
\setcounter{figure}{0}
\renewcommand{\thetable}{A\arabic{table}}
\renewcommand{\thefigure}{A\arabic{figure}}

\clearpage
# Appendix 1

```{r tbl-cb-wb}
#| include: true
```


```{r tbl-country-prop-wb}
#| include: true
```


```{r tbl-sex1-wb}
#| include: true
```

\clearpage

```{r tbl-cb-mh}
#| include: true
```

```{r tbl-country-prop-mh}
#| include: true
```


```{r tbl-sex1-mh}
#| include: true
```