appendix.Rmd

\renewcommand{\figurename}{Supplementary Figure}
\renewcommand{\tablename}{Supplementary Table}
\renewcommand{\thefigure}{\arabic{figure}} \setcounter{figure}{0}
\renewcommand{\thetable}{\arabic{table}} \setcounter{table}{0}
\renewcommand{\theequation}{S\arabic{equation}} \setcounter{equation}{0}

# **Supplementary Information**

\setcounter{page}{1}
\centering
\noindent \small Cross-national analyses require additional controls to account for the non-independence of nations \newline
\hspace*{1cm} \small Scott Claessens\textsuperscript{1}, Thanos Kyritsis\textsuperscript{1}, & Quentin D. Atkinson\textsuperscript{1,2} \newline
\raggedright
\noindent \footnotesize \textsuperscript{1} School of Psychology, University of Auckland, New Zealand \newline
\noindent \footnotesize \textsuperscript{2} School of Anthropology and Museum Ethnography, University of Oxford, United Kingdom
\normalsize
\newpage

## Supplementary Methods

### Bayesian models for reanalysis

We provide model formulae for our reanalyses of cross-national correlations, for a general bivariate case with standardised outcome $Y$ and predictor $X$ variables. In the naive regression model without controls for non-independence:

\begin{align}
Y_{i} &\sim \text{Normal}(\mu_{i},\sigma) \\
\mu_{i} &= \alpha + \beta X_{i} \nonumber \\
\alpha &\sim \text{Normal}(0, 0.4) \nonumber \\
\beta &\sim \text{Normal}(0, 0.4) \nonumber \\
\sigma &\sim \text{Exponential}(5) \nonumber
\end{align}

The priors in this model were arrived at by prior predictive checks, with wider priors making predictions beyond the scale of standardised outcome variables and narrower priors being too informative.

To control for spatial non-independence, we add a Gaussian process to this model and feed it a scaled geographic distance matrix $D$ based on Euclidean distances between latitude and longitude coordinates. This distance matrix is computed internally by the _brms_ R package^1^. The Gaussian process uses an exponentiated quadratic covariance kernel, the only covariance kernel currently supported by _brms_. The model formula is:

\begin{align}
Y_{i} &\sim \text{Normal}(\mu_{i},\sigma) \\
\mu_{i} &= \alpha + \kappa_{\text{NATION}[i]} + \beta X_{i} \nonumber \\
\begin{pmatrix}
\kappa_{1} \\ \kappa_{2} \\  ... \\ \kappa_{n}
\end{pmatrix} &\sim \text{MVNormal}
\begin{pmatrix}
\begin{pmatrix}
0 \\ 0 \\ ... \\ 0
\end{pmatrix},\textbf{K}
\end{pmatrix} \nonumber \\
\textbf{K}_{ij} &= sdgp^2 \text{exp} \big (-D_{ij}^2 / (2 lscale^2) \big ) \nonumber \\
\alpha &\sim \text{Normal}(0, 0.4) \nonumber \\
\beta &\sim \text{Normal}(0, 0.4) \nonumber \\
\sigma &\sim \text{Exponential}(5) \nonumber \\
sdgp &\sim \text{Exponential}(5) \nonumber \\
lscale &\sim \text{InverseGamma}(?,?) \nonumber 
\end{align}

