Week 6 Worksheet Notes

Hierarchies: Multilevel models

Aims

By the end of the session you will learn how to:

  • use survey weights in single-level linear regression models
  • estimate variance in hierarchical/clustered data using robust standard errors in a single-level modelling framework
  • fit random intercept models to hierarchical/clustered data

R packages

# Just in case we get errors asking that the package repositories be explicitly set when installing new packages:

options(repos = list(CRAN = "http://cran.rstudio.com/"))

# Install and load required packages

# We can use the {pacman} package to easily install and load several packages:
# ({pacman} itself may need installing if not yet installed)

pacman::p_load(
  tidyverse, easystats, sjlabelled,
  ggformula, ggeffects, marginaleffects, 
  modelsummary, gtsummary,
  survey, sandwich, lmtest, lme4)            # new modelling packages

Example 1. Österman (2021)

In Worksheet 4 we fitted versions of Models 1-3 reported in Table 3 of Österman (2021) (see also Table A.3. in their Appendix for a more complete reporting on the models).

Comparing the model we fitted to the one reported by the author, our results were very close, but not totally equivalent on several coefficients and standard errors. The main reasons for the divergence had to do with two aspects of the published model that we had disregarded: (1) we did not include survey weights to correct for sampling errors, and (2) we did not allow for intra-group correlation in the standard errors among respondents from the same country (and the same age cohort). We will first implement these two additional steps and compare our final results again to those reported in Österman (2021). Then, we will refit the model in a multilevel/hierarchical framework.

Data preparation

As a first step, let’s import the osterman dataset that underpins the Österman (2021) article (see the Data page on the course website for information about the datasets available in this course). We can download the data to a local folder and load it from there, or we can load it directly from the web (if it works…):

# Import the data
osterman <- data_read("https://cgmoreh.github.io/HSS8005-24/data/osterman.dta")

It’s always a good idea to inspect the dataset after importing, to identify any issues. One option is to check a codebook, for example:

# `View` the codebook:
data_codebook(osterman) |> View()

Remember that some “tagged” NA values that were imported from the Stata data format need changing (see Week 2 worksheet):

osterman <- sjlabelled::zap_na_tags(osterman)

Models 1 and 2

In Worksheet 4 we have already fit Models 1 and 2 and got results very similar - but not identical - to those reported in the published article. Here, let’s first refit these models:

osterman_m1 <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d +                        
                 fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) +       
                 agea + yrbrn + I(agea^2) + I(yrbrn^2) +                                
                 factor(reform_id_num) +                                                
                 yrbrn:factor(reform_id_num) +                                          
                 agea:factor(reform_id_num) +  I(agea^2):factor(reform_id_num),         
              data = osterman)                                                         



osterman_m2 <- lm(trustindex3 ~ reform1_7 + paredu_a_high + reform1_7*paredu_a_high + female + blgetmg_d + 
                 fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) +       
                 agea + yrbrn + I(agea^2) + I(yrbrn^2) +                                
                 factor(reform_id_num) +                                                
                 yrbrn:factor(reform_id_num) +                                          
                 agea:factor(reform_id_num) +  I(agea^2):factor(reform_id_num),         
              data = osterman)

Problem 1: Survey weights

One of the differences between our coefficients and the ones in the published table is due to the fact that we did not account for survey weights, while in Österman (2021), p. 220:

The data are weighted using ESS design weights

To understand in more detail what this means, we need some understanding of how weights are constructed in the European Social Survey (ESS). In a nutshell, ESS data are distributed containing 3 types of weights:

  1. dweight: These are the so-called design weights. Quoting the ESS website: “the main purpose of the design weights is to correct for the fact that in some countries respondents have different probabilities to be part of the sample due to the sampling design used.” These weights have been built to correct for the coverage error, that is, the error created by the different chances that individuals from the target population are covered in the sample frame.

  2. pspwght: These are the post-stratification weights. According the the ESS website, these “are a more sophisticated weighting strategy that uses auxiliary information to reduce the sampling error and potential non-response bias.” These weights have been computed after the data has been collected, to correct from differences between population frequencies observed in the sample and the “true” population frequencies (i.e. those provided by the Labour Force Survey funded by the EU and available on Eurostat). Unlike the design weights, which are based on the probability of inclusion of different groups of individuals in the sample frames, these have been calculated starting from variables that are there in the data, and are really an “adjustment” of the design weight made to reach observed frequencies that match those of the target population.

  3. pweight: These are the population size weights. These weights have the purpose to match the numbers of observations collected in each country to the country populations. They are to be used only when we calculate statistics on multiple countries (for instance, unemployment in Scandinavia). Their value is the same for all observations within the same country.

