A Bayesian Data Science Approach To A Multilevel Problem In Hrm

Literature review

For a state of art, wide and exhaustive introduction to multilevel thinking and modelling in HRM research, we refer to Bliese (2002), Renkema, Meijerink and Bondarouk (2017), Shen et al; (2018), Shen (2016), Aguinis et al; (2013), and the references therein. Gelman and Hill (2007) is a good, perhaps dated from a computational point of view, introduction to general Bayesian statistical methods. Gelman et al. (2013) is a quite complete, but mathematically and technically demanding, presentation of modern Bayesian methods with an appendix on using Stan and R. McElreath (2016) is a gentle introduction to Bayesian methods, including multilevel linear models, with R and Stan. Fox and Weisberg (2019) is a good reference for applied regression with R, including multilevel (mixed) models, with a freely available web appendix on Bayesian estimation with Stan. For the Stan program, the state of art in Bayesian statistical software, we manly refer to Carpenter et al; (2017). Gabry and Goodrich (2018) present the rstanarm package, a Stan based R package. Gabry, Mahr, Bürkner, Modrák, and Barret (2018) is the reference for the bayesplot package for plotting Bayesian models. Muth et al. (2018) is an excellent introduction to Bayesian hierarchical linear models fitting with rstanarm package and the user-friendly exploration of results with the shinystan package. References therein can guide the interested readers towards a deeper knowledge of the subject.

The Bayesian approach to model fitting with Stan and rstanarm.

The statistical models presented later in this paper can all be estimated and fitted to data with a frequentist approach. However when data presenting a multilevel structure asks for hierarchical modelling, the Bayesian statistical approach is the most suitable one. The starting point of Bayesian statistics is Bayes’ theorem. In a regression model setting its general form states that the posterior probability distribution $p\left(\theta \right|y,x)$ of parameters $\theta$ given data $y$ and $x$ is proportional to the product of the likelihood function, $p\left(y\right|\theta ,x)$, of the observed response data $y$ given the parameters $\theta$ and the regressor data $x$, and the prior distribution $p\left(\theta \right)$ of the parameters:

$$p\left(\theta |y,x\right)\propto p\left(y|\theta ,x\right)p\left(\theta \right).$$

Here we have assumed that, as it is the case in standard regression setting, $x$ and $\theta$ are independent, so that the prior $p\left(\theta \right)$ doesn’t depend on the regressor $x$.

The most compelling applications of Bayesian methods typically aren't to standard regression models, such as the normal linear model or generalized linear models, but to problem where we want to specify a customized probability model for the data such as a Hierarchical Linear Model, in short HLM (Fox & Weisberg, 2019; Muth et al;, 2018). HLM is an extension of the classical linear model that is well suited for data like ours that present a nested structure where individuals are grouped in teams. Relationships are expected to be similar in the same group and groups are not generally exchangeable (Gelman & Hill, 2007; Hox, Moerbeek, & van de Schoot, 2018). Modeling in the Bayesian statistical framework is mainly composed of four key steps (Muth et al., 2018): data model and prior specification, parameters estimation, sampling quality and model fit checking, results summarization and interpretation. Here, because of space limitations, we will not deeply explore all the four steps. The Stan estimation engine has been used to fit the statistical models described in the next sections via the rstanarm package. We used the default weakly informative priors of rstanarm. The rstanarm package allows the user to easily specify the models considered (Carpenter et al., 2017). Several fundamental tools for MCMC diagnostic are also implemented in the package as well as useful numerical device for results summarization. Graphical representation tools used in this paper are provided by the bayesplot package (Gabry et al., 2018).

We encourage the reader to further explore the literature on Bayesian methods, for example via Gelman and Hill (2007), Gelman et al. (2013), Kruschke (2015), McElreath (2016).

The data

Our data set is part of a publicly available larger data set originated by Aguinis et al., (2013) and patterned after a study by Chen et al; (2007), composed of individual observations of several variables on 630 employee. We have selected 4 columns from the original data set, named as: quality of leader-member exchange ( LMX ), individual empowerment ( IE ), leadership climate ( LC ), and Team . So our data set is composed of 630 rows (observations) and 4 columns (variables). The 630 employees are nested in 105 teams, identified by the Team variable, each composed of 6 individuals. The leadership climate is a Team dependent variable. The goals of the previously cited studies were to investigate whether the quality of leader-member exchange ( LMX ) predicts individual empowerment ( IE ) and the role of the team dependent leadership climate ( LC ). Chen et al.’s model predicted that individuals who report a better relationship with their leader (i.e., a higher LMX value) has the autonomy and capability to perform meaningful works that positively affect their organization (i.e., they feel more empowered). The model assumed also that the team-level variable LC would affect the individual-level empowerment IE positively. It was also supposed that the relationship between LMX and IE would be stronger for teams with higher LC . See Aguinis et al. (2013) for more information on the data.

