Baby steps in Bayes: Recoding predictors and homing in on specific comparisons
Interpreting models that take into account a host of possible
interactions between predictor variables can be a pain, especially
when some of the predictors contain more than two levels.
In this post, I show how I went about fitting and then making sense
of a multilevel model containing a three-way interaction between
its categorical fixed-effect predictors. To this end, I used
which makes it relatively easy to fit Bayesian models using a notation
that hardly differs from the one used in the popular
I won’t discuss the Bayesian bit much here (I don’t think it’s too important),
and I will instead cover the following points:
- How to fit a multilevel model with
R’s default way of handling categorical predictors (treatment coding).
- How to interpret this model’s fixed parameter estimates.
- How to visualise the modelled effects.
- How to recode predictors to obtain more useful parameter estimates.
- How to extract information from the model to home in on specific comparisons.
For a longitudinal project, 328 children wrote narrative and argumentative texts in Portuguese at three points in time. About a third of the children hailed from Portugal, about a third were children of Portuguese heritage living in the French-speaking part of Switzerland, and about a third were children of Portuguese heritage living in the German-speaking part of Switzerland. Not all children wrote both kinds of texts at all three points in time, and 1,040 texts were retained for the analysis. For each text, we computed the Guiraud index, which is a function of the number of words (tokens) and the number of different words (types) in the texts. Higher values are assumed to reflect greater diversity in vocabulary use.
If you want to know more about this project, check out Bonvin et al. (2018), Lambelet et al. (2017a,b) and Vanhove (forthcoming); you’ll find the references at the bottom of this page.
Read in the data:
Let’s also plot the data. Incidentally, and contrary to popular belief,
I don’t write
ggplot code such as this from scratch. What you see is
the result of drawing and redrawing (see comments).
Figure 1. The texts’ Guiraud values by time of data collection, text type, and language background.
A multilevel model with treatment coding
Our data are nested: Each child wrote up to 6 texts, and
the data were collected in classes, with each child belong to one class.
It’s advisable to take such nesting into account since you may end
up overestimating your degree of certainty about the results otherwise.
I mostly use
glmer() functions to handle such data,
but as will become clearer in a minute,
brm() function offers
some distinct advantages. So let’s load that package:
Fitting the model
We’ll fit a model with a three-way fixed-effect interaction between
Group as well as with by-
Child and by-
random intercepts. In order to take into account the possibility that
children vary in the development of their lexical diversity, we add
a random slope of
Child, and in order to take into account
the possibility that their lexical diversity varies by text type, we
do the same for
TextType. Similarly, we add by-
Class random slopes
Interpreting the parameter estimates
Family: gaussian Links: mu = identity; sigma = identity Formula: Guiraud ~ Time * TextType * Group + (1 + TextType + Time | Class) + (1 + TextType + Time | Child) Data: d (Number of observations: 1040) Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; total post-warmup samples = 4000 [[...]] Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 5.25 0.13 5.01 5.51 2151 1.00 TimeT2 0.57 0.13 0.30 0.82 2017 1.00 TimeT3 0.91 0.14 0.63 1.17 2509 1.00 TextTypenarrative -0.27 0.17 -0.61 0.07 1884 1.00 GroupPortugueseMFrench -1.07 0.17 -1.39 -0.73 2131 1.00 GroupPortugueseMGerman -1.30 0.16 -1.63 -0.98 2422 1.00 TimeT2:TextTypenarrative -0.23 0.16 -0.55 0.08 2295 1.00 TimeT3:TextTypenarrative -0.06 0.17 -0.38 0.27 2489 1.00 TimeT2:GroupPortugueseMFrench -0.55 0.19 -0.91 -0.18 2508 1.00 TimeT3:GroupPortugueseMFrench -0.38 0.19 -0.73 0.01 2790 1.00 TimeT2:GroupPortugueseMGerman -0.63 0.18 -0.99 -0.29 2469 1.00 TimeT3:GroupPortugueseMGerman -0.22 0.19 -0.61 0.16 2612 1.00 TextTypenarrative:GroupPortugueseMFrench 0.14 0.24 -0.33 0.63 2027 1.00 TextTypenarrative:GroupPortugueseMGerman 0.19 0.23 -0.28 0.66 2348 1.00 TimeT2:TextTypenarrative:GroupPortugueseMFrench 0.36 0.24 -0.10 0.82 2703 1.00 TimeT3:TextTypenarrative:GroupPortugueseMFrench 0.25 0.24 -0.22 0.71 2652 1.00 TimeT2:TextTypenarrative:GroupPortugueseMGerman 0.54 0.25 0.04 1.04 4000 1.00 TimeT3:TextTypenarrative:GroupPortugueseMGerman 0.26 0.25 -0.23 0.74 3065 1.00 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sigma 0.60 0.02 0.55 0.65 459 1.00 [[...]]
The output looks pretty similar to what we’d obtain when using
but let’s review what these estimates actually refer to.
R uses treatment coding. This entails that the
refers to a specific combination of factors: the combination of all
reference levels. Again by default, the reference levels are chosen
Timeconsists of three levels (
T3); for alphabetical reasons,
T1is chosen as the default reference level.
Groupalso consists of three levels (
monolingual Portugueseis chosen as the default level.
TextTypeconsists of two levels (
argumentativeis the default reference level.
Intercept, then, shows the modelled mean Guiraud value of
monolingual Portuguese children at
If you’re unsure which factor level was used as the reference level, you can use the
The reference level is the one in whose rows only zeroes occur.
Crucially, all other estimated effects are computed with respect to this
intercept. That is,
TimeT2 (0.57) shows the difference between T1 and T2
for monolingual Portuguese children writing argumentative texts.
TimeT3 (0.91) shows the difference between T1 and T3
for monolingual Portuguese children writing argumentative texts,
TextTypenarrative (-0.27) shows the difference between
the mean Guiraud values of argumentative and narrative texts written
by monolingual Portuguese children writing at T1.
The texts written by the Portuguese-German and Portuguese-French bilinguals
don’t enter into these estimates.
Now, it’s possible to piece together the mean values associated with each combination of predictor values, but questions such as the following remain difficult to answer with just these estimates at hand:
- What’s the overall difference between T2 and T3 and its uncertainty?
- What’s the overall difference between the Guiraud values of texts written by Portuguese-French and Portuguese-German children and its uncertainty?
We’ll tackle these questions in a minute; for now, the point is merely that the estimated parameters above all refer to highly specific comparisons that may not be the most relevant.
Plotting the fitted values and the uncertainty about them
When working with
brms, it’s relatively easy to obtain
the modelled average outcome value for each combination
of the predictor variables as well as a measure of the uncertainty associated
First construct a small data frame containing the unique combinations of predictor variables in our dataset:
If you feed the model (here:
the data frame we’ve just created (
d_pred) to the
it outputs the modelled mean estimate for each combination of
predictor values (
Estimate), the estimated error of this mean estimate (
and a 95% uncertainty interval about the estimate (
One more thing: The
re_formula = NA line specifies that we do not want the
variability associated with the by-
Class and by-
Child random effects to affect
the estimates and their uncertainty. This is what I typically want.
Group Time TextType Estimate Est.Error Q2.5 Q97.5 1 monolingual Portuguese T1 argumentative 5.254257 0.1275708 5.007925 5.511725 2 monolingual Portuguese T1 narrative 4.988935 0.1736867 4.632209 5.333172 3 monolingual Portuguese T2 argumentative 5.819746 0.1370042 5.546781 6.096290 4 monolingual Portuguese T2 narrative 5.320470 0.1830628 4.948748 5.670897 5 monolingual Portuguese T3 argumentative 6.166700 0.1551590 5.855077 6.460323 6 monolingual Portuguese T3 narrative 5.843974 0.1869116 5.465407 6.218007 7 Portuguese-French T1 argumentative 4.189003 0.1138625 3.969089 4.417871 8 Portuguese-French T1 narrative 4.066947 0.1600089 3.759097 4.394976 9 Portuguese-French T2 argumentative 4.206893 0.1211772 3.973083 4.454998 10 Portuguese-French T2 narrative 4.215198 0.1547604 3.912582 4.526929 11 Portuguese-French T3 argumentative 4.725361 0.1280175 4.483346 4.985827 12 Portuguese-French T3 narrative 4.796162 0.1575997 4.487911 5.103481 13 Portuguese-German T1 argumentative 3.949633 0.1023994 3.748789 4.147540 14 Portuguese-German T1 narrative 3.870226 0.1504619 3.580507 4.159204 15 Portuguese-German T2 argumentative 3.883007 0.1130342 3.658081 4.099607 16 Portuguese-German T2 narrative 4.108290 0.1569674 3.804663 4.411843 17 Portuguese-German T3 argumentative 4.639510 0.1297738 4.387543 4.887834 18 Portuguese-German T3 narrative 4.759506 0.1483870 4.470633 5.046205
So where do these estimates and uncertainty intervals come from?
In the Bayesian approach, every model parameter hasn’t got just
one estimate but an entire distribution of estimates. Moreover,
everything that depends on model parameters also has an entire
distribution of estimates associated with it. The mean modelled
outcome values per cell depend on the model parameters, so they, too,
have entire distributions associated with them.
fitted() function summarises these distributions for us:
it returns their means as
Estimate, their standard deviations
Est.Error and their 2.5th and 97.5 percentiles as
If so inclined, you can generate these distributions yourself
This returns matrix of 4000 rows and 18 columns.
4000 is the number of ‘post-warmup samples’ (see the output of
18 is the number of combinations of predictor values in
The first column of
posterior_fit contains the distribution associated with the first
d_pred. If you compute its mean, standard deviation and 2.5th and 97.5th
percentiles, you end up with the same numbers as above:
Or similarly for the 10th row of
d_pred (Portuguese-French, T2, narrative):
At the moment, using
posterior_linpred() has no added value, but it’s good to know
where these numbers come from.
Let’s draw a graph showing these modelled averages and the uncertainty about them. 95% uncertainty intervals are typically used, but they may instill dichotomous thinking. To highlight that such an interval highlights but two points on a continuum, I’m tempted to add 80% intervals as well:
Group Time TextType Estimate Est.Error Q2.5 Q10 Q90 Q97.5 1 monolingual Portuguese T1 argumentative 5.254257 0.1275708 5.007925 5.093502 5.414252 5.511725 2 monolingual Portuguese T1 narrative 4.988935 0.1736867 4.632209 4.769862 5.207548 5.333172 3 monolingual Portuguese T2 argumentative 5.819746 0.1370042 5.546781 5.649066 5.992942 6.096290 4 monolingual Portuguese T2 narrative 5.320470 0.1830628 4.948748 5.091014 5.553054 5.670897 5 monolingual Portuguese T3 argumentative 6.166700 0.1551590 5.855077 5.966155 6.361068 6.460323 6 monolingual Portuguese T3 narrative 5.843974 0.1869116 5.465407 5.607795 6.079387 6.218007 7 Portuguese-French T1 argumentative 4.189003 0.1138625 3.969089 4.042663 4.332786 4.417871 8 Portuguese-French T1 narrative 4.066947 0.1600089 3.759097 3.862814 4.270945 4.394976 9 Portuguese-French T2 argumentative 4.206893 0.1211772 3.973083 4.058347 4.360160 4.454998 10 Portuguese-French T2 narrative 4.215198 0.1547604 3.912582 4.020499 4.407168 4.526929 11 Portuguese-French T3 argumentative 4.725361 0.1280175 4.483346 4.567181 4.886452 4.985827 12 Portuguese-French T3 narrative 4.796162 0.1575997 4.487911 4.596159 4.994741 5.103481 13 Portuguese-German T1 argumentative 3.949633 0.1023994 3.748789 3.815955 4.079486 4.147540 14 Portuguese-German T1 narrative 3.870226 0.1504619 3.580507 3.673275 4.059239 4.159204 15 Portuguese-German T2 argumentative 3.883007 0.1130342 3.658081 3.738352 4.023601 4.099607 16 Portuguese-German T2 narrative 4.108290 0.1569674 3.804663 3.908339 4.309186 4.411843 17 Portuguese-German T3 argumentative 4.639510 0.1297738 4.387543 4.471301 4.802937 4.887834 18 Portuguese-German T3 narrative 4.759506 0.1483870 4.470633 4.572786 4.951090 5.046205
And now for the graph:
Figure 2. The modelled mean Guiraud values and their uncertainty (thick vertical lines: 80% interval; thin vertical lines: 95% interval).
A model with more sensible coding
Tailoring the coding of categorical predictors to the research questions
summary() output for
m_default was difficult to interpret because
treatment coding was used. However, we can override this default behaviour to
end up with estimates that are more readily and more usefully interpretable.
The first thing we can do is to override the default refence level.
Figure 1 showed that the Guiraud values at T2 tend to be somewhere midway
between those at T1 and T3, so we can make the intercept estimate more
representative of the dataset as a whole by making T2 the reference level
Time rather than T1. A benefit of doing so is that we will now have
TimeT3 that specify the difference
between T1-T2 and T2-T3, respectively. In other words, the estimated parameters
will directly reflect the progression from data collection to data collection.
(Before, the parameter estimates specified the differences between T1-T2 and T1-T3,
so a direct estimate for T2-T3 was lacking.)
Second, there’s no reason for preferring
as the reference level. If we sum-code this predictor, the intercept reflects
the grand mean of the argumentative and narrative texts (at T2), and the estimated
parameter then specifies how far the mean Guiraud value of each text type
is removed from this mean:
Similarly, there are a couple of reasonable ways to choose the reference
Group when using treatment coding. But you can also sum-code
this predictor so that the intercept reflects the grand mean of the Guiraud
values of texts written by monolingual Portuguese and bilingual Portuguese-French
and Portuguese-German kids (at T2).
Refitting the model
No difference here.
Interpreting the parameter estimates
Family: gaussian Links: mu = identity; sigma = identity Formula: Guiraud ~ Time * TextType * Group + (1 + TextType + Time | Class) + (1 + TextType + Time | Child) Data: d (Number of observations: 1040) Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; total post-warmup samples = 4000 [[...]] Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 4.59 0.06 4.47 4.72 1809 1.00 TimeT1 -0.21 0.06 -0.32 -0.09 2735 1.00 TimeT3 0.56 0.06 0.44 0.68 2411 1.00 TextType1 0.05 0.05 -0.05 0.14 1617 1.00 Group1 0.98 0.10 0.78 1.17 1426 1.00 Group2 -0.38 0.09 -0.55 -0.19 1332 1.00 TimeT1:TextType1 0.03 0.05 -0.07 0.14 2685 1.00 TimeT3:TextType1 -0.02 0.05 -0.12 0.07 2595 1.00 TimeT1:Group1 -0.24 0.09 -0.41 -0.07 1801 1.00 TimeT3:Group1 -0.13 0.09 -0.31 0.06 1986 1.00 TimeT1:Group2 0.12 0.09 -0.05 0.29 1899 1.00 TimeT3:Group2 -0.01 0.09 -0.18 0.17 2104 1.00 TextType1:Group1 0.21 0.07 0.07 0.34 1446 1.00 TextType1:Group2 -0.04 0.06 -0.17 0.08 1495 1.00 TimeT1:TextType1:Group1 -0.15 0.07 -0.29 -0.01 2201 1.00 TimeT3:TextType1:Group1 -0.07 0.07 -0.20 0.06 2206 1.00 TimeT1:TextType1:Group2 0.03 0.07 -0.11 0.17 2437 1.00 TimeT3:TextType1:Group2 -0.01 0.07 -0.14 0.12 2353 1.00 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sigma 0.60 0.02 0.55 0.65 433 1.01 [[...]]
Intercept reflects the grand mean of the Guiraud values for
both argumentative and narrative texts for all three groups written at T2.
TimeT1 estimate (-0.21) shows the difference between T1 and T2
averaged over all text types and all groups (0.21 points worse at T1);
TimeT3 estimate (0.56) shows the difference between T2 and T3
averaged over all text types and all groups (0.56 points better at T3).
TextType1 (0.05) shows that the mean Guiraud value of one text type (still written at T2!)
averaged over all groups is 0.05 points higher than the grand mean; and by
implication that the mean Guiraud value of the other text type is 0.05 lower than the grand mean.
To find out which text type is which, use
argumentative is coded as 1, it’s the argumentative texts that have
the higher Guiraud values at T2.
Group1 (0.98) shows that one group has higher-than-average Guiraud values
averaged across text types at T2, whereas
Group2 (-0.38) shows that another group
has a mean Guiraud value that lies 0.38 points below the average at T2. By implication,
the third group’s mean Guiraud value lies 0.60 points below average ((0.98-0.38-0.60)/3 = 0).
To see which group is which, use
monolingual Portuguese is ‘1’ for the purposes of
the purposes of
Portuguese-German is the third group.
We can double-check these numbers by generating the modelled mean values for each predictor value combination:
Group Time TextType Estimate Est.Error Q2.5 Q97.5 1 monolingual Portuguese T1 argumentative 5.256561 0.1478317 4.972034 5.555039 2 monolingual Portuguese T1 narrative 4.989113 0.1653303 4.658477 5.300454 3 monolingual Portuguese T2 argumentative 5.825202 0.1429521 5.548448 6.111764 4 monolingual Portuguese T2 narrative 5.317109 0.1593537 4.999305 5.632874 5 monolingual Portuguese T3 argumentative 6.168232 0.1693523 5.840185 6.509393 6 monolingual Portuguese T3 narrative 5.843890 0.1743146 5.499537 6.185465 7 Portuguese-French T1 argumentative 4.186961 0.1300522 3.932579 4.443590 8 Portuguese-French T1 narrative 4.057777 0.1628389 3.740767 4.376938 9 Portuguese-French T2 argumentative 4.210701 0.1241935 3.974560 4.464941 10 Portuguese-French T2 narrative 4.205977 0.1419996 3.941287 4.500929 11 Portuguese-French T3 argumentative 4.730486 0.1375478 4.471109 5.013081 12 Portuguese-French T3 narrative 4.790799 0.1532786 4.492967 5.097257 13 Portuguese-German T1 argumentative 3.949099 0.1160846 3.718169 4.177300 14 Portuguese-German T1 narrative 3.874875 0.1450871 3.590572 4.162969 15 Portuguese-German T2 argumentative 3.876068 0.1162319 3.643204 4.101122 16 Portuguese-German T2 narrative 4.111929 0.1464768 3.823828 4.393865 17 Portuguese-German T3 argumentative 4.630554 0.1345938 4.361076 4.890592 18 Portuguese-German T3 narrative 4.750127 0.1408441 4.470444 5.026297
Some sanity checks:
(1) Intercept = 4.59 = grand mean at T2:
(2) TimeT3 = 0.56 = T2/T3 difference across texts and groups:
(3) Portuguese-German lies 0.60 below average at T2 across texts:
I won’t plot the modelled averages and their uncertainty, because the result will be the same as before: Recoding the predictors in this way doesn’t affect the modelled averages per cell; it just makes the summary output easier to parse.
Homing in on specific comparisons
Finally, let’s see how we can target some specific comparisons without having to refit the model several times. A specific comparison you might be interested in could be “How large is the difference in Guiraud scores for narrative texts written by Portuguese-French bilinguals between T1 and T2?” Or a more complicated one: “How large is the difference in the progression from T1 to T3 for argumentative texts between Portuguese-French and Portuguese-German children?”
To answer such questions, we need to generate the distribution of the modelled averages per predictor value combination:
Question 1: Progression T1-T2 for narrative texts, Portuguese-French bilinguals?
This question requires us to compare the modelled average
for narrative texts written by Portuguese-French bilinguals at T2
to that of the narrative texts written by Portuguese-French bilinguals at T1.
The first combination of predictor values can be found in row 10 in
so the corresponding estimates are in column 10 in
The second combination of predictor values can be found in row 8 in
so the corresponding estimates are in column 8 in
Now compute and plot the pairwise differences:
Figure 3. Estimate of the progression in Guiraud values for narrative texts by Portuguese-French bilinguals from T1 to T2.
The mean progression is easily calculated:
The estimated error for this estimate is:
And its 95% uncertainty interval is:
According to the model, there’s about a 84% chance that there’s indeed some progression going from T1 to T2.
Question 2: T1-T3 progression for argumentative texts, Portuguese-French vs. Portuguese-German?
This question requires us to take into consideration the modelled average
for argumentative texts written by Portuguese-French bilinguals at T1,
that for argumentative texts written by Portuguese-French bilinguals at T3,
and the same for the texts written by Portuguese-German bilinguals.
We need the following columns in
- 7 (Portuguese-French, T1, argumentative)
- 11 (Portuguese-French, T3, argumentative)
- 13 (Portuguese-German, T1, argumentative)
- 17 (Portuguese-German, T3, argumentative)
We compute the progression for the Portuguese-French bilinguals and that for the Portuguese-German bilinguals. Then we compute the difference between these progressions:
The mean progression for the Portuguese-French bilinguals was 0.54 compared to 0.68 for the Portuguese-German bilinguals:
The mean difference between these progressions, then, is 0.14 in favour of the Portuguese-German bilinguals:
However, there is considerable uncertainty about this difference:
The probability that the Portuguese-German bilinguals make more progress than the Portuguese-French bilinguals is 77%, and according to the model, there’s a 95% chance its size is somewhere between -0.25 and 0.52 points.
By investing some time in recoding your predictors, you can make the parameter estimates more relevant to your questions. Any specific comparisons you may be interested in can additionally be addressed by making use of the entire distribution of estimates. You can also use these estimate distributions to draw effect plots.
- R Library Contrast Coding Systems for categorical variables
- Matthew Kay. 2018. Extracting and visualizing tidy draws from brms models
- Daniel J. Schad, Sven Hohenstein, Shravan Vasishth and Reinhold Kliegl. 2018. How to capitalize on a priori contrasts in linear (mixed) models: A tutorial.
- In simpler models, you can use bootstrapping to generate distributions of estimates. This would be possible for these data, too, but it would take a considerable amount of time.
Audrey Bonvin, Jan Vanhove, Raphael Berthele and Amelia Lambelet. 2018. Die Entwicklung von produktiven lexikalischen Kompetenzen bei Schüler(innen) mit portugiesischem Migrationshintergrund in der Schweiz. Zeitschrift für Interkulturellen Fremdsprachenunterricht 23(1). 135-148. Data and R code available from figshare.
Amelia Lambelet, Raphael Berthele, Magalie Desgrippes, Carlos Pestana and Jan Vanhove. 2017a. Chapter 2: Testing interdependence in Portuguese Heritage speakers in Switzerland: the HELASCOT project. In Raphael Berthele and Amelia Lambelet (eds.), Heritage and school language literacy development in migrant children: Interdependence or independence?, pp. 26-33. Multilingual Matters.
Amelia Lambelet, Magalie Desgrippes and Jan Vanhove. 2017b. Chapter 5: The development of argumentative and narrative writing skills in Portuguese heritage speakers in Switzerland (HELASCOT project). In Raphael Berthele and Amelia Lambelet (eds.), Heritage and school language literacy development in migrant children: Interdependence or independence?, pp. 83-96. Multilingual Matters.
Jan Vanhove, Audrey Bonvin, Amelia Lambelet and Raphael Berthele. Forthcoming. Predicting perceptions of the lexical richness of short French, German, and Portuguese texts. Journal of Writing Research. Technical report, data (including texts), elicitation materials, and R code available from the Open Science Framework.