There is a lot to be said about survey weights and options for dealing with them, which we will not cover in more detail in this course. But as part of this exercise we will get to know some functions that can help with including survey weights and can be extended to include more complex design weights as well. Specifically, we will look at the {survey} package

We start by creating a weighted data object using the svydesign() function from {survey}, which includes the dweight as used in the original article:

## Create weighted data
osterman_w <- svydesign(id = ~1,                 # specifying cluster IDs is needed; ~0 or ~1 means no clusters
                        weights = ~dweight,      # we apply the design weights
                        data = osterman)

We then fit the model using the svyglm() function from the {survey} package and save the model object; note that we specify a design = option with the weighted data object we created earlier:

osterman_m1_w <- svyglm(trustindex3 ~ reform1_7 + female + blgetmg_d +                        
                   fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) +       
                   agea + yrbrn + I(agea^2) + I(yrbrn^2) +                                
                   factor(reform_id_num) +                                                
                   yrbrn:factor(reform_id_num) +                                          
                   agea:factor(reform_id_num) +  I(agea^2):factor(reform_id_num),
                 design = osterman_w, data = osterman) 


osterman_m2_w <- svyglm(trustindex3 ~ reform1_7 + paredu_a_high + reform1_7*paredu_a_high + female + blgetmg_d + 
                   fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) +       
                   agea + yrbrn + I(agea^2) + I(yrbrn^2) +                                
                   factor(reform_id_num) +                                                
                   yrbrn:factor(reform_id_num) +                                          
                   agea:factor(reform_id_num) +  I(agea^2):factor(reform_id_num),    
                 design = osterman_w, data = osterman)

To compare the results from the weighted model to the one we produced earlier, we can check them in a modelsummary() table; we can reuse the list of coefficients to display that we created earlier:

# List and name the models
models <- list(
  "Unweighted M1" = osterman_m1,
  "Unweighted M2" = osterman_m2,
  "Weighted M1" = osterman_m1_w,
  "Weighted M2" = osterman_m2_w)

# Select and name coefficients to be tabulated
coefs <- c("reform1_7" = "Reform",
           "paredu_a_high" = "High parental edu",
           "reform1_7:paredu_a_high" = "Ref x High par edu",
           "femaleFemale" = "Female",
           "blgetmg_d" = "Ethnic minority",
           "fbrneur" = "Foreign-born father, Europe",
           "mbrneur" = "Foreign-born mother, Europe",
           "fnotbrneur" = "Foreign-born father, outside Europe",
           "mnotbrneur" = "Foreign-born mother, outside Europe")

# Tabulate the models
modelsummary::modelsummary(models, stars = TRUE, coef_map = coefs)
Warning in logLik.svyglm(x): svyglm not fitted by maximum likelihood.

Warning in logLik.svyglm(x): svyglm not fitted by maximum likelihood.
Unweighted M1  Unweighted M2 Weighted M1  Weighted M2
Reform 0.063* 0.081** 0.063* 0.083**
(0.027) (0.029) (0.029) (0.031)
High parental edu 0.349*** 0.340***
(0.020) (0.021)
Ref x High par edu −0.049+ −0.043
(0.029) (0.029)
Ethnic minority −0.241*** −0.207*** −0.261*** −0.226**
(0.054) (0.056) (0.066) (0.069)
Foreign-born father, Europe −0.111** −0.106* −0.090* −0.084+
(0.042) (0.044) (0.046) (0.047)
Foreign-born mother, Europe −0.108* −0.115* −0.092* −0.103*
(0.044) (0.046) (0.046) (0.048)
Foreign-born father, outside Europe −0.065 −0.047 −0.053 −0.040
(0.073) (0.076) (0.079) (0.080)
Foreign-born mother, outside Europe −0.110 −0.135+ −0.087 −0.103
(0.078) (0.081) (0.082) (0.083)
Num.Obs. 68796 64960 68796 64960
R2 0.200 0.209 0.206 0.215
R2 Adj. 0.198 0.208 0.204 0.214
AIC 268913.8 252985.1
BIC 270056.1 254138.4 305027.4 286775.7
Log.Lik. −134331.877 −126365.537 −151817.519 −142684.156
RMSE 1.71 1.69 1.71 1.69
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