The statistical models

The first model considered in this paper is a simple linear model explaining the relationship between the level-1 variables individual empowerment ( IE ) and leader-member exchange ( $LMX$):

$${IE}_i=\mathrm{b}\mathrm{0}+\mathrm{b}\mathrm{1}\mathrm{}{LMX}_i+{\epsilon }_i,\mathrm{}\mathrm{}{\epsilon }_i\sim Normal(\mathrm{0},s^{\mathrm{2}})\mathrm{}.$$

We will refer to this model as model1. This model could be considered a conceptually wrong one, because it ignores the nested (hierarchical) structure of our data. However we will fit it to data to derive an overall relationship between IE and LMX. We will later compare model1 to each team-specific fitted regression line considered in model2 (see Figure 01 ), as defined below.

The second considered model, referred as model2, is a hierarchical two-levels model , with random intercept and random slope:

$${IE}_{ij}=\mathrm{b}\mathrm{0}+\mathrm{b}\mathrm{0}\mathrm{j}+\left(\mathrm{b}\mathrm{1}+\mathrm{b}\mathrm{1}\mathrm{j}\right)\mathrm{}{LMX}_{ij}+{\epsilon }_{ij},$$

$$\left(\begin{array}{c}\mathrm{b}\mathrm{0}\mathrm{j}\\ \mathrm{b}\mathrm{1}\mathrm{j}\end{array}\right)\mathrm{}\sim Normal\left(\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\end{array}\right),\mathrm{}\left(\begin{array}{cc}{\mathrm{s}\mathrm{0}}^{\mathrm{2}}& r\mathrm{}\mathrm{s}\mathrm{0}\mathrm{}\mathrm{s}\mathrm{1}\\ r\mathrm{}\mathrm{s}\mathrm{0}\mathrm{}\mathrm{s}\mathrm{1}& {\mathrm{s}\mathrm{1}}^{\mathrm{2}}\end{array}\right)\right)\mathrm{},\mathrm{}\mathrm{}{\epsilon }_{ij}\sim Normal(\mathrm{0},s^{\mathrm{2}})\mathrm{}.$$

This model improves model1 by taking into account the nested structure of the data: employees are grouped in teams. It is expected that people in the same team share similar IE-LMX relationship. The strength of this relationship is supposed to vary among teams. In model2 , parameters $\mathrm{b}\mathrm{0}$ and $\mathrm{b}\mathrm{1}$ are the overall intercept and slope parameters, while $\mathrm{b}\mathrm{0}+\mathrm{b}\mathrm{0}\mathrm{j}$ and $b1+b1j$ are the intercept and slope parameters for the $j$th team. The error ${\epsilon }_{ij}$ reflect individual-level variation within each team, and its variance $s^2$ quantifies the within-team variation. The variance of $\mathrm{b}\mathrm{0}\mathrm{j}$, denoted by ${\mathrm{s}\mathrm{0}}^2$, quantifies the degree of heterogeneity in intercept across teams. The variance of $\mathrm{b}\mathrm{1}\mathrm{j}$, denoted by ${\mathrm{s}\mathrm{1}}^2$, quantifies the degree of heterogeneity in slope across teams. The parameter $r$ represents the correlation between intercepts and slopes and it should be interpreted in the following way: a positive value for $r$ means that teams with stronger IE-LMX relationship tend to have higher IE level (Aguinis et al., 2013).

The third model, model3, is a cross-level interaction hierarchical model , with random intercept and random slope:

$${IE}_{ij}=\mathrm{b}\mathrm{0}+\mathrm{b}\mathrm{0}\mathrm{j}+\left(\mathrm{b}\mathrm{1}+\mathrm{b}\mathrm{1}\mathrm{j}\right)\mathrm{}{LMX}_{ij}+\mathrm{b}\mathrm{2}\mathrm{}{LC}_j+\mathrm{b}\mathrm{3}\mathrm{}{LMX}_{ij}\times {LC}_j+{\epsilon }_{ij},$$

