# Import the data
osterman <- data_read("https://cgmoreh.github.io/HSS8005/data/osterman_t3.dta")Week 5 Application Lab
Hierarchies: Hierarchical data structures and multilevel modelling
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
Setup
In Week 1 you set up R and RStudio, and an RProject folder (we called it “HSS8005_labs”) with an .R script and a .qmd or .Rmd document in it (we called these “Lab_1”). Ideally, you saved this on a cloud drive so you can access it from any computer (e.g. OneDrive). You will be working in this folder.
Create a new Quarto markdown file (.qmd) for this session (e.g. “Lab_5.qmd”) and work in it to complete the exercises and report on your final analysis.
R packages
Install (if needed) and load the following packages
Packages used previously:
New packages:
Data source
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…):
It’s always a good idea to inspect the dataset after importing, to identify any issues. One option is to check a codebook using the datawizard::data_codebook() function; to allow for more levels of categorical variables to be tabulated, we can change the number of maximum values shown from the default 10 to something larger:
# `View` the codebook:
data_codebook(osterman, max_values = 20) |> View()We do notice a few strange values for some variables. The so-called “tagged” NA/missing values are causing the ordering of the values to be unexpected: because of letters appearing as part of a numeric value scale, the scale is automatically following “alphabetic sorting” rather than “natural sorting”. R is not designed to handle “tagged” NA/missing values, but these are commonly used in other statistical software, and the original dataset was imported from a Stata/.dta format. In most cases this will not pose an issue, but it’s safer to convert all “tagged” missing values to standard NAs. We can easily do that with the zap_na_tags() function from the sjlabelled package:
osterman <- sjlabelled::zap_na_tags(osterman)We can test the results with another call to datawizard::data_codebook(). The ordering of the variables should be okay now.
Problem 1: Polynomials
In the footnotes to Table 3 and Table A.3 in the Appendix, the author tells us:
All models include reform FEs [Fixed Effects], a general quadratic birth year trend, and reform-specific trends for birth year (linear) and age (quadratic). Additional controls: foreign-born parents and ESS round dummies.
The tables in the Appendix (Table A.3) “report all of the coefficients for the control variables, except for the age and birth year controls. Since the latter variables are interacted with all of the reform fixed effects, they are not interpretable as single coefficients”.
In the main text, the author further explains some of the reasoning behind their modelling choices:
One dilemma for the design is that there has been a trend of increasing educational attainment throughout the studied time period, which means that the reform-windows of treated and non-treated cohorts will also pick up the effects of this trend. To counter this, [the list of covariates] includes a general quadratic birth year trend, reform-specific linear controls for birth year and reform-specific quadratic age trends. The quadratic terms are included to allow sufficient flexibility in absorbing possible non-linear trends of increasing education within the reform-window of seven treated and seven untreated birth year cohorts. … The reform-fixed effects are also essential because they imply that only the within-reform-window variation is used to estimate the effects and between-reform differences are factored out, such as pre-reform differences in social trust. (Österman 2021:221–22)
Before we fit the model, some of the concepts in the quotation need unpacking. A quadratic term is a second-degree polynomial term: put simply, it’s the square of the variable concerned. The quadratic of a variable such as \(age\) is therefore \(age^2\), or \(age \times age\). In other words, it is like a variable’s “interaction” with itself. Because of this, there are several ways in which the quadratic terms to be included in a model can be specified:
- We could create the quadratic terms as new variables, and include those in the model. Effectively, we would be creating new variables/columns in the dataset and using those in the model. We could do that as shown below:
# Create new quadratic variables by multiplication with themselves and add them to the dataset, saving it as a new data object:
osterman <- osterman |>
mutate(agea_quad = agea*agea, # quadratic age variable
yrbrn_quad = yrbrn*yrbrn) |> # quadratic birth-year variable
sjlabelled::var_labels(agea_quad = "Age (quadratic)", # we can label the variable if we want
yrbrn_quad = "Birth year (quadratic)")
# We now have two additional variables in the dataset; we can fit a model using those:
m1_prequad <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d + # main covariates reported in Table 3
fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + # additional controls for foreign-born parents and ESS Round dummies (we treat `round` as a factor for this)
agea + yrbrn + agea_quad + yrbrn_quad + # general quadratic birth year trend and quadratic age
factor(reform_id_num) + # reform fixed effects dummies
yrbrn:factor(reform_id_num) + # reform-specific birth year trend
agea:factor(reform_id_num) + agea_quad:factor(reform_id_num), # reform-specific quadratic age trend
data = osterman) # the new expanded dataset- We can get the same results by creating the quadratic terms directly as part of the modelling function. The one thing we should keep in mind is that if we want to include standard mathematical operations within a formula function, we need to isolate or insulate the operation from
R’s formula parsing code using theI()function, which returns the contents ofI(...)“as.is”. The model formula would then be:
# Create quadratic terms internally as part of the modelling function:
m1_funquad <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d + # main covariates reported in Table 3
fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) + # additional controls for foreign-born parents and ESS Round dummies (we treat `round` as a factor for this)
agea + yrbrn + I(agea^2) + I(yrbrn^2) + # general quadratic birth year trend and quadratic age
factor(reform_id_num) + # reform fixed effects dummies
yrbrn:factor(reform_id_num) + # reform-specific birth year trend
agea:factor(reform_id_num) + I(agea^2):factor(reform_id_num), # reform-specific quadratic age trend
data = osterman) # the original dataset- In the two previous options, the quadratic terms will be correlated with the original variables. To avoid this by relying on so-called orthogonal polynomials we should use the
poly()function. We can also fit the same correlated polynomial model as the ones above by adding theraw = TRUEoption to thepoly()function. In the code below, we will fit the correlated version first, then the orthogonal version (This stackoverflow discussion explains in more detail the difference between the two options):
m1_polyraw <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d +
fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) +
poly(agea, 2, raw = TRUE) + poly(yrbrn, 2, raw = TRUE) +
factor(reform_id_num) +
yrbrn:factor(reform_id_num) +
agea:factor(reform_id_num) + poly(agea, 2, raw = TRUE):factor(reform_id_num),
data = osterman)
m1_polyorth <- lm(trustindex3 ~ reform1_7 + female + blgetmg_d +
fbrneur + mbrneur + fnotbrneur + mnotbrneur + factor(essround) +
poly(agea, 2) + poly(yrbrn, 2) +
factor(reform_id_num) +
yrbrn:factor(reform_id_num) +
agea:factor(reform_id_num) + poly(agea, 2):factor(reform_id_num),
data = osterman)For a side-by-side overview comparison of the models we fitted so far, we can use one of the many model summary tabulation functions that exist in various packages. Some popular options are stargazer::stargazer(), jtools::export_summs(), sjPlot::tab_model(), modelsummary::modelsummary() or parameters::compare_models from the easystats ecosystem. Below we will use modelsummary(), which produces publication-ready summary tables in HTML format, but which can easily be exported to other formats, such as Microsoft Word or PDF:
# It's cleaner to first make a list of the models we want to summarise; we can even name them:
models <- list(
"Pre-calculated quadratic" = m1_prequad,
"Within-function quadratic" = m1_funquad,
"poly() with raw coding" = m1_polyraw,
"poly() with default orthogonal coding" = m1_polyorth
)
# modelsummary table with stars for p-values added
modelsummary::modelsummary(models, stars = TRUE)The results from the modelsummary() are not shown here because it’s a long and ugly table, but it’s useful for perusing to compare the results across the different models. We do notice some differences in the affected variables between the orthogonal-codes version and the other versions. It’s worth noting, however, that the Stata routine used by the author fitted correlated/raw coded polynomials, so the orthogonal version from the output above is just for a comparison and we will not use it going forward. We generally want our transformed (polynomial) variables to be correlated with the original variables, as they are in fact measuring the same thing.
For a cleaner table showing only the results included in Table A.3 in the Appendix to Österman (2021), we can use the coef_map or the coef_omit option in modelsummary() and only include m1_funquad, which will be the polynomial fitting routine that we will use going forward.
# It's cleaner to first make a vector of the coefficients we wish to include; we can name the coefficients as they appear in Table A.3; note that we also leave out the Intercept, as in the published table:
coefs <- c("reform1_7" = "Reform",
"female" = "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",
"factor(essround)2" = "ESS Round 2",
"factor(essround)3" = "ESS Round 3",
"factor(essround)4" = "ESS Round 4",
"factor(essround)5" = "ESS Round 5",
"factor(essround)6" = "ESS Round 6",
"factor(essround)7" = "ESS Round 7",
"factor(essround)8" = "ESS Round 8",
"factor(essround)9" = "ESS Round 9")
# Then we pass the vector to coef_map to select the coefficients to print
modelsummary::modelsummary(list(m1_funquad), stars = TRUE, coef_map = coefs)| (1) | |
|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |
| Reform | 0.063* |
| (0.027) | |
| ESS Round 2 | 0.059 |
| (0.045) | |
| ESS Round 3 | 0.162* |
| (0.075) | |
| ESS Round 4 | 0.243* |
| (0.108) | |
| ESS Round 5 | 0.360* |
| (0.144) | |
| ESS Round 6 | 0.397* |
| (0.179) | |
| ESS Round 7 | 0.449* |
| (0.212) | |
| ESS Round 8 | 0.655** |
| (0.246) | |
| ESS Round 9 | 0.816** |
| (0.283) | |
| Num.Obs. | 68796 |
| R2 | 0.200 |
| R2 Adj. | 0.198 |
| AIC | 268913.8 |
| BIC | 270056.1 |
| Log.Lik. | -134331.877 |
| RMSE | 1.71 |
We can compare our table to the published one; there will, of course, be some differences, as our model specification is still not the same as that of the original paper. One source of the inconsistencies are the use of “survey weights” in the original article.
Problem 2: Survey weights
First, let’s refit Models 1 and 2:
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) 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:
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.
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.
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)| Unweighted M1 | Unweighted M2 | Weighted M1 | Weighted M2 | |
|---|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||||
| Reform | 0.063* | 0.081** | 0.063* | 0.083** |
| (0.027) | (0.029) | (0.029) | (0.031) | |
| Female | 0.058*** | 0.059*** | 0.061*** | 0.063*** |
| (0.013) | (0.013) | (0.014) | (0.014) | |
| 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 |
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 3: Variance estimation using clustered robust standard 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.
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.
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.
In the original article the clustering to which correction was applied is at the birth-year level by country. We can create such a variable and call it cntry_cohort as in the original. One approach, using the datawizard::data_unite() function, is to create a new string variable that pastes together the values of the two variables of interest:
osterman <- osterman |>
data_unite(new_column = "cntry_cohort", # a new variable called `cntry_cohort`
select = c("cntry", "yrbrn"), # will be created from the values of `cntry` and `yrbrn`
append = TRUE) # and appended to the existing datasetWe now have the needed clustering variable. 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)
# We can check that the two have the same result:
all.equal(osterman_m1_cl, osterman_m1_cl2)[1] TRUE# The output should be TRUELet’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)| Unweighted model | Weighted model | Unweighted model with cluster-robust errors | |
|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||
| Reform | 0.063* | 0.063* | 0.063* |
| (0.027) | (0.029) | (0.029) | |
| Female | 0.058*** | 0.061*** | 0.058*** |
| (0.013) | (0.014) | (0.015) | |
| 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 | |
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)| Unweighted model | Weighted model | Unweighted model with cluster-robust errors | Weighted with cluster IDs | |
|---|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||||
| Reform | 0.063* | 0.063* | 0.063* | 0.063* |
| (0.027) | (0.029) | (0.029) | (0.030) | |
| Female | 0.058*** | 0.061*** | 0.058*** | 0.061*** |
| (0.013) | (0.014) | (0.015) | (0.015) | |
| 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 | |
This final model takes us very close to the model reported in Österman (2021).
Problem 4: Random intercept model
Ö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.71random_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 68796Below 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)
# 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)
# 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)| Unweighted model | Weighted model | Unweighted model with cluster-robust errors | Weighted with cluster IDs | Random intercepts (cntry_cohort) | Random intercepts (cntry) | |
|---|---|---|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||||||
| 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) | |
| Female | 0.058*** | 0.061*** | 0.058*** | 0.061*** | 0.058*** | 0.058*** |
| (0.013) | (0.014) | (0.015) | (0.015) | (0.013) | (0.013) | |
| 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 | |
Problem 5: Random coefficients/slopes model
As a take-home task, think about a potential (i.e. theoretically substantiated) random coefficients/slopes model that could fit these data. Try fitting the model and compare the results to those from the random intercepts model above.
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.