We can check our results again to those published in Table 3 and Table A.3 in the Appendix to Österman (2021), and we will see that all the coefficients are now the same, with some minor divergences still in the standard errors, which are due to the fact that we did not use cluster robust standard errors to account for intra-group correlations within countries.

Problem 2: Variance estimation using robust clustered errors

Österman (2021:222) write that:

All models are estimated with OLS, … and apply country-by-birth year clustered robust standard errors. [Footnote: An alternative would be to cluster the standard errors on the country level but such an approach would risk to lead to biased standard errors because of too few clusters]

When estimating observations that are structured within groups/clusters (e.g. individual respondents within countries; students within schools/classes; test results or pertaining to the same pupil over time), one of the central assumptions of ordinary least squares regression (OLS) - that our variables are “independent and identically distributed” (iid)) - is breached. Ignoring this hierarchical structure of the data can lead to misleadingly precise coefficient estimates, with inaccurately low standard errors, too narrow confidence intervals, very low p-values and therefore possibly wrong conclusions. Applying a statistical correction for clustered standard errors is a common approach to dealing with this problem. It is a more elementary (but less flexible) way to account for breaches of the iid assumption than fitting a multilevel/hierarchical model.

We will first fit our model using clustered robust standard errors, as done by Österman (2021), then we will check the results against those we would obtain from a multilevel model design to see how the risk related to the low number of country clusters affects the results. For simplicity, we will look at Model 1.

Unfortunately, while estimating clustered robust errors in simple linear (OLS) regression models with software such as Stata is very easy, the lm() command in R was not designed to accommodate this procedure, and as such we need to rely on additional packages and functions. The standard way of applying error corrections is by using the {sandwich} and {lmtest} packages, but they do not handle well the summary objects produced by the {survey} package for weighted estimates. Using the unweighed model we fitted earlier, we could do the following:

# extract variance-covariance matrix with clustered correction
vc_cl <- sandwich::vcovCL(osterman_m1, type='HC1', cluster=~cntry_cohort)

# get coefficients
osterman_m1_cl <- coeftest(osterman_m1, vc_cl)

## Or, the two steps above can be combined into one call:

osterman_m1_cl2 <- coeftest(osterman_m1, 
                          vcov = vcovCL, 
                          type='HC1', 
                          cluster = ~cntry_cohort)

# And we can check that the two have the same result:
all.equal(osterman_m1_cl, osterman_m1_cl2)
[1] TRUE

Let’s tabulate the results for comparison:

# List and name the models
models <- list(
  "Unweighted model" = osterman_m1,
  "Weighted model" = osterman_m1_w,
  "Unweighted model with cluster-robust errors" = osterman_m1_cl)

# Tabulate the models
modelsummary::modelsummary(models, stars = TRUE, coef_map = coefs)
Warning in logLik.svyglm(x): svyglm not fitted by maximum likelihood.
Unweighted model Weighted model Unweighted model with cluster-robust errors
Reform 0.063* 0.063* 0.063*
(0.027) (0.029) (0.029)
Ethnic minority −0.241*** −0.261*** −0.241***
(0.054) (0.066) (0.060)
Foreign-born father, Europe −0.111** −0.090* −0.111*
(0.042) (0.046) (0.046)
Foreign-born mother, Europe −0.108* −0.092* −0.108*
(0.044) (0.046) (0.047)
Foreign-born father, outside Europe −0.065 −0.053 −0.065
(0.073) (0.079) (0.077)
Foreign-born mother, outside Europe −0.110 −0.087 −0.110
(0.078) (0.082) (0.086)
Num.Obs. 68796 68796 68796
R2 0.200 0.206
R2 Adj. 0.198 0.204
AIC 268913.8 406007.8
BIC 270056.1 305027.4 1033594.4
Log.Lik. −134331.877 −151817.519
RMSE 1.71 1.71
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

