Multilevel modeling for longitudinal data: concepts and applications

Purpose – This paper aims to discuss multilevel modeling for longitudinal data, clarifying the circumstances in which they can be used. Design/methodology/approach – The authors estimate three-level models with repeated measures, offering conditions for their correct interpretation. Findings – From the concepts and techniques presented, the authors can propose models, in which it is possible to identify the fixed and random effects on the dependent variable, understand the variance decomposition of multilevel random effects, test alternative covariance structures to account for heteroskedasticity and calculate and interpret the intraclass correlations of each analysis level. Originality/value – Understanding how nested data structures and data with repeated measures work enables researchers and managers to define several types of constructs from which multilevel models can be used.


Introduction
Regression models for longitudinal data are very useful when the researcher wishes to study the behavior of a given phenomenon in the presence of nested data structures with repeated, or longitudinal, measures.
While in nested structures of clustered data certain explanatory variables do not present variation between observations (representing a level of analysis) coming from a given group (representing another level of analysis), in data structures with repeated measures there is also the temporal evolution, a fact that enables the researcher to investigate the individual reasons that may lead each of the observations to present different behaviors of the dependent variable, for the same group or for distinct groups, over time (Fávero, 2010;Martins & Terra, 2015;Misangyi, Lepine, Algina, & Goeddeke, 2006).
For example, certain school data that does not vary among students, such as location and size, can be compared with data from other schools; and certain student data, such as sex and religion, that do not vary over time, can be compared with data from other students, which allows the different influences in the dependent variable to be analyzed. In all of these situations (nested data without or with repeated measures), datasets provide structures from which hierarchical models can be estimated.
Multilevel regression models have become considerably important in several fields of knowledge, and the publication of papers that use estimations related to these models has become more and more frequent (Goldstein, 2011). The reason for the importance of multilevel modeling is due mainly to the determination of research constructs that consider the existence of nested data structures, in which certain variables show variation between distinct units that represent groups but do not assess variation between observations that belong to the same group. In addition, the computational development and investments that data analysis software developers have made in the processing capacity to estimate multilevel models have also provided support to researchers who are increasingly interested in this type of approach (Gelman & Hill, 2007;Hough, 2006;Santos, Fávero, & Distadio, 2016;Serra & Fávero, 2018).
Theoretically, researchers can define a construct with a greater number of levels of analysis, even if the interpretation of model parameters is not something trivial. For instance, imagine the study of school performance, throughout time, of students nested into schools, these nested into municipal districts, these into municipalities, and these into states of the federation. In this case, we would be working with six analysis levels (temporal evolution, students, schools, municipal districts, municipalities and states).
The main advantage of multilevel models over traditional regression models estimated, for instance, by ordinary least squares (OLS), is the possibility of considering a natural nesting of data (Steenbergen & Jones, 2002), that is, multilevel models enable us to identify and analyze individual heterogeneities, and heterogeneities between the groups, to which these individuals belong, making it possible to specify random components in each analysis level (Heck & Thomas, 2009).
Multilevel models correct for the fact that observations in the same group are not independent and thus, compared to OLS models, lead to unbiased estimates of standard errors (SEs). But one could say that the same can be obtained with clustered standard errors in OLS. Indeed, if the number of clusters is plentiful (i.e. above 20), clustered SEs in OLS models and multilevel models are equally adequate for precision estimates of group-level effects. On the other hand, if there are less than 20 clusters, researchers should avoid using clustered SEs and adopt multilevel modeling. Furthermore, if researchers are also interested in testing whether group-level covariates moderate individual-level effects, multilevel models may prove to be the most appropriate choice (Arceneaux & Nickerson, 2009;Steenbergen & Jones, 2002).
According to Courgeau (2003), within a model structure with a single equation, there seems to be no connection between individuals and the society in which they live. In this sense, the use of level equations enables the researcher to "jump" from one science to another: students and schools, families and neighborhoods, firms and countries. Ignoring this relationship means to elaborate incorrect analyzes about the behavior of the individuals and, equally, about the behavior of the groups. Only the recognition of these reciprocal influences allows the correct analysis of the phenomena. This is in line with what is called by Mathieu and Chen (2011) the multilevel paradigm, which refers to a way of thinking: considering management phenomena in context and looking for driving variables not only from the focal unit of analysis but also from levels RAUSP 54,4 above and below. Such an approach often implies the development of multidisciplinary theories and investigations, what is the spirit articulated by Hitt, Beamish, Jackson, and Mathieu (2007) when discussing the built of theoretical and empirical bridges across levels through multilevel modeling. Most modern-day multilevel studies seek to associate relationships across proximal layers, such as team attributes and members' attitudes or environmental conditions and variable performance.
In an effort to make multilevel modeling more accessible, we provide the syntax for the mixed procedures in Stata for each step and show how to test and compare these designs in the model-building process. Previous discussions involving multilevel data have illustrated the use of multilevel modeling in programs such as MLn (Kreft & de Leeuw, 1998), R (Bliese & Ployhart, 2002), HLM (Raudenbush, Bryk, Cheong, Congdon, & Du Toit, 2004), and SAS (Littell, Milliken, Stroup, & Wolfinger, 2004;Singer, 1998).
In this paper, our focus will be on hierarchical linear models (HLM), also known as linear mixed models (LMM). According to West, Welch, and Gałecki (2015), the term "linear mixed models" comes from the fact that these models present linear specification and the explanatory variables include a mix of fixed and random effects. That is, they can be inserted into components with fixed effects, as well as into components with random effects. While the estimated fixed effects parameters indicate the relationship between explanatory variables and the metric dependent variable, the random effects components can be represented by the combination of explanatory variables and non-observed random effects.
Our main objectives are: to introduce the concepts of nested data structures; to define the type of model to be estimated based on the characteristics of the data; to estimate parameters through several methods in Stata; to interpret the results obtained through several types of existing estimations for multilevel models; and to define the most suitable estimation for diagnosing and forecasting effects in each of the cases studied.