$$\left(\begin{array}{c}\mathrm{b}\mathrm{0}\mathrm{j}\\ \mathrm{b}\mathrm{1}\mathrm{j}\end{array}\right)\mathrm{}\sim Normal\left(\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\end{array}\right),\mathrm{}\left(\begin{array}{cc}{\mathrm{s}\mathrm{0}}^{\mathrm{2}}& r\mathrm{}\mathrm{s}\mathrm{0}\mathrm{}\mathrm{s}\mathrm{1}\\ r\mathrm{}\mathrm{s}\mathrm{0}\mathrm{}\mathrm{s}\mathrm{1}& {\mathrm{s}\mathrm{1}}^{\mathrm{2}}\end{array}\right)\right)\mathrm{},\mathrm{}\mathrm{}{\epsilon }_{ij}\sim Normal(\mathrm{0},{\sigma }^{\mathrm{2}})\mathrm{}.$$

This is the model considered in Step 4 of Aguinis et al., 2013. Note that model3 presents the same hierarchical structure of model2 , but it is also hypothesized, as in Chen et al. (2007), that the team-level variable LC would affect individual-level IE and the relationship between IE and LMX.

For more information about models as model2 and model3 , readers may refer to Aguinis et al. (2013), Hox et al., (2018), Snijders & Bosker (2012).

The software

We have intensively used the Stan program (Carpenter et al., 2017) via the rstanarm package (Gabry & Goodrich, 2018), which provides an R software interface to Stan programming. The bayesplot R package (Gabry et al., 2018) has complemented the rstanarm tools providing an extensive library of function for plotting posterior draws, visual MCMC diagnostic, as well as graphical posterior predictive checking.

Results for model1 fitting.

The estimated posterior mean, standard deviation and 95% Credible Intervals (CI’s) for the intercept b0, the slope b1 and the residual standard error s for model1 are reported on Table 01 .

Results for model2 fitting.

Results for model2 has 213 more estimated parameters than results for model1. This is because model2 estimates intercept and slope for each of the 105 teams, and it estimates also the level-2 spread in terms of st. dev.s s0, s1, and correlation r. In Table 02 are reported the estimated mean, standard deviation and 95% C.I. for the overall population intercept b0, the overall population slope b1, the level-1 residual standard error s, the between-team intercept standard deviation s0, the between-team slope standard deviation s1, the between-team intercept-slope correlation r. For space saving, Table 02 is truncated to exclude the 210 estimated intercept and slope increments for all teams. One of the best features of Bayesian estimation is that it returns (posterior) samples for estimated parameters, from which we can derive estimated probability density functions.

In Figure 01 , such estimated density are represented for relevant parameters of model2 , with 80% CI’s. Note that all CI’s don’t contain 0. This can be interpreted in a Bayesian framework as evidence that all parameters should be considered significantly different from 0.

In Figure 02 a graphical comparison of model1 and model2 is represented for the first 4 teams. Dots are the observed data. The dashed line is the estimated posterior regression line from model1. The dark lines are the estimated posterior regression lines from model2 . The light lines represent uncertainty in posterior intercept and slope parameters for each team.

See Full Size >

Results for model3 fitting.

Results for model3 comprise 2 more estimated parameters than results for model2. This is because model3 estimates also coefficients for leadership climate and for leader-member exchange:leadership climate interaction; The estimated mean, standard deviation and 95% C.I. for the relevant parameters of model3 are reported in Table 03 . For space saving, Table 03 is truncated to exclude the 210 estimated intercepts and slopes for all teams. In Figure 03 are represented the posterior medians and 80% and 90% credible intervals for coefficients of LMX, LC, and LMX:LC interaction. All intervals don’t contain 0. This means that model3 should be considered a better model than model2.

Results of a more analytic model comparison are presented in the next subsection.

See Full Size >

Model comparison.

We have evaluated our models with the leave-one-out (loo) cross-validation and compared them on the basis of the difference in expected log predictive density ( elpd-diff). When comparing model1 and model2, we have found elpd-diff=35.2 with an estimated standard error st.err.= 7.3. The positive difference in elpd (more than twice the estimated standard error) indicates that model2 is expected to have a better predictive performance than model1. When comparing model2 and model3, we have found elpd-diff=9.4 with an estimated standard error st.err.= 5.0. It follows that model3 should be preferred to model2 for the predictive performance. For computation we have used the loo R package. For more information on the loo package and leave-one-out cross-validation and model comparison please refer to Vehtari, Gelman and Gabry (2017).

See Full Size >