where $n$ is the number of nations, and $D^2_{ij}$ reflects the squared Euclidean distances between latitude and longitude coordinates for the $i$-th and $j$-th nations. Notice that the inverse gamma prior on $lscale$ is left undetermined. This is because the _brms_ package intelligently tunes the prior for this parameter based on the covariates of the Gaussian process (see https://betanalpha.github.io/assets/case_studies/gp_part3/part3.html).

To control for cultural phylogenetic non-independence, we manually specify the covariance structure for nation random intercepts using a pre-computed linguistic proximity matrix $L$ (see previous section). The covariance between two nations is assumed to be linearly proportional to the linguistic proximity between those nations. This assumption is justified if we assume that cultural traits evolve neutrally via Brownian motion along a language phylogeny. Although this is a conservative assumption, we follow this approach here rather than another Gaussian Process for several reasons: (1) to avoid model non-convergence when later including a spatial Gaussian Process in the same model, (2) to work around the lack of a coordinate system for linguistic distances, as opposed to latitude and longitude values for geographic distance, and (3) to showcase different ways that researchers can allow nations to covary to control for non-independence.

The non-centered parameterisation of this model is:

\begin{align}
Y_{i} &\sim \text{Normal}(\mu_{i},\sigma) \\
\mu_{i} &= \alpha + z_{\text{NATION}[i]}\sigma_{\alpha}L + \beta X_{i} \nonumber \\
\alpha &\sim \text{Normal}(0, 0.4) \nonumber \\
\beta &\sim \text{Normal}(0, 0.4) \nonumber \\
z_{j} &\sim \text{Normal}(0, 1) \nonumber \\
\sigma_{\alpha} &\sim \text{Exponential}(5) \nonumber \\
\sigma &\sim \text{Exponential}(5) \nonumber 
\end{align}

Finally, we can control for spatial and cultural phylogenetic non-independence simultaneously by including both a Gaussian process over latitude and longitude coordinates _and_ nation random intercepts that covary according to linguistic proximity. The resulting model is as follows:

\begin{align}
Y_{i} &\sim \text{Normal}(\mu_{i},\sigma) \\
\mu_{i} &= \alpha + \kappa_{\text{NATION}[i]} + z_{\text{NATION}[i]}\sigma_{\alpha}L + \beta X_{i} \nonumber \\
\begin{pmatrix}
\kappa_{1} \\ \kappa_{2} \\ ... \\ \kappa_{n}
\end{pmatrix} &\sim \text{MVNormal}
\begin{pmatrix}
\begin{pmatrix}
0 \\ 0 \\ ... \\ 0
\end{pmatrix},\textbf{K}
\end{pmatrix} \nonumber \\
\textbf{K}_{ij} &= sdgp^2 \text{exp} \big (-D_{ij}^2 / (2 lscale^2) \big ) \nonumber \\
\alpha &\sim \text{Normal}(0, 0.4) \nonumber \\
\beta &\sim \text{Normal}(0, 0.4) \nonumber \\
z_{j} &\sim \text{Normal}(0, 1) \nonumber \\
\sigma_{\alpha} &\sim \text{Exponential}(5) \nonumber \\
\sigma &\sim \text{Exponential}(5) \nonumber \\
sdgp &\sim \text{Exponential}(5) \nonumber \\
lscale &\sim \text{InverseGamma}(?,?) \nonumber 
\end{align}

\newpage

## Supplementary Figures

(ref:dagCaption) _A causal directed acyclic graph of spatial and cultural phylogenetic non-independence in cross-national studies._ We are interested in estimating the direct effect of national-level exposure $X$ on national-level outcome $Y$. But these variables are confounded by their common unobserved cause $U$. $U$ is a stand-in for shared environmental, ecological, and geographic causes (e.g. climate, biodiversity, physical topography) and cultural and institutional causes (e.g. cultural norms, technologies, and institutions). In this causal model, we need to condition on $U$ to estimate the direct path from $X$ to $Y$, but we cannot since it is unobserved. However, geographic $G$ and linguistic $L$ relationships between societies influence $U$, since changing a nation's spatial distance to or shared cultural ancestry with other nations will change its environmental and cultural traits. We can thus use $G$ and $L$ to model the covariation between $X$ and $Y$ induced by $U$. Failing to do this and simply estimating the bivariate correlation between $X$ and $Y$ will produce spurious relationships and residuals that are spatially and culturally non-independent around the world.

```{r dag, fig.width=6, fig.height=3.5, fig.cap="(ref:dagCaption)", fig.pos='H', out.width="\\textwidth"}
tar_read(dag)
```

\newpage

(ref:plotReview2Caption) _Analysis-level results from literature review of 100 highly-cited cross-national studies of economic development (red) and cultural values (blue)._ (a) Proportion of analyses accounting for non-independence, overall and split by common methods of controlling for non-independence. Points represent posterior median proportions and ranges represent equal-tailed 95% credible intervals (n = 4308 observations). (b) The association between journal impact factor and the probability that an analysis accounts for non-independence. (c) Estimated trend over time for the probability that an analysis accounts for non-independence. Lines and shaded areas are posterior median regression lines and equal-tailed 50% credible intervals from Bayesian multilevel models (n = 4308 observations). Histograms represent relative counts for individual analyses that did (top) or did not (bottom) account for non-independence. Region FEs = region fixed effects.

```{r plotReview2, fig.width=6.5, fig.height=6, fig.cap="(ref:plotReview2Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotReview2)
```

\newpage

(ref:plotSim1Caption) _Distribution of cross-national correlations from simulation study under strong spatial autocorrelation._ In these simulations, the strength of spatial autocorrelation is set to 0.8 for both outcome and predictor variables, and the true correlation is set to 0. For frequentist regression models, points represent correlation estimates and ranges represent two-tailed 95% confidence intervals (n = 236 observations). For Bayesian regression models, points represent posterior means and ranges represent equal-tailed 95% credible intervals (n = 236 observations). Correlations are ordered by effect size independently in each panel. Red point ranges indicate that the 95% confidence / credible interval excludes zero. Black point ranges indicate that the 95% confidence / credible interval includes zero. SEs = standard errors.

```{r plotSim1, fig.width=9.5, fig.height=5.5, fig.cap="(ref:plotSim1Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotSim1)
```

\newpage

(ref:plotSim2Caption) _Distribution of cross-national correlations from simulation study under strong cultural phylogenetic autocorrelation._ In these simulations, the strength of cultural phylogenetic autocorrelation is set to 0.8 for both outcome and predictor variables, and the true correlation is set to 0. For frequentist regression models, points represent correlation estimates and ranges represent two-tailed 95% confidence intervals (n = 236 observations). For Bayesian regression models, points represent posterior means and ranges represent equal-tailed 95% credible intervals (n = 236 observations). Correlations are ordered by effect size independently in each panel. Red point ranges indicate that the 95% confidence / credible interval excludes zero. Black point ranges indicate that the 95% confidence / credible interval includes zero. SEs = standard errors.

```{r plotSim2, fig.width=9.5, fig.height=5.5, fig.cap="(ref:plotSim2Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotSim2)
```

\newpage

(ref:plotPowerGeo01Caption) _Statistical power estimates for different methods of controlling for spatial non-independence in our simulation study, assuming a small true effect size (r = 0.1)._ For simulated outcome and predictor variables, we systematically varied the strength of spatial autocorrelation, from weak (0.2) to moderate (0.5) to strong (0.8). We simulated 100 datasets per parameter combination assuming a small true correlation between variables (_r_ = 0.1) and fitted different models to each dataset. Statistical power was operationalised as the proportion of models that estimated a slope with a two-tailed 95% confidence / credible interval excluding zero. Points represent raw proportions of models with slope 95% confidence / credible intervals excluding zero, ranges represent two-tailed 95% bootstrap confidence intervals (n = 1000 bootstrap samples), and dashed lines indicate 80% power. Colours indicate whether the strength of autocorrelation for the predictor variable is 0.2 (red), 0.5, (green) or 0.8 (blue). SEs = standard errors.

```{r plotPowerGeo01, fig.width=7, fig.height=5.5, fig.cap="(ref:plotPowerGeo01Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotPowerGeo0.1)
```

\newpage

(ref:plotPowerGeo03Caption) _Statistical power estimates for different methods of controlling for spatial non-independence in our simulation study, assuming a medium true effect size (r = 0.3)._ For simulated outcome and predictor variables, we systematically varied the strength of spatial autocorrelation, from weak (0.2) to moderate (0.5) to strong (0.8). We simulated 100 datasets per parameter combination assuming a medium true correlation between variables (_r_ = 0.3) and fitted different models to each dataset. Statistical power was operationalised as the proportion of models that estimated a slope with a two-tailed 95% confidence / credible interval excluding zero. Points represent raw proportions of models with slope 95% confidence / credible intervals excluding zero, ranges represent two-tailed 95% bootstrap confidence intervals (n = 1000 bootstrap samples), and dashed lines indicate 80% power. Colours indicate whether the strength of autocorrelation for the predictor variable is 0.2 (red), 0.5, (green) or 0.8 (blue). SEs = standard errors.

```{r plotPowerGeo03, fig.width=7, fig.height=5.5, fig.cap="(ref:plotPowerGeo03Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotPowerGeo0.3)
```

\newpage

(ref:plotPowerGeo05Caption) _Statistical power estimates for different methods of controlling for spatial non-independence in our simulation study, assuming a large true effect size (r = 0.5)._ For simulated outcome and predictor variables, we systematically varied the strength of spatial autocorrelation, from weak (0.2) to moderate (0.5) to strong (0.8). We simulated 100 datasets per parameter combination assuming a large true correlation between variables (_r_ = 0.5) and fitted different models to each dataset. Statistical power was operationalised as the proportion of models that estimated a slope with a two-tailed 95% confidence / credible interval excluding zero. Points represent raw proportions of models with slope 95% confidence / credible intervals excluding zero, ranges represent two-tailed 95% bootstrap confidence intervals (n = 1000 bootstrap samples), and dashed lines indicate 80% power. Colours indicate whether the strength of autocorrelation for the predictor variable is 0.2 (red), 0.5, (green) or 0.8 (blue). SEs = standard errors.

```{r plotPowerGeo05, fig.width=7, fig.height=5.5, fig.cap="(ref:plotPowerGeo05Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotPowerGeo0.5)
```

\newpage

(ref:plotPowerLin01Caption) _Statistical power estimates for different methods of controlling for cultural phylogenetic non-independence in our simulation study, assuming a small true effect size (r = 0.1)._ For simulated outcome and predictor variables, we systematically varied the strength of cultural phylogenetic autocorrelation, from weak (0.2) to moderate (0.5) to strong (0.8). We simulated 100 datasets per parameter combination assuming a small true correlation between variables (_r_ = 0.1) and fitted different models to each dataset. Statistical power was operationalised as the proportion of models that estimated a slope with a two-tailed 95% confidence / credible interval excluding zero. Points represent raw proportions of models with slope 95% confidence / credible intervals excluding zero, ranges represent two-tailed 95% bootstrap confidence intervals (n = 1000 bootstrap samples), and dashed lines indicate 80% power. Colours indicate whether the strength of autocorrelation for the predictor variable is 0.2 (red), 0.5, (green) or 0.8 (blue). SEs = standard errors.

```{r plotPowerLin01, fig.width=7, fig.height=5.5, fig.cap="(ref:plotPowerLin01Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotPowerLin0.1)
```

\newpage

(ref:plotPowerLin03Caption) _Statistical power estimates for different methods of controlling for cultural phylogenetic non-independence in our simulation study, assuming a medium true effect size (r = 0.3)._ For simulated outcome and predictor variables, we systematically varied the strength of cultural phylogenetic autocorrelation, from weak (0.2) to moderate (0.5) to strong (0.8). We simulated 100 datasets per parameter combination assuming a medium true correlation between variables (_r_ = 0.3) and fitted different models to each dataset. Statistical power was operationalised as the proportion of models that estimated a slope with a two-tailed 95% confidence / credible interval excluding zero. Points represent raw proportions of models with slope 95% confidence / credible intervals excluding zero, ranges represent two-tailed 95% bootstrap confidence intervals (n = 1000 bootstrap samples), and dashed lines indicate 80% power. Colours indicate whether the strength of autocorrelation for the predictor variable is 0.2 (red), 0.5, (green) or 0.8 (blue). SEs = standard errors.

```{r plotPowerLin03, fig.width=7, fig.height=5.5, fig.cap="(ref:plotPowerLin03Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotPowerLin0.3)
```

\newpage

(ref:plotPowerLin05Caption) _Statistical power estimates for different methods of controlling for cultural phylogenetic non-independence in our simulation study, assuming a large true effect size (r = 0.5)._ For simulated outcome and predictor variables, we systematically varied the strength of cultural phylogenetic autocorrelation, from weak (0.2) to moderate (0.5) to strong (0.8). We simulated 100 datasets per parameter combination assuming a large true correlation between variables (_r_ = 0.5) and fitted different models to each dataset. Statistical power was operationalised as the proportion of models that estimated a slope with a two-tailed 95% confidence / credible interval excluding zero. Points represent raw proportions of models with slope 95% confidence / credible intervals excluding zero, ranges represent two-tailed 95% bootstrap confidence intervals (n = 1000 bootstrap samples), and dashed lines indicate 80% power. Colours indicate whether the strength of autocorrelation for the predictor variable is 0.2 (red), 0.5, (green) or 0.8 (blue). SEs = standard errors.

```{r plotPowerLin05, fig.width=7, fig.height=5.5, fig.cap="(ref:plotPowerLin05Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotPowerLin0.5)
```

\newpage

(ref:plotReplications2Caption) _Reanalysis models fitted to raw data, for economic development (a) and cultural values (b) studies._ Data points are labelled using ISO 3166-1 alpha-2 letter country codes. In all reanalyses, outcome and predictor variables are standardised, making regression slopes comparable to Pearson's correlation coefficients. Lines and shaded areas represent posterior median regression lines and equal-tailed 95% credible intervals. Colours indicate predictions from models with no control (red), spatial control only (green), cultural control only (blue), or both spatial and cultural control (purple). For models with covariates (Adamzyck and Pitt 2009; Gelfand et al. 2011), marginal effects are presented holding all covariates at zero or their reference categories. GDP = gross domestic product; SME = small and medium-sized enterprise.

```{r plotReplications2, fig.width=8, fig.height=8, fig.cap="(ref:plotReplications2Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotReplications2)
```

\newpage

(ref:plotReplications3Caption) _Posterior estimates of Gaussian process functions mapping spatial autocorrelation onto geographic distance from our reanalyses of economic development (a) and cultural values (b) studies._ The y-axis represents the amount of spatial autocorrelation between data points with increasing distance between those points on the x-axis (logged distance in kilometres). Lines and shaded areas represent median posterior spatial autocorrelation functions and equal-tailed 50% and 95% credible intervals.

```{r plotReplications3, fig.width=8.2, fig.height=6, fig.cap="(ref:plotReplications3Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotReplications3)
```

\newpage

(ref:plotReplications4Caption) _Posterior estimates of cultural phylogenetic signal from our reanalyses._ Cultural phylogenetic signal is operationalised as the proportion of national-level variance explained by linguistic proximity between nations. Ridges are full posterior distributions, points are posterior medians, and lines represent equal-tailed 50% and 95% credible intervals. Numbers of observations from the models are as follows, from top to bottom: n = 69, n = 45, n = 103, n = 98, n = 63, n = 89, n = 33, n = 75, n = 67, n = 28, n = 38, and n = 28. 

```{r plotReplications4, fig.width=5, fig.height=5, fig.cap="(ref:plotReplications4Caption)", fig.pos='H', out.width="\\textwidth"}
tar_read(plotReplications4)
```

\newpage

(ref:plotReplications5Caption) _The estimated degree of spatial and cultural phylogenetic non-independence predicts reductions in effect size in our reanalysis._ (a) Higher estimated degrees of spatial autocorrelation at 1,000 km distance predict more pronounced reductions in effect sizes when controlling for non-independence. (b) Higher estimated levels of cultural phylogenetic signal predict more pronounced reductions in effect sizes when controlling for non-independence. In both panels, the y-axis represents the ratio of the effect size when controlling for spatial and cultural non-independence to the original effect size (from naive regression model), and the x-axis represents posterior median model estimates. Regression lines are plotted with equal-tailed 95% credible intervals.

```{r plotReplications5, fig.width=8, fig.height=4, fig.cap="(ref:plotReplications5Caption)"}
tar_read(plotReplications5)
```

\newpage

## Supplementary Tables

(ref:tableSignalCaption) Geographic and cultural phylogenetic signal results for economic development and cultural values variables. Signal estimates reflect the proportion of national-level variance explained by geographic and linguistic covariance matrices. Bayes Factors (BF) reflect support for the hypothesis that the signal estimate differs from zero. HDI = Human Development Index; GDPpc = gross domestic product per capita.

```{r tableSignal, echo=F, warning=F, message=F}
apa_table(
  tar_read(tableSignal),
  booktabs = TRUE,
  placement = "H",
  font_size = "small",
  caption = "(ref:tableSignalCaption)"
)
```

\newpage

(ref:tableSim1Caption) False positive rates in simulation with varying levels of spatial autocorrelation for the predictor variable ($\rho$) and for the outcome variable ($\lambda$), assuming that the true correlation is 0. Numbers represent the total number of analyses, out of 100, that estimated a slope with a two-tailed 95% confidence / credible interval excluding zero, falsely inferring a relationship when none is present.

```{r tableSim1, echo=F, warning=F, message=F}
apa_table(
  tar_read(tableSim1),
  booktabs = TRUE,
  landscape = TRUE,
  escape = FALSE,
  col_spanners = list(`$\\rho$ = 0.2` = c(2,4),
                      `$\\rho$ = 0.5` = c(5,7),
                      `$\\rho$ = 0.8` = c(8,10)),
  col.names = c("Model", rep(paste0("$\\lambda$ = ", 
                                    c(0.2, 0.5, 0.8)), times = 3)),
  placement = "H",
  font_size = "small",
  caption = "(ref:tableSim1Caption)"
)
```

\newpage

(ref:tableSim2Caption) False positive rates in simulation with varying levels of cultural phylogenetic autocorrelation for the predictor variable ($\rho$) and for the outcome variable ($\lambda$), assuming that the true correlation is 0. Numbers represent the total number of analyses, out of 100, that estimated a slope with a two-tailed 95% confidence / credible interval excluding zero, falsely inferring a relationship when none is present.

```{r tableSim2, echo=F, warning=F, message=F}
apa_table(
  tar_read(tableSim2),
  booktabs = TRUE,
  landscape = TRUE,
  escape = FALSE,
  col_spanners = list(`$\\rho$ = 0.2` = c(2,4),
                      `$\\rho$ = 0.5` = c(5,7),
                      `$\\rho$ = 0.8` = c(8,10)),
  col.names = c("Model", rep(paste0("$\\lambda$ = ", 
                                    c(0.2, 0.5, 0.8)), times = 3)),
  placement = "H",
  font_size = "small",
  caption = "(ref:tableSim2Caption)"
)
```

\newpage

(ref:tableReplicationsCaption) Estimated cross-national correlations from our reanalysis of twelve previous cross-national analyses. Numbers are posterior median slopes (equivalent to Pearson's r correlation coefficients) with equal-tailed 95% credible intervals.

```{r tableReplications, echo=F, warning=F, message=F}
apa_table(
  tar_read(tableReplications),
  booktabs = TRUE,
  placement = "H",
  font_size = "footnotesize",
  align = "r",
  landscape = TRUE,
  caption = "(ref:tableReplicationsCaption)"
)
```

\newpage

(ref:tableDeviationsCaption) List of sample size deviations between original analyses and reanalyses.

```{r tableDeviations, echo=F, warning=F, message=F}
apa_table(
  tar_read(tableDeviations),
  booktabs = TRUE,
  placement = "H",
  font_size = "scriptsize",
  longtable = TRUE,
  align = c("m{2cm}", "m{2cm}", "m{2cm}", "m{9cm}"),
  caption = "(ref:tableDeviationsCaption)"
)
```

\newpage

## Supplementary References

1. Bürkner, P.-C. brms: An R package for Bayesian multilevel models using Stan. _Journal of Statistical Software_ **80**, 1–28 (2017).