Nested data structures
Models that take into account the presence of nested structures in the data offer benefits to researchers since they make possible the study of the sources of variance, in different levels, of an outcome variable. Raudenbush and Bryk (2002) discuss the applications of multilevel modeling from nested data structures in various areas of knowledge, particularly education. In this field, works of Aitkin and Longford (1986), Raudenbush and Bryk (1986), Garner and Raudenbush (1991), Raudenbush (1993), Rumberger andThomas (1993), O'Connell andMcCoach (2008), and Goldstein (2011) deserve mention.
Multilevel modeling is also extensively used in strategy literature to compare existing variances at firm and industry levels for firm performance composition, such as Schmalensee (1985), McGahan and Porter (1997), Brush and Bromiley (1997), Mauri andMichaels (1998), Brush, Bromiley, andHendrickx (1999), Chang and Singh (2000), Bowman and Helfat (2001), McGahan and Porter (2002), Ruefli and Wiggins (2003), Short, Ketchen, Palmer, and Hult (2007) and Short, McKelvie, Ketchen, and Chandler (2009). Other authors have analyzed the country-of-origin effect on performance, notably Christmann, Day, and Yip (1999), Lee (2003), Hawawini, Subramanian, and Verdin (2004), Makino, Beamish, and Multilevel modeling for longitudinal data Zhao (2004a), Makino, Isobe, and Chan (2004b), Misangyi et al. (2006), Goldszmidt, Brito, andVasconcelos (2007), Fávero (2008), and Holcomb, Combs, Sirmon, and Sexton (2010). Therefore, multilevel regression models enable us to formally investigate the behavior of a certain dependent variable Y, which represents the phenomenon we are interested in, based on the behavior of explanatory variables, whose changes may occur for clustered data, between observations and between groups to which these observations belong, and for data with repeated measures throughout time. In other words, there must be variables that have data that change between individuals that represent a certain level. But these variables remain unchanged for certain groups of individuals, and these groups represent a higher level.
First, imagine a dataset with data on n individuals, and each individual i = 1, . . ., n belongs to one of the j = 1, . . ., J groups, obviously n > J. Therefore, this dataset can have certain explanatory variables X 1 , . . ., X Q that refer to each individual i, and other explanatory variables W 1 , . . ., W S that refer to each group j; but they are invariable for the individuals of a certain group. Table I shows the general model of a dataset with a two-level clustered/nested data structure (individual and group).
Based on Table I, we can see that X 1 , . . ., X Q (columns 4 to 6) are level-1 variables (data change between individuals), and W 1 , . . ., W S (columns 7 to 9) are level-2 variables (data change between groups; however, not for the individuals in each group). Furthermore, quantities of individuals in groups 1, 2, . . ., J (column 2) are equal, respectively, to n 1 , n 2 À n 1 , . . ., n À n JÀ1 (column 1). Figure 1 shows the existing nesting between the level-1 units (individuals) and the level-2 units (groups), which characterizes the existence of clustered data.
Imagine another dataset in which, in addition to the nesting presented for clustered data, there is temporal evolution. That is, data with repeated measures. Thus, besides the individuals that will now belong to level 2 and therefore will be called j = 1, . . ., J, nested into Table I. General model of a dataset with a twolevel clustered/nested data structure Note: The dataset will have a balanced nested data structure if n 1 = n 2 À n 1 = . . . = n À n JÀ1 RAUSP 54,4 k = 1, . . ., K groups (which now belong to level 3), we will also have t = 1, . . ., T j periods in which each j individual is monitored. Consequently, this new dataset can have the same explanatory variables X 1 , . . ., X Q that refer to each j individual. But now they are invariable for each j individual during the periods of monitoring. Moreover, the dataset can also have the same explanatory variables W 1 , . . ., W S that refer to each group k. But W 1 , . . ., W S are also invariable throughout time for each group k. Table II provides the logic to describe a dataset with a three-level nested data structure with repeated measures (time, individual and group). Based on Table II, we can now see the variable that corresponds to the period is a level-1 explanatory variable (column 1), since the data change is in each row of the dataset, and that X 1 , . . ., X Q (columns 5 to 7) become level-2 variables (data change between individuals, but not for the same individual throughout time), and that W 1 , . . ., W S (columns 8 to 10) become level-3 variables (data change between K groups (column 3), but not for the same group throughout time). Furthermore, quantities of periods in which individuals 1, 2, . . ., J (column 2) are monitored are equal, respectively, to T 1 , T 2 À T 1 , . . ., T J À T JÀ1 (column 1).
Similar to what was shown for the case with two levels, Figure 2 enables us to see the existing nesting between the level-1 units (temporal variation), the level-2 units (individuals), and the level-3 units (groups), which characterizes a data structure with repeated measures.
Through Tables I and II, as well as through the corresponding Figures 1 and 2, we can see that the data structures present absolute nesting. That is, a certain observation can be nested into only one group, and this group into only another higher-level group, and so on.
In the next section, we will estimate multilevel models with repeated measures in Stata, whose econometric development is in Appendix 1. Appendix 2 is intended for the presentation of the commands in Stata.