The error-correction does not have an impact on the coefficients, but it did affect the standard errors. Notice that the standard errors came closer to those in the weighted model without any error correction, showing that the weighting procedure already applies an error correction, albeit not specifically on the cluster variable cntry_cohort.

One deficiency of coeftest()/vcovCL() is that it doesn’t automatically correct for cases when the number of clusters is small. In our case, the number of cluster-cohorts is sufficiently large, so the estimates are good. If that were not the case - if we only had a handful of countries as clusters in the dataset, for example - then one option would be to manually specify the “degrees of freedom” within the coeftest() function, which is the number of the clusters minus 1 (e.g. with five clusters, the DF would be 4, and we would add an additional specification saying df = 4 within coeftest().

Another option is to use the feols() function from the {fixest} package, which fits OLS models that have many fixed effects (i.e. indicator variables representing clusters/categories of a cluster/grouping variable). The function uses the same formula syntax as lm(), but with the extension that allows to specify the variables representing the “fixed effects” (indicator variables for the cluster variable) after a | sign. Apart from the advantage that we can do all the estimations in one step - as opposed to how coeftest() works, requiring an lm() model object as an input first -, another advantage is that the output automatically hides the “fixed effects”, so we do not need to be troubled about potentially very long lists of coefficients that we are not directly interested in when running a model summary(). The function includes a vcov argument for specifying how you want to handle the standard errors, just like in the case of coeftest() earlier. There are also a number of options to choose from, which may be worth reading up about in the function documentation.

An easy way to include information on the clustering directly in the svydesign() function is to add cluster IDs:

# Re-specify the survey design
osterman_w_cl <- svydesign(id = ~cntry_cohort,        # This time we add cntry_cohort as id
                         weights = ~dweight, 
                         data = osterman)

# Re-fit the model with the new survey design
m1_w_cl <- svyglm(trustindex3 ~ reform1_7 + female + blgetmg_d + 
               fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + 
               agea + yrbrn + I(agea^2) + I(yrbrn^2) +
               factor(reform_id_num) +           
               yrbrn:factor(reform_id_num) + agea:factor(reform_id_num) +  
               I(agea^2):factor(reform_id_num),
           design = osterman_w_cl, data = osterman) 

# Tabulate the models

# Add the latest model to the existing list using the append() function to save typing
models <- append(models, 
                 list("Weighted with cluster IDs" = m1_w_cl)
                 )

# Tabulate the models
modelsummary::modelsummary(models, stars = TRUE, coef_map = coefs)
Warning in logLik.svyglm(x): svyglm not fitted by maximum likelihood.

Warning in logLik.svyglm(x): svyglm not fitted by maximum likelihood.
Unweighted model Weighted model Unweighted model with cluster-robust errors  Weighted with cluster IDs
Reform 0.063* 0.063* 0.063* 0.063*
(0.027) (0.029) (0.029) (0.030)
Ethnic minority −0.241*** −0.261*** −0.241*** −0.261***
(0.054) (0.066) (0.060) (0.067)
Foreign-born father, Europe −0.111** −0.090* −0.111* −0.090+
(0.042) (0.046) (0.046) (0.047)
Foreign-born mother, Europe −0.108* −0.092* −0.108* −0.092+
(0.044) (0.046) (0.047) (0.047)
Foreign-born father, outside Europe −0.065 −0.053 −0.065 −0.053
(0.073) (0.079) (0.077) (0.078)
Foreign-born mother, outside Europe −0.110 −0.087 −0.110 −0.087
(0.078) (0.082) (0.086) (0.081)
Num.Obs. 68796 68796 68796 68796
R2 0.200 0.206 0.206
R2 Adj. 0.198 0.204 −209.204
AIC 268913.8 406007.8
BIC 270056.1 305027.4 1033594.4 7093411.5
Log.Lik. −134331.877 −151817.519 −3546009.588
RMSE 1.71 1.71 1.71
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

This final model takes us very close to the model reported in Österman (2021).

Example 3. Österman (2021): Random intercept model

There is a vast literature assessing whether one approach is better suited than the other in different contexts. To get a feel for the issues in question and for a deeper understanding of how these methods are useful, compare the analysis of Cheah (2009), whose results “suggest that modeling the clustering of the data using a multilevel methods is a better approach than fixing the standard errors of the OLS estimate”, to that of McNeish, Stapleton, and Silverman (2017), who discuss a number of cases (focusing on psychology literature) where cluster-robust standard error may be more advantageous than multilevel/hierarchical modelling.

Österman (2021:222) write that “An alternative [to the model above] would be to cluster the standard errors on the country level but such an approach would risk to lead to biased standard errors because of too few clusters”. The number of clusters is an important consideration when choosing to go for a multilevel model. For our aims, we will disregard this warning and go ahead and fit a multilevel model on the osterman data regardless, and we can think about what the differences in the results mean.

We can fit the model in a multilevel framework using the {lme4} package and the lmer() function.

Below we begin by looking at the distribution of the cntry variable that codes the countries where respondents are from. Since the European Social Survey is a cross-national survey, it employs a random sampling approach within each participant country and the resulting cross-national dataset is therefore clustered at the country level.

data_tabulate(osterman, cntry)
Country (cntry) <character>
# total N=68796 valid N=68796

Value |    N | Raw % | Valid % | Cumulative %
------+------+-------+---------+-------------
AT    | 2387 |  3.47 |    3.47 |         3.47
BE    | 5220 |  7.59 |    7.59 |        11.06
DE    | 3035 |  4.41 |    4.41 |        15.47
DK    | 3979 |  5.78 |    5.78 |        21.25
ES    | 6062 |  8.81 |    8.81 |        30.06
FI    | 2505 |  3.64 |    3.64 |        33.71
FR    | 6942 | 10.09 |   10.09 |        43.80
GB    | 5494 |  7.99 |    7.99 |        51.78
GR    | 2065 |  3.00 |    3.00 |        54.78
HU    | 3101 |  4.51 |    4.51 |        59.29
IE    | 6086 |  8.85 |    8.85 |        68.14
IT    | 1502 |  2.18 |    2.18 |        70.32
NL    | 6063 |  8.81 |    8.81 |        79.13
PL    | 5372 |  7.81 |    7.81 |        86.94
PT    | 5987 |  8.70 |    8.70 |        95.65
SE    | 2996 |  4.35 |    4.35 |       100.00
<NA>  |    0 |  0.00 |    <NA> |         <NA>

There are 16 countries in the data.

Initially, it is advisable to first fit some simple, preliminary models, in part to establish a baseline for evaluating larger models. Since in the multilevel framework we are interested in understanding the effect that the clusters have in the data, the simplest model we can fit is a so-called random intercepts model or null model, which models the outcome without any predictors. In the single-level framework, this is equivalent to estimating the overall mean of the outcome variable; in the hierarchical framework, we are modelling the mean of the outcome within each category/cluster of the grouping variable. With lmer() this means a simple addition of a specification of the form (1 | cluster) to the model code we are already familiar with, where cluster is the name of the clustering variable, in our case cntry:

mod_null <- lmer(trustindex3 ~ 1 + (1 | cntry), data = osterman)

The result is the following:

model_parameters(mod_null)
# Fixed Effects

Parameter   | Coefficient |   SE |       95% CI | t(68793) |      p
-------------------------------------------------------------------
(Intercept) |        5.23 | 0.23 | [4.78, 5.69] |    22.47 | < .001

# Random Effects

Parameter             | Coefficient
-----------------------------------
SD (Intercept: cntry) |        0.93
SD (Residual)         |        1.71

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation. Uncertainty intervals for
  random effect variances computed using a Wald z-distribution
  approximation.
random_parameters(mod_null)
# Random Effects

Within-Group Variance        2.94 (1.71)
Between-Group Variance
  Random Intercept (cntry)   0.87 (0.93)
N (groups per factor)
  cntry                        16
Observations                68796

Below we fit a random intercept model with cntry_cohort as the single grouping variable and disregarding survey weights:

# Fit a random intercept model with `cntry_cohort` as the grouping variable
m2 <- lmer(trustindex3 ~ reform1_7 + female + blgetmg_d + 
               fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + 
               agea + yrbrn + I(agea^2) + I(yrbrn^2) +
               factor(reform_id_num) +           
               yrbrn:factor(reform_id_num) + agea:factor(reform_id_num) +  
               I(agea^2):factor(reform_id_num) +
               (1|cntry_cohort),
           data = osterman)
Warning: Some predictor variables are on very different scales: consider
rescaling
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00384146 (tol = 0.002, component 1)
# Fit a random intercept model with only `cntry` as the grouping variable
m2_cntry <- lmer(trustindex3 ~ reform1_7 + female + blgetmg_d + 
               fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + 
               agea + yrbrn + I(agea^2) + I(yrbrn^2) +
               factor(reform_id_num) +           
               yrbrn:factor(reform_id_num) + agea:factor(reform_id_num) +  
               I(agea^2):factor(reform_id_num) +
               (1|cntry),
           data = osterman)
Warning: Some predictor variables are on very different scales: consider
rescaling
# Add the latest models to the existing list using the append() function
models <- append(models, 
                 list("Random intercepts (cntry_cohort)" = m2,
                      "Random intercepts (cntry)" = m2_cntry)
                 )

# Tabulate and compare the models
modelsummary::modelsummary(models, stars = TRUE, coef_map = coefs)
Warning in logLik.svyglm(x): svyglm not fitted by maximum likelihood.
Warning in logLik.svyglm(x): svyglm not fitted by maximum likelihood.
Unweighted model Weighted model Unweighted model with cluster-robust errors  Weighted with cluster IDs Random intercepts (cntry_cohort) Random intercepts (cntry)
Reform 0.063* 0.063* 0.063* 0.063* 0.064* 0.063*
(0.027) (0.029) (0.029) (0.030) (0.028) (0.027)
Ethnic minority −0.241*** −0.261*** −0.241*** −0.261*** −0.241*** −0.241***
(0.054) (0.066) (0.060) (0.067) (0.054) (0.054)
Foreign-born father, Europe −0.111** −0.090* −0.111* −0.090+ −0.111** −0.111**
(0.042) (0.046) (0.046) (0.047) (0.042) (0.042)
Foreign-born mother, Europe −0.108* −0.092* −0.108* −0.092+ −0.108* −0.108*
(0.044) (0.046) (0.047) (0.047) (0.044) (0.044)
Foreign-born father, outside Europe −0.065 −0.053 −0.065 −0.053 −0.064 −0.065
(0.073) (0.079) (0.077) (0.078) (0.073) (0.073)
Foreign-born mother, outside Europe −0.110 −0.087 −0.110 −0.087 −0.110 −0.110
(0.078) (0.082) (0.086) (0.081) (0.078) (0.078)
Num.Obs. 68796 68796 68796 68796 68796 68796
R2 0.200 0.206 0.206
R2 Adj. 0.198 0.204 −209.204
R2 Marg. 0.199 0.126
R2 Cond. 0.200 0.494
AIC 268913.8 406007.8 269877.8 269878.2
BIC 270056.1 305027.4 1033594.4 7093411.5 271029.3 271029.7
ICC 0.0 0.4
Log.Lik. −134331.877 −151817.519 −3546009.588
RMSE 1.71 1.71 1.71 1.70 1.71
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

References

Cheah, B. 2009. “Clustering Standard Errors or Modeling Multilevel Data ?”
McNeish, Daniel, Laura M. Stapleton, and Rebecca D. Silverman. 2017. “On the Unnecessary Ubiquity of Hierarchical Linear Modeling.” Psychological Methods 22(1):114–40. doi: 10.1037/met0000078.
Österman, Marcus. 2021. “Can We Trust Education for Fostering Trust? Quasi-experimental Evidence on the Effect of Education and Tracking on Social Trust.” Social Indicators Research 154(1):211–33. doi: 10.1007/s11205-020-02529-y.