Multilevel modelling ( American Psychological Association ):-

Multilevel models (MLMs, also known as linear mixed models, hierarchical linear models or mixed-effect models) have become increasingly popular in psychology for analyzing data with repeated measurements or data organized in nested levels (e.g., students in classrooms). Although mathematically sophisticated, MLMs are easy to use once familiar with some basic concepts. In this primer, I use a simple example to demonstrate the fundamental principles of MLMs.

A psychologist recruits three bilingual participants (Adam, Betty and Carol) who are equally fluent in English and Spanish. She shows them three color words printed in a nonmatching color (e.g., the word “red” printed in green) and asks the participants to name the color of the ink instead of the word presented. The color words are always presented in English. However, participants are cued to respond in English for three trials and Spanish for the other three. The psychologist carefully records the time for the participant to begin speaking in milliseconds (ms). She is interested in whether a matching or mismatching response language will affect the time it takes to begin the utterance. The data is presented to the right. Note that the data is in “long” format, with one observation per row (i.e., no averaging of data).

In ordinary least squares (OLS) regression, we would model this data using the following formula: 

Timei ~ B0 + B1 * Languagei + ei

For each case, the “Time” score can be separated into fixed and random effects. The fixed effects include the intercept (B0) and the slope (B1) for the dichotomous independent variable "Language." These are considered fixed because they take on a predetermined set of values. In the presented example, we find an estimated intercept of 622.11 and a (nonsignificant) slope of -24.89, and these numbers apply to all cases. The general error term (ei) captures the error in prediction for case i. These are considered “random” because we assume there is some idiosyncratic variation in each case that cannot be accounted for by the fixed effects. We estimate the variance of this error term to find our standard error.

The issue with this approach is that we repeatedly sampled from the same individual, and this violates the OLS assumption that observations are independent from each other. Instead, observations from the same individual are likely to be “clustered,” and a MLM adds some extra parameters that control for this clustering. In particular, MLMs add structure to the error term from OLS. Instead of one general random-effect that captures how each observation deviates from the predicted fixed-effects, there will be multiple random-effects that capture how observations deviate within a cluster, and how each cluster deviates from the overall group. We estimate the variability for each random-effect and use that to control for the variance when estimating the significance our fixed-effects. Thus, we can model our data at the observation level (micro-level) and at the cluster level (macro-level). This combination of different “levels” of analysis gives rise to the term multi-level modeling.

Returning to our example, we are going to add a random-effect for “Subject” that characterizes each participant’s idiosyncratic variation from the fixed-effect estimates. We may expect each participant has some individual differences that affects the intercept but not the slope. In our example, this would indicate variability in the time participants took to begin speaking in English, but a constant relationship between the response language and the time. This is known as a random-intercepts model, because we are going to estimate a different intercept for each subject. Our regression equation now looks like: 

Timei = B0 + B1 * Languagei + (B0|Subject) + ei

The expression (B0│Subject) means to estimate a different intercept for each participant within the “Subject” variable. While this accounts for the individual by-subject variation, we still have our general error term ei to capture the random deviations within the same subject. In other words, this formula indicates that there are multiple responses per subject, and that each response will depend on that subject’s baseline. Our model will control for this within-cluster variance and test the fixed-effect estimates against the remaining between-cluster variance. In this way, we resolve the issue of having multiple responses per subject. Note, we would not use a MLM for this exact data in practice as estimates of random effects with only three groups are very unreliable. In the left-hand portion of the figure below, we can see a graphical representation of this random-intercepts model. A different intercept is estimated for each participant (dotted lines), assuming the same slope for all participants. In addition, there is also the fixed-effect regression (solid line) that captures the overall group effect.

We may also have reason to believe that the relationship between time to speak and the response language might differ for each individual, in addition to the baseline differences. That is, we might also wish to include random slopes in addition to the random intercepts. We can easily update our model by adding a term for the regression slope: 

Timei = B0 + B1 * Languagei + (B0 + B1 * Language|Subject) + ei

The expression (B0 + B1 * Language|Subject) tells the model to estimate different intercepts and slopes for each individual, as shown in the right-hand portion of the figure below.