Estimation of multilevel models with repeated measures in Stata
This section gives researchers the opportunity to estimate multilevel models through Stata Statistical Software®. For our example, we will use the step-up strategic multilevel analysis proposed by Raudenbush and Bryk (2002), and Snijders and Bosker (2011). That is, we first studied the variance decomposition from the definition of a null model (non-conditional model), so that afterwards, a random intercepts model and a random intercepts and slopes model could be estimated. Finally, from the definition of the random nature of the error terms, we estimated the complete model by including level-2 variables into the analysis. We, therefore, estimate a three-level hierarchical linear model, in which the nesting of data will be characterized due to the presence of repeated measures. Thus, there is temporal evolution in the behavior of the dependent variable.  3.1 Hypotheses and data A dataset was constructed by a professor interested in monitoring students' school performance for a certain period of time, in order to investigate if there is variability in this performance throughout time between students within the same school, and between those from different schools. In addition, if yes, if there are certain student and school characteristics that explain this variability. This dataset follows the logic of the seminal work developed by Raudenbush, Rowan, and Kang (1991). A total of 15 schools volunteered to provide data on their students' school performance (scores from 0 to 100) in the last four years, a total of 610 students. In addition, the professor also obtained each student's gender in the dataset in order to verify if there are differences in school performance resulting from this variable. The variable regarding professors' years of teaching experience, for each school, also was included in the study. The dataset PerformanceTimeStudentSchool.dta can be found in Fávero and Belfiore (2019).
It is important to mention that, although traditional maximum likelihood estimation methods for multilevel modeling have been shown to provide biased estimates when the number of clusters is below 30, methods such as restricted maximum likelihood (REML) estimation have shown potential to perform well with ten clusters or fewer (McNeish & Stapleton, 2016). In this paper, as discussed below, we estimate all models through REML.
First, we test three fundamental hypotheses regarding the nature of student's school performance over time, following the logic proposed by Raudenbush and Bryk (2002) and Short, Ketchen, Bennett, and Du Toit (2006): H1. There is significant variance in student's school performance over time, both within students and schools.
H2. Student's school performance follows a linear trend over time.
H3. There is significant variance around a linear performance trend.
Student and schools' characteristics may influence student's school performance over time. Thus, we test the following: H4. There is a significant relationship between students' characteristics and student's school performance over time. H5. There is a significant relationship between schools' characteristics and student's school performance over time.
We have a balanced longitudinal data structure since all 610 students are monitored in the four periods. Figure 3 enables us to analyze the temporal evolution of the school performance of the first 50 students in the sample. From the trends in the lines, we can see that the temporal evolutions of the school performance have different intercepts and slopes between students. These different intercepts and slopes justify the use of multilevel modeling and provide reasons to include intercept and slope random effects in Level 2 of the models that will be estimated. Figure 4 shows the temporal evolutions of the average school performance. The increasing trend over time provides further justification for estimating a three-level   Figure 4 also shows the linear adjustment through OLS of the school performance behavior over time for each school. In addition, the figure enables us to display the intercept and slope random effects in Level 3 of the models that will be estimated, since the temporal evolutions of the school performance present different intercepts and slopes between the schools.