To assess whether the addition of the random-slopes improved the fit of our model, we can use a goodness of fit test known as a likelihood ratio (LR) chi-square difference test (i.e., nested model test; Snijders & Bosker, 2012). LR tests are used to compare nested models, wherein a “base model” is compared to another model with additional parameters of interest. Thus, we will compare whether the addition of the random slopes to the random-intercepts model results in an improved fit. Here is a snippet of what the output would look like in R:

Random Intercepts: Timei = B0 + B1 * Languagei + (B0|Subject) + ei

Random Slopes: Timei = B0 + B1 * Languagei + (B0 + B1 * Language|Subject) + ei

The result is nonsignificant, suggesting that addition of the random slopes did not improve the fit of the model. Thus, the random-intercept model appears to be the best fit.

Now that we specified our random effects, we can test whether the fixed effect of language is significant. This is accomplished using an LR test too. To do so, we create a base model that is identical to the original model minus the effect of interest. Thus, we will create a model with only an intercept as a fixed effect, and compare it to our random-intercepts model.

Base Model: Timei = B0 + (B0|Subject) +ei

Random Intercepts: Timei = B0 + B1 * Languagei + (B0|Subject) +ei

📷

The result is significant ((X(1) = 3.96, p < .05), suggesting that the addition of the fixed-effect of language improves model fit.

With the fixed-effects and random-effects specified, we can interpret the fixed-effects similarly to an OLS regression. We find an intercept of 522.11 and a slope of -24.89, suggesting that participants took about 522 ms to begin speaking when responding in English, but began speaking 25 ms faster when responding in Spanish. Observant readers will recognize that this is the same numbers we obtained OLS regression. This is not a mistake. Violating the assumption of independence does not bias coefficient estimates, but does tend to underestimate standard errors. Thus, while this 25 ms difference was not significant using OLS regression, it is significant in a MLM. By controlling for the variability in the time to speak due to each subject’s different baselines (i.e., random-intercepts), the remaining variability indicated a significant fixed-effect of language across individuals.

From this simple example, we can easily scale up our model. For example, since we only used three color words out of the entire population of possible color words, we might use the items as another random-effect to control for “by-item” variability. Our model would then include a (B0│Subject)+(B0 |Item) for random intercepts by-subject and by-item. We could also include random slopes for one or both random-effects. These are called crossed-random effects because the random effects are orthogonal to each other. We can also have nested random effects. For example, if we measured students nested within classrooms, our random-intercept would take the form (B0 |Classrooms/Students), which indicates to the model to estimate a random-intercept for each student and each classroom. We can include more predictors at both the micro-level and macro-level, which can be used as random slopes, and can even interact (a cross-level interaction). This flexibility makes MLM an attractive option for researchers.

The purpose of this primer was to demonstrate some of the fundamentals of MLMs and to encourage researchers to consider them further. Thus, I highly encourage you to continue learning about MLMs. Introductions to using MLMs are available for R (Winters, 2013), SPSS (Hayes, 2006), and SAS (Peugh, 2010). Baayen et al. (2008) and West et al. (2011) provide excellent in-depth discussions. Finally, I highly recommend the textbook chapter on multilevel models in Cohen, Cohen, West and Aiken (2013) and the entire textbook by Gelman and Hill (2006). Happy modeling.

References

Baayen, R.H., Davidson, D.J., & Bates, D.M. (2008). Mixed-effects modeling with crossed random effects for subjects and items. Journal of Memory and Language, 59(4), 390-412.

Cohen, J., Cohen, P., West, S.G., & Aiken, L.S. (2013). Applied multiple regression/correlation analysis for the behavioral sciences (3rd Ed.). New York: Routledge.

Gelman, A., & Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models. Cambridge, England: Cambridge University Press.

Hayes, A.F. (2006). A primer on multilevel modeling. Human Communication Research, 32(4), 385-410.

Peugh, J.L. (2010). A practical guide to multilevel modeling. Journal of School Psychology, 48(1), 85-112.

Snijders, T.A.B., & Bosker, R.J. (2012). Multilevel analysis: An introduction to basic and advanced multilevel modeling (2nd Ed.). London: Sage.

West, S.G., Ryu, E., Kwok, O.M., & Cham, H. (2011). Multilevel modeling: Current and future applications in personality research. Journal of Personality, 79(1), 2-50.

Winter, B. (2013). Linear models and linear mixed effects models in R with linguistic applications. arXiv:1308.5499. Retrieved from http://arxiv.org/pdf/1308.5499.pdf (PDF, 2.69MB).

Hello Marieke,

1) You should first run a model without predictors in order to consider if ICC is different from 0 (otherwise you may consider running a one-level regression analysis) (see Sommet and Morselli, 2017). This should also help you to decide the hierarchical structure of your data.

2) A stepwise regression model can be appropriate especially in case of linear regressions. AIC and BIC values allow you to evaluate step by step how well fit a model.

I hope it helps.

Best regards,

Mattia

I find it useful to add variables one at a time and in blocks to see what happens to the badness of fit and the size of the estimated relations - I try to do that in a way so as to answer meaningful questions - the final chapter of this sets out a modelling strategy

Book Developing multilevel models for analysing contextuality, he...

If you leave out school but include pupils in classes the unexplained between school variance will accrue to the highest included level that is classes. But this may be misleading if it is really the schools that are different and not the classes. Fitting a three-level model will reveal what is going on, and then you may find that you do not need the extra complexity of level 3.

This site gives resources for the three-level model

http://www.bristol.ac.uk/cmm/software/mlwin/mlwin-resources.html#hiernormal

Hello Marieke,

I agree with Mattia. It can be helpful to first calculate a random intercept model (aka 'fully unconditional model') where you do not specify any fixed factors. This will give you the plain variance estimates which you then can use to calculate the intra-class correlation (ICC) for each level. The ICC basically tells you how much of the total variance in your model results from a certain level. On this basis, you can start specifying more sophisticated but also more economical models.

Best,

Steffen

Hello Marieke,

In the book Mixed Effects Models and Extensions in Ecology with R (page 121), the authors suggests a strategy called "top-down" for model building. This strategy is also recommended by Diggle et al (2002). It has four main steps:

1. Start with a model where the fixed component contains all explanatory variables and as many interactions as possible. This is called the beyond optimal model.

2. Using the beyond optimal model, find the optimal structure of the random component. This way, hopefully, the random component will not contain any information that we would like to have in the fixed component. In this case we have to use the REML fit.

3. Once the optimal random structure has been found, it is time to find the optimal fixed structure, using ML.

4. Present the final model using REML estimation.

The book then continues with examples of the "bad approach" and the "good approach" to multilevel model building.

In this case you can apply all variable together and then use point scale and analyse in SPSS for psychometric analysis it will be very helpful.

Mattia Vacchiano : Thanks for the clarification on the ICC. I have a question concerning the ICC, how do you know when the ICC is 'too low' to include it as a level in the model? If the ICC is almost zero but you still include it as a level in the model, it can change results drastically. How can this happen?

Hester

an ICC is not a level but the ratio of variance between units to variance between units plus variance within units) that is variance between units divided by total unexplained variance.

For example with pupils (level 1) nested in schools (level 2) some 10% of the variability in pupil achievement typically lies at the school level (ICC = 0.10). It is a measure of how similar pupils in a school are on the outcome variable. So if ICC = 0.10 and we pick hundreds of pairs of pupils, where a pair is defined as belong to the same school, the correlation will be 0.10. So the ICC is a measure of dependency and the further the ICC goes from 0 (i.e. no dependency) then the greater the correlation/ dependency. That means you have not got as much information as you think you have got from an observed sample of say 1000 pupils in 20 schools. As this effective sample size affects the standard error of estimates , there is an increased risk of a Type 1 error ; that is overstatement of significance.

Kish's (1965) approximate formula for an ICC of 0.10 with 25 pupils in each of 50 schools (a balanced sample) is that there are not 25*50 equals 1250 pupils (that is ICC = 0) but only 368 - a really substantially smaller sample than it would first appear.

Article What Is an Intracluster Correlation Coefficient? Crucial Con...

The multilevel model both aims to estimate the ICC and automatically corrects the estimated standard errors to give better estimates of uncertainty. But of course thus is not just a technical matter for we should be interested substantively in how much a difference a school makes.

see the arguments in

Article Multilevel approaches to modelling contexuality: from nuisan...

Hello Hester,