3.2
Step-up strategic multilevel modeling Having characterized the temporal nesting of the students from different schools in the data with repeated measures, we can initially estimate a null model (non-conditional model) that enables us to determine if there is variability in the school performance between students from the same school and between those from different schools. No explanatory variable is inserted into the modeling, which only considers the existence of one intercept and of error terms u 00k , r 0jk and e tjk , with variances, respectively, equal to t u000 , t r000 and s 2 . The model to be estimated has the following expression: Null Model: which results in[1]: At the top of Figure 5, we can initially demonstrate that we have a balanced longitudinal data structure since for each student we have minimum and maximum quantities of periods of monitoring equal to four, with a mean also equal to four. About the fixed effects component, we can see that the estimation of parameter g 000 is equal to 68.714, which corresponds to the average of students' expected annual school performance of the (horizontal line estimated in the null model, or general intercept). Moreover, at the bottom of Figure 5, the estimations of the variances of error terms t u000 = 180.194, t r000 = 325.799 and s 2 = 41.649 are presented.
We can, therefore, define two intraclass correlations, given the existence of two variance proportions. The first one refers to the correlation between the data of variable performance in t and in t 0 (t = t 0 ) of a certain student j from a certain school k (Level-2 intraclass correlation). The other one refers to the correlation between the data of variable performance in t and in t 0 (t = t 0 ) of different students j and j 0 (j = j 0 ) from a certain school k (Level-3 intraclass correlation). Therefore, we have: Level-2 intraclass correlation: Hence, the correlation between annual school performances is equal to 32.9 per cent (rho school ) for the same school, and the correlation between annual school performances is equal to 92.4 per cent (rho student|school ) for the same student of a certain school. Therefore, for the model without explanatory variables, while the annual school performance is slightly correlated between schools, the same becomes strongly correlated when the calculation is carried out for the same student from a certain school. In this last case, we estimate that students and schools random effects representing approximately 92 per cent of the total variance of the residuals[2]. Regarding the statistical significance of these variances, the fact that the estimated values of t u000 , t r000 and s 2 are considerably higher than their respective standard errors suggests that there is significant variation in the annual school performance between students and between schools.
This information is essential to underpin the choice of the multilevel modeling, instead of a simple and traditional regression model through OLS. At the bottom of Figure 5, we can verify this fact by analyzing the result of the likelihood-ratio test. Given Sig. x 2 = 0.000, we can reject the null hypothesis that the random intercepts are equal to zero (H 0 : u 00k = r 0jk = 0), which makes the estimation of a traditional linear regression model be ruled out for the data with repeated measures. Even though researchers frequently ignore the estimation of null models, analyzing the results may help to reject the research hypotheses or not. It may even provide adjustments in relation to the constructs proposed. For our data, the results of the null model allow us to state that there is significant variability in the school performance throughout the four years under analysis. Furthermore, there is significant variability in the school performance over time between students of the same school, and there is significant variability in the school performance over time between students from different schools.
H1 can be supported, i.e. there is significant variance in student's school performance over time within both students and schools. Since our main objective is to verify if there are student and school characteristics that would explain the variability in the school performance between students from the same school and between those from different schools, we will continue with the next modeling steps, respecting the step-up strategic multilevel analysis.
Let us insert level-1 variable year into the analysis, aiming at investigating if the temporal variable has a relationship to students' school performance behavior and, more than this if the school performance has a linear behavior throughout time.
Linear Trend Model with Random Intercepts: which results in the following expression[1]: First, we can see that the mean of the annual increase in school performance is statistically significant, with an estimated parameter of g 100 = 4.348, ceteris paribus.
Thus, we can also support H2, since student's school performance over time statistically follows a linear trend. Regarding the random effects components, we have also verified that there is statistical significance in the variances of u 00k , r 0jk and e tjk , because the estimations of t u000 , t r000 and s 2 are considerably higher than the respective standard errors. Therefore, new intraclass correlations can be calculated, as follows: Level-2 intraclass correlation: Both variance proportions are higher than the ones obtained in the estimation of the null model, which shows the importance of including the variable that corresponds to the repeated measure in level 1. Besides, the result of the likelihood-ratio test at the bottom of Figure 6 allows us to prove that the estimation of a simple traditional linear regression model (performance based on year) only with fixed effects must be ruled out (H3 supported). Therefore, now, our model starts to have the following specification: performance tjk ¼ 57:844 þ 4:348 Á year jk þ u 00k þ r 0jk þ e tjk Figures 7 and 8 provide better visualization of the random intercepts per school and per student. Therefore, we are able to conclude that students' school performance follows a linear trend throughout time. In addition, there is a significant variance of intercepts between those who study at the same school and between those who study at different schools. Thus, we also need to verify if there is a significant variance of the school performance slopes throughout time between the different students. Therefore, let us insert slope random effects into Levels 2 and 3 of our multilevel model that, by maintaining the intercept random effects, will have the following expression: Linear Trend Model with Random Intercepts and Slopes: which results in[1]: performance tjk ¼ g 000 þ g 100 Á year jk þ u 00k þ u 10k Á year jk þ r 0jk þ r 1jk Á year jk þ e tjk Note that the variable year is present in the fixed effects component and in the Level-3 random effects components (by multiplying error term u 10k ), and in the Level-2 ones (by multiplying error term r 1jk ). Figure 9 shows the results obtained through this estimation. We can see that, even though the fixed-effects parameter estimations do not change considerably in relation to the previous model, the variance estimations are different, which generates new intraclass correlations, as follows: Level-2 intraclass correlation: Level-3 intraclass correlation: Therefore, for this model, we estimate that the students and schools random effects represent approximately 99 per cent of the total variance of the residuals. Figure 10 shows a likelihood-ratio test applied to compare the estimations of the linear trend models with random intercepts and with random intercepts and slopes. By using the values of the restricted likelihood functions obtained in Figures 5 and 9, we arrive at the x 2 statistic for the test, with 2 degrees of freedom: x 2 2 = (À2·LL r-randomintercept À (À2·LL r-randomslope )) = {À2·(À7,801.420) À [À2·(À7,464.819)]} = 673.20 which results in a Sig. x 2 2 = 0.000 < 0.05 and ends up favoring the linear trend model with and random intercepts and slopes. It is important to mention that this likelihood-ratio test is only valid when a comparison of the estimations obtained through REML of two models is carried out with identical fixed effects specification. Given that in our case, both models, that were estimated through REML, present the same fixed effects specification, g 000 þ g 100 year jk , the test is considered valid. Hence, our model starts to have the following specification: In the current situation, we are able to state that students' school performance follows a linear trend throughout time. In addition, there is a significant variance of intercepts and slopes between those students who study at the same school and between those who study at different schools. Therefore, let us insert Level-2 variable gender into the analysis to determine if this characteristic explains the variation in the annual school performance between students. Linear Trend Model with Random Intercepts and Slopes and with Level-2 Variable gender: which results in the following expression[1]: performance tjk ¼ g 000 þ g 100 Á year jk þ g 010 Á gender jk þ g 110 Á gender jk Á year jk þ u 00k þu 10k Á year jk þ r 0jk þ r 1jk Á year jk þ e tjk Figure 10.

Multilevel modeling for longitudinal data
This model shows significant estimations for the fixed effects parameters, as well as for the variances of the random effects terms, at a significance level of < 0.05 ( Figure 11). Moreover, at this moment of the modeling, we are able to state that students' school performance follows a linear trend throughout time, and there is a significant variance of intercepts and slopes between those who study at the same school and between those who study at different schools. Additionally, the fact that a certain student is female or male is part of the reason why there is this variation in school performance (H4 supported).
Finally, let us investigate if Level-3 variable texp (professors' years of teaching experience), also explains the variation in the annual school performance between the students. After some intermediate analyses, let us move on to estimate the three-level hierarchical model with the following specification: Linear Trend Model with Random Intercepts and Slopes, Level-2 Variable gender and Level-3 Variable texp (Complete Model): Figure 11. Linear trend model with random intercepts and slopes and Level-2 variable gender RAUSP 54,4 which results in the following expression[1]: Even though the estimations of the fixed effects parameters and random effects variances are significant, at a significance level of 0.05, it is necessary to study the structure of the random effects (u 00k , u 10k and r 0jk , r 1jk ) variance-covariance matrix. Based on the outputs found in Figure 12, we have: Since we did not specify any covariance structure for these error terms, we are assuming that this structure is independent, that is, both cov(u 00k , u 10k ) = 0 and cov(r 0jk , r 1jk ) = 0. Nevertheless, we can generalize the structure of these matrices by allowing u 00k and u 10k to be correlated, and r 0jk and r 1jk to be correlated too. Thus, following Short et al. (2006), one additional contribution of this paper is the testing of alternative covariance structures to account for heteroskedasticity. In our sample, we found significant differences based on the model assumptions and covariance structure specified in our empirical test. Figure 13.
Outputs of the linear trend model with random intercepts and slopes and Level-2 variable gender and Level-3 variable texp, with correlated random effects (u 00k , u 10k ) and (r 0jk , r 1jk ) RAUSP 54,4 The results obtained through the estimation considering correlated random effects (u 00k , u 10k ) and (r 0jk , r 1jk ) are shown in Figure 13. The fixed effects parameter estimations are close to those obtained when estimating the model that considers the existence of a structure that is independent from the random-effects variance-covariance matrices (Figure 12). Regarding the random-effects parameters, except for the estimations of u 10k and cov(u 00k , u 10k ), which are statistically significant at a significance level of 0.10, all the other estimations are significant at a significance level of 0.05. Considering that cov(u 00k , u 10k ) and cov(r 0jk , r 1jk ) are statistically different from zero, based on the outputs in Figure 13 we can write that: Random effects variance-covariance matrix for level school: Even statistically different from zero, the estimations of the random effects covariances in both levels of the analysis, if researchers wish to prove the better suitability of this last model over the one that considers the matrix with independent error terms, they just need to run a likelihood-ratio test to compare both estimations ( Figure 14). Figure 15 presents the result of the likelihood-ratio test applied to compare the estimations of the complete models with independent and correlated random effects (u 00k , u 10k ) and (r 0jk , r 1jk ).

Figure 14.
Variance-covariance matrices with correlated random effects (u 00k , u 10k ) and (r 0jk , r 1jk ) Figure 15. Likelihood-ratio test Multilevel modeling for longitudinal data As Sig. x 2 2 = 0.000 < 0.05, we can state that the structure of the random effects variancecovariance matrices can be considered unstructured. That is, we can consider that error terms u 00k and u 10k are correlated (cov(u 00k , u 10k ) = 0) and that error terms r 0jk and r 1jk are correlated too (cov(r 0jk , r 1jk ) = 0).

Á year jk
We have seen that students' school performance follows a linear trend over time. Moreover, there is a significant variance of intercepts and slopes between those who study at the same school and between those who study at different schools, and students' gender is significant to explain part of this variation. H5 can be supported since professors' years of teaching experience at each school (Level-3 variable) itself also explains part of the discrepancies in the annual school performance between students from different schools. Figure 16 shows the predicted values of school performance throughout time for the first 50 students in the sample (yhatstudent) and, through which, we can see different intercepts and slopes throughout time for different students.
Finally, a more inquisitive researcher, aiming at questioning the superiority of multilevel models in relation to traditional regression models estimated through OLS, whenever there are datasets with nested structures, decides to construct a graph (Figure 17). Through this graph, it is possible to compare the predicted school performance values generated by this three-level hierarchical modeling (HLM3) to those generated by an estimation through OLS, for all the students in the sample, in each of the periods analyzed, using the same explanatory variables year, gender, texp, genderyear and texpyear. Obviously, there are only fixed effects in the estimation through OLS.
The gray line at 45°shows the observed school performance values of each one of the students in the sample in each of the periods analyzed (performance x performance). By RAUSP 54,4 using Figure 17, we can clearly see the superiority of the linear trend model with explanatory variables and random intercepts and slopes in levels 2 and 3 (complete HLM3 model) over the multiple linear regression model estimated through OLS with the same explanatory variables in this case. The absence of random terms for each contextual effect in traditional regression models, such as OLS, prevents greater adherence between predicted and observed values of the outcome variable in phenomena where is possible to directly or indirectly identify hierarchies, or levels, in the data structure (Raudenbush & Bryk, 2002). This demonstrates the importance of considering the random effects components whenever there are nested data structures.

Final remarks
This paper provides a brief discussion about the concepts, processes, stages, tasks, and the types of methods and techniques it can employ. It enables researchers and managers to assess the relationship between a certain performance variable and one or more predictor variables, which characterize different levels of analysis. Moreover, as well as contributed and discussed by Short et al. (2006), this study offers progress toward resolving the ambiguity related to the structure of the random effects variance-covariance matrices by applying a multilevel model with random intercepts and slopes, by explicitly accounting for and modeling heteroskedasticity in the data analysis, and by articulating the importance of and offering interpretations for different specifications to testing the effects of time on performance. Time studies in business management research have not reflected enough attention to these issues, and multilevel modeling provides a tool to ameliorate such issues when using a longitudinal design.
Multilevel modeling is a broad theme that can often be explored in depth in the field of business management. Each level is formed by individuals or groups nested into other groups and so on. Since variables from a certain group are invariable between groups or individuals that correspond to lower levels that are nested into that group, it is natural for many researchers and constructs to use such models (Zhang, Li, & Song, 2014).
Many can be the characteristics of the datasets with nested data structures. The most common are those with absolute nesting, in which there are clustered data or data with repeated measures. In this paper, we chose to use a dataset to estimate three-level hierarchical linear models with repeated measures. Nonetheless, from which, we believe researchers will have the conditions to estimate, for example, for three-level models with clustered data or even to consider a higher number of analysis levels, resulting from more complex nesting structures.
Multilevel models enable us to identify and analyze individual heterogeneities and the heterogeneities between the groups to which these individuals belong, making it possible to specify random components in each analysis level. This fact represents the main difference of the traditional regression models estimated through OLS, which cannot consider the natural nesting of data and, consequently, generate biased parameter estimators (Fávero & Belfiore, 2017;Lazega & Snijders, 2016;Pinheiro & Bates, 2009).
Although many papers use multilevel models only to estimate null models to investigate the variance decomposition of the phenomenon being studied in the different analysis levels, the possibility of including explanatory variables that correspond to the different levels in the fixed and random effects components enables us to investigate possible relationships between these variables and the dependent variable. This makes it possible to establish and examine new research objectives and interesting constructs (Gelman, 2006).
The models we studied in this paper are part of what we call generalized linear latent and mixed models (GLLAMM), which encompass the hierarchical linear models (HLM). In this sense, as discussed by Fávero, Santos, and Serra (2018a), a researcher can also propose hierarchical nonlinear models, that refer to the situations in which, in a nested data structure, the dependent variable is a categorical variable or a variable with count data. Thus, logistic, Poisson or negative binomial multilevel models could also be estimated.
Discovering implicit and contextual standards from larger and larger volumes of data becomes an essential condition for organizations to become successful in competitive environments, and multilevel modeling contributes in a considerable way to our ability to understand phenomena. RAUSP 54,4 Notes 1. All parameters are defined in the Appendix 1.
2. One might notice that the sum of variance percentages is not equal to 1. This is the default for displaying Stata outputs. In order to obtain the variance percentages of each level on the outcome variable, we must proceed with the following calculations: Level-1 intraclass correlation: where g pqs (s = 0, 1, . . ., S pq ) refer to the level-3 coefficients, W sk is the s-th level-3 explanatory variable for unit k, and u pqk are the level-3 random effects, assuming that, for each unit k, the vector formed by terms u pqk follows a multivariate normal distribution with each element having mean zero and variance t up pp , which results in a variance-covariance matrix T b with a maximum dimension equal to: which depends on the number of Level-3 coefficients specified with random effects. Let's imagine a single Level-1 explanatory variable that corresponds to the periods in which the data of the dependent variable are monitored. In other words, Level-2 units j nested into Level-3 units k are monitored for a period t (t = 1, . . ., T j ), which makes the dataset have j time series, as shown in Table II. The main objective is to verify if there are discrepancies in the temporal evolution of the data of the dependent variable and, if yes, if these occur due to characteristics of the Level-2 and Level-3 units. This temporal evolution is what characterizes the term repeated measures.
In this regard, Expression (A1) can be rewritten as follows, in which subscripts i become subscripts t: where p 0jk represents the intercept of the model that corresponds to the temporal evolution of the dependent variable of Level-2 unit j nested into Level-3 unit k, and p 1jk corresponds to the average evolution (slope) of the dependent variable for the same unit throughout the period analyzed. The substructures that correspond to Levels 2 and 3 remain with the same specifications as those respectively presented in Expressions (A2) and (A3). Figure A1 shows the plotting of the set of models represented by Expression (A5) in a conceptual way. Through the plotting of the models, we can see that the individual models that represent Level-2 units j can present different intercepts and slopes throughout period t. The fact that this may occur is due to certain characteristics of the Level-2 units j themselves or due to characteristics of the Level-3 units k. Figure A1. Individual models that represent the temporal evolution of the dependent variable for each of the J Level-2 units Multilevel modeling for longitudinal data Besides that, g 010 represents the increase in the expected value of the dependent variable at the initial moment for a certain unit jk when there is a unit alteration in the characteristic X of j, ceteris paribus, and g 011 represents the increase in the expected value of the dependent variable at the initial moment for a certain unit jk when there is a unit alteration in W.X, ceteris paribus. Moreover, u 00k and u 01k represent the error terms that indicate there is randomness in the intercepts, and the last one impacts the alterations in variable X.
Additionally, g 100 represents the alteration in the expected value of the dependent variable when there is a unit alteration in the analysis period (change in the slope due to a unit temporal evolution), ceteris paribus, g 101 represents the alteration in the expected value of the dependent variable due to a unit temporal evolution for a certain unit jk when there is a unit alteration in the characteristic W, ceteris paribus.
Finally, g 110 represents the alteration in the expected value of the dependent variable due to a unit temporal evolution for a certain unit jk when there is a unit alteration in the characteristic X, ceteris paribus, and g 111 represents the alteration in the expected value of the dependent variable due to a unit temporal evolution for a certain unit jk when there is a unit alteration in W.X, ceteris paribus. Terms u 10k and u 11k represent errors that indicate there is randomness in the slopes, and the last one impacts the alterations in variable X.
Expression (A13) facilitates the visualization that the intercept and slope can be influenced by random effects resulting from different behaviors of the dependent variable throughout time for each of the level-2 units (different time series), and this phenomenon can be a result of these units' characteristics, as well as of characteristics of the groups to which such units belong.
If researchers wish to elaborate an analysis about the fixed and random effects components that can influence the behavior of the dependent variable, given that this procedure even facilitates the insertion of the commands to estimate multilevel models in Stata, as we will see below we just need to rearrange the terms of Expression (A13) as follows: Y tjk ¼ g 000 þ g 001 Á W k þ g 010 Á X jk þ g 011 Á W k Á X jk þ g 100 Á period jk þ g 101 Á W k Á period jk þ g 110 Á X jk Á period jk þ g 111 Á W k Á X jk Á period jk ) Fixed Effects þ u 00k þ u 01k Á X jk þ u 10k Á period jk þ u 11k Á X jk Á period jk þ r 0jk þ r 1jk Á period jk þ e tjk |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Random Effects As multilevel models allow interactions between variables in the fixed effects component and, more than that, allow interactions between error terms and variables in the random effects component, the model can be naturally heteroskedastic. Regarding Expression (A14), if the variances of the random terms u 10k , u 11k , r 0jk and r 1jk are statistically different from zero, traditional parameter estimations, such as OLS, will not be adequate. In three-level hierarchical models, we can define two intraclass correlations given the existence of two variance proportions. One corresponds to the behavior of the data that belong to the same Level-2 units j and the same Level-3 units k (Level-2 intraclass correlation), and the other corresponds to the behavior of the data that belong to the same Level-3 units k. But the data are from different Level-2 units j (Level-3 intraclass correlation).
While fixed effects parameters are estimated traditionally th rough maximum likelihood, the variance components of the error terms can be estimated through maximum likelihood or through REML.

Appendix 2. Commands in stata
We now present all Stata commands (Version 15) used throughout the paper to generate tables and graphs and to estimate models: Figure 3: graph twoway connected performance year if student <= 50, connect(L) Figure 4: statsby intercept=_b [_cons] slope=_b [year], by(school) saving (ols, replace): reg performance year sort school merge school using ols drop _merge gen yhat_ols= intercept 1 slope*year sort school year separate performance, by(school) separate yhat_ols, by(school) graph twoway connected yhat_ols1-yhat_ols15 year || lfit performance year, clwidth(thick) clcolor(black) legend(off) ytitle(performance at school) Figure 5: xtmixed performance || school: || student:, var nolog reml This command shows two random effects components, one that corresponds to level 3 (school) and another to level 2 (student). The order in which the random effects components are inserted into the command xtmixed is decreasing when there are more than two levels. That is, we must begin with the highest data nesting level and continue until the lowest level (level 2): Obtention of intraclass correlation: estat icc right after the estimation of the corresponding model.