As stated by Kelvyn Jones, ICC is not a level but a measure to understand if there is (or not) a group-effect. That is, through the ICC we want to understand if observations at LEVEL 1 behave more or less homogeneously depending on the group they belong, that we call LEVEL 2.

To the best of my knowledge there is not an ICC value of reference in order to consider it relevant or not. The further the ICC goes from 0 than higher (and more relevant) will be the group-effect. If the ICC is almost 0 then you could consider your LEVEL 2 has no effect on observations and so that would be a reason not to perform a multilevel analysis.

I hope it helps.

Best,

Mattia

On adding/removing/shrinking variables from your model, if you can it is preferred to add/remove/shrink the variables based on your knowledge of how the variables relate. If you don't know how the variables will likely relate then you can build the model based on the sample data and for multilevel models probably use something based on the likelihood (like information criteria) or cross-validation to help determine which to include. I including the "shrink" verb above because in the linear model selection problem people often shrink the model coefficients using, for example, ridge or lasso regression, and there are packages for this for the multilevel situation (https://cran.r-project.org/web/packages/lmmlasso/lmmlasso.pdf). The importance of model selection is related to the covariance of the predictors. It is worth examining these first to think is any should be transformed (e.g., per capita variables, change variables).

I would prefer step wise

The problem with automated stepwise based on some decision rule involving p values (or the like) is that there may be several theoretically interesting models that will be ignored ; there is often not a single best model but several alternatives that can be insightful; and these approaches often capitalize on chance results - and we all now know of the diffuclity of repilicating results.

https://en.wikipedia.org/wiki/All_models_are_wrong

Well it does depend on your variables as well. I would rather suggest to do it stepwise which would be more easy to analyze as well.

It doesn't seem like either stepwise or putting everything into the model at once are particularly easy to analyze. Muhammad, you are advocating a stepwise approach. Do you recommend the questioner does anything with the p-values output for the coefficients?

If you are still interested in this question, stepwise regression is very unstable. See the attachment and its references. Best, David Booth

I have a hierarchical dataset with observations on teachers within schools. Observations on teachers are repeated (one observation per subject they teach). So I am modelling this as a three-level logit model using melogit in stata (only random intercept). I have two explanatory variables, X1 and X2, both at the teacher level. I estimate four specifications of this model:

1. unconditional

2. with X1 only

3. with X2 only

4. with X1 and X2

In each of these models the estimate of the variance at the school level is greater than zero.

However the estimate of the variance at the teacher level is greater than zero only in the first three specifications. Its estimate is zero in the fourth specification (the standard error is also 0).

How can one explain the change in the variance to 0? What causes the variance to drop to 0?

Thanks

I suggest two further unconditional multilevel models, one with x1 as the response, and separately x2 as the response. I expect that your lower level variables will have a higher level component, and that could account for reduction of the school variation in the original response to zero.

Thanks Kelvin. I am encountering convergence issues when implementing your suggestion. X1and X2 are both nominal variables with more than three categories. Aa melogit is not appropriate I used the gsem module in stata to estimate.the unconditional models. Also the problem I described has the variance at the lower level (teachers) which reduces to absolute 0.

Back to your Chandra Shah original question, what are your levels and how many observations at each level.

Hi Kelvyn

There are 722 school clusters and 5137 teachers. Some teachers have repeated observations, between 1 and 3. Total number of observations is 6834 which includes the repeated observations for some teachers. I am happy to share the actual data file and the Stata code, if required.

------------------------------------------------------------

| No. of Observations per Group

Group Variable | Groups Minimum Average Maximum

----------------+--------------------------------------------

cntschid | 722 1 9.5 24

cnttchid | 5,137 1 1.3 3

Hi Kelvyn

The correlation structure between the variables is as below. The dependent variable is of.

| nqual nsubje~s of

-------------+---------------------------

nqual | 1.0000

nsubjects | 0.4077 1.0000

of | -0.1482 0.1997 1.0000

Hi Kelvyn

Could the problem with my model be that the repeated measures are on some teachers (level 2) and not all? For more than third of all teachers there is a single measure.

Could it could be that the random effect at level 2 is not identified and Stata outputs the result with 0 rather than alerting that there is a lack of identification?

Well, it depends on your variables too. I prefer to suggest doing this step by step, which would be easier to analyze as well.