Chapter 1 · Multilevel Data, Multilevel Questions

Suppose a researcher wants to know whether pupils who spend more time on mathematics homework obtain higher mathematics scores. The data come from a national survey in which schools were sampled first and pupils were then sampled within each school. The researcher runs an ordinary regression of mathematics score on homework time using all pupils, obtains a large and highly significant slope, and writes it up. A reviewer returns the paper with three questions. Are pupils in the same school really independent observations? Is the slope a statement about pupils, about schools, or a mixture of both? And could the answer differ from school to school? This chapter is about why these three questions arise whenever data have a nested structure, and why they are not technicalities but part of the substance of the research question.

Learning objectives

After working through this chapter, you should be able to:

recognize nested, longitudinal, cross-classified, and multiple-membership data structures and name their levels and units;
state at which level each variable in a study is measured, and draw the structure with plates;
explain why observations in the same group tend to resemble each other, and define the intraclass correlation;
explain how ignoring dependence distorts standard errors and significance tests, and compute a design effect;
distinguish within-group, between-group, and total relations, and recognize the ecological and atomistic fallacies;
formulate a research question as a micro-level, macro-level, macro-to-micro, cross-level, or micro-to-macro proposition;
describe what multilevel models offer and what they cost.

1.1What Multilevel Analysis Is#

Snijders & Bosker (2012) describe multilevel analysis as a methodology for data whose variability has a complex pattern, with a particular focus on sources of variability that are nested within one another. Pupils vary around the mean of their class, classes vary around the mean of their school, and schools vary around the mean of the population. Each of these sources of variation can be of interest in its own right, and each can be related to explanatory variables measured at the corresponding level. A multilevel model is a statistical model that represents all of these sources of variation at once, together with the explanatory variables that are thought to account for them.

Multilevel models were developed in several disciplines at roughly the same time, and they carry several names (Table 1.1). In education and psychology the most common names are hierarchical linear model (Raudenbush & Bryk, 2002) and multilevel model (Goldstein, 2011; Hox et al., 2018; Snijders & Bosker, 2012). In statistics and biostatistics the same models are called linear mixed models or mixed-effects models (Laird & Ware, 1982; Verbeke & Molenberghs, 2000), because they combine fixed effects, which are constant across groups, with random effects, which vary across groups. Other names are random coefficient models (Longford, 1993), variance components models (Searle et al., 1992), growth curve models for longitudinal data (Singer & Willett, 2003), and hierarchical Bayesian models (Gelman & Hill, 2007). The names differ in what they emphasize, but for continuous outcomes they refer to one family of models; the hierarchical Bayesian version adds prior distributions for the parameters and is estimated by Bayesian methods (Chapter 13). The generalized linear mixed model extends the same family to binary, ordinal, and count outcomes and is the subject of Part VI of this book.

Table 1.1 Names Under Which Multilevel Models Appear

Name	What the name emphasizes	Typical fields
Hierarchical linear model (HLM)	Separate equations for each level of a nested structure	Education, psychology
Multilevel model	Variables and variation at several levels	Education, sociology, psychology
Linear mixed (effects) model	A single equation with fixed and random effects	Statistics, biostatistics, linguistics
Random coefficient model	Regression coefficients that vary over groups	Statistics, economics
Variance components model	Partitioning outcome variance into level-specific parts	Statistics, genetics, measurement
Growth curve model	Person-specific trajectories over time	Developmental and clinical psychology
Hierarchical Bayesian model	Parameters that follow a distribution with its own parameters	Statistics, cognitive science
Generalized linear mixed model (GLMM)	Random effects in models for binary, ordinal, or count outcomes	All of the above

Note. The names overlap; for continuous outcomes they refer to one family of models written in different ways (see Chapter 2).

This book uses the terms multilevel model, hierarchical linear model, and linear mixed model interchangeably, and it chooses among them according to which way of writing the model is most helpful at a given point. Chapter 2 shows that one model can be written in three equivalent ways (as a set of level-specific equations, as a single mixed-model equation, and as a marginal model for the mean and covariance of the outcome), and that moving among the three ways is the key to reading the output of any software package.

1.2Hierarchies Everywhere#

1.2.1Examples of nested data#

Data in psychology, education, health, and the social sciences are very often grouped. Pupils are taught in classes, and classes belong to schools. Patients are treated by therapists, who work in clinics. Employees work in teams within firms. Survey respondents are interviewed by interviewers and live in households within sampling areas. Table 1.2 lists examples from several fields. In each row, the units in one column are nested within the units in the next column, except where the table marks a classification as crossed.

Table 1.2 Examples of Hierarchical Structures in Psychology and the Social Sciences

Field	Level 1	Level 2	Level 3 (if any)
Education	Pupils	Classes	Schools
Developmental psychology	Children	Families	Neighborhoods
Clinical psychology	Patients	Therapists	Clinics
Organizational psychology	Employees	Work teams	Firms
Survey research	Respondents	Interviewers or households	Sampling areas
Longitudinal research	Measurement occasions	Persons	Groups (e.g., schools)
Experience sampling	Daily or momentary reports	Persons	Couples or families
Experimental psychology	Trials	Participants (crossed with stimuli)	Laboratories
Measurement	Item responses	Persons (crossed with items)	Schools or countries
Meta-analysis	Study participants	Studies	Research groups
Cross-national research	Individuals	Country-by-year surveys	Countries

Note. Each row lists units from the lowest level (level 1) to higher levels. “Crossed with” marks a second classification that is not nested in the first (see Figure 1.1).

Two entries in the table deserve comment. First, repeated measurements on the same persons also form a hierarchy: measurement occasions (level 1) are nested within persons (level 2). Longitudinal and repeated-measures data are therefore multilevel data, and much of the modern analysis of change is multilevel analysis (Hoffman, 2015; Singer & Willett, 2003). Second, some structures are not strictly nested. In a reaction-time experiment in which every participant responds to every stimulus, trials are classified both by participant and by stimulus, and neither classification is nested in the other (Baayen et al., 2008; Judd et al., 2012).

1.2.2Five data structures#

Figure 1.1 draws five common structures as unit diagrams, in which each box or circle is one unit and each line means “belongs to” (see also Browne et al., 2001). In a two-level structure (panel a) every pupil belongs to exactly one school. In a three-level structure (panel b) pupils belong to classes and classes belong to schools. In a repeated-measures or longitudinal structure (panel c) occasions belong to persons, and persons may be observed on different numbers of occasions. In a cross-classified structure (panel d) pupils belong to a neighborhood and to a school, and because pupils from one neighborhood attend different schools and each school draws pupils from several neighborhoods, neither classification is nested in the other (Rasbash & Goldstein, 1994). In a multiple-membership structure (panel e) some pupils belong to more than one school, for example because they changed schools during the year, and each membership receives a weight.

Figure 1.1: Five Multilevel Data Structures Drawn as Unit Diagrams — Figure 1.1 Five Multilevel Data Structures Drawn as Unit Diagrams

Parts II to V of this book deal mainly with the nested structures in panels (a) to (c). Cross-classified and multiple-membership structures are treated in Chapter 26, and crossed designs with persons and items return in the measurement chapter (Chapter 30).

1.2.3Levels, units, and terminology#

The level of a unit is its position in the hierarchy. By convention the lowest level, where the outcome is usually measured, is level 1. Units at level 1 are also called micro-level units, elementary units, or individuals; units at level 2 are called macro-level units, clusters, or groups. This book mostly uses level-1 unit and cluster (or the substantive name, such as pupil and school), and uses group as a neutral word when the kind of cluster does not matter. The number of level-2 units is written $J$ , the number of level-1 units in cluster $j$ is written $n_{j}$ , and the total number of level-1 units is $N = \sum_{j} n_{j}$ .

A variable also has a level: the level of the units on which it is measured. Mathematics score and homework time are pupil-level (level-1) variables. School sector (public or private) and pupil–teacher ratio are school-level (level-2) variables. A school-level variable can also be constructed from pupil data by aggregation; the mean socioeconomic status of the pupils in a school is a school-level variable, even though it is computed from pupil-level scores. Conversely, a school-level variable is often attached to every pupil in the data file (by disaggregation), so the same value appears on many rows. A column that varies within clusters is a level-1 variable, but a column that is constant within clusters may be either a level-2 variable or a level-1 variable that happens not to vary in these data; the level of a variable follows from the design and from the definition of the variable, not from the pattern of values in the file.

1.2.4Where variables live: a first look at plates#

This book draws every model with one graphical grammar, the Unified Model Graph (UMG; Yu, 2026a, 2026b), which is introduced fully in Chapter 2. One element of that grammar is useful from the start. A plate is a rounded box that means “repeated for every unit at this level” (Buntine, 1994). Plates are nested when units are nested. A variable is drawn inside every plate whose index it carries and outside every other plate.

Figure 1.2: Where Variables Live: Three Data Structures Drawn With Plates — Figure 1.2 Where Variables Live: Three Data Structures Drawn With Plates

Figure 1.2 shows three examples. In panel (a) the pupil plate ( $i = 1, \dots, n_{j}$ ) sits inside the school plate ( $j = 1, \dots, J$ ). The pupil variables $x_{ij}$ and $y_{ij}$ are inside the pupil plate because they take a different value for every pupil; the school variable $z_{j}$ is inside the school plate but outside the pupil plate because it takes one value per school. In panel (b) occasions ( $t$ ) are nested in persons ( $i$ ): age or time $a_{t i}$ and the outcome $y_{t i}$ change from occasion to occasion, whereas a stable characteristic such as sex, $w_{i}$ , does not. In panel (c) persons and items are crossed: the response $y_{j k}$ lies where the person plate and the item plate overlap, because there is one response for each combination of person and item, while person variables such as $w_{j}$ and item variables such as $d_{k}$ lie in only one of the two plates. The subscripts in the equations and the plates in the figure carry the same information. A reader who can place every variable of a study in the correct plate has already understood most of what distinguishes a multilevel analysis from an ordinary one.

1.3Why Members of a Group Resemble Each Other#

Children with the same biological parents tend to be more alike than children chosen at random from the population. They are more alike because they share genes, because they share an environment, and usually because of both. The same is true, for different reasons, of pupils in the same school. Pupils in one school share the same school environment and often the same teachers; they influence each other through direct contact and shared norms; and they tend to come from the same neighborhood, so they were already similar before they entered the school. In general, members of a group resemble each other for four kinds of reasons: they share a common context (school resources, a therapist’s skill, a manager’s style); they were selected into the group in ways that depend on their characteristics (families choose neighborhoods and schools; patients choose clinics); they interact and influence one another (contagion, norms, competition); and they are measured in a common way (the same interviewer, the same rater, the same test administration).

The consequence is the same whatever the mechanism: observations in the same group are statistically dependent. Knowing the score of one pupil tells us something about the likely score of another pupil in the same school, beyond what the population mean tells us. Table 1.3 makes the idea concrete with a simple measurement, blood pressure, taken under four designs. As the design moves from unrelated persons (Design I) to repeated measurements on members of many families (Design IV), new sources of variability are added, and each new source that several measurements share introduces a new kind of dependence among the measurements.

Table 1.3 Sources of Variability in Four Designs for Measuring Blood Pressure

Design	Measurements of blood pressure taken …	Sources of variability
I	once on each of many unrelated persons	Measurement error; differences between persons
II	once on each member of many families	Measurement error; differences between members; differences between families
III	repeatedly on each person, over time or conditions	Measurement error; serial (time-related) variation; differences between persons
IV	repeatedly on each member of many families	Measurement error; serial variation; differences between members; differences between families

Note. Designs II and III each add a source of variability to Design I, and Design IV combines both additions. Every added source that is shared by several measurements creates dependence among those measurements.

Dependence is sometimes treated as a nuisance: a violation of the independence assumption that must be “corrected.” It is more useful to treat it as information. The amount of variation that lies between schools, and the way it changes when explanatory variables are added, tells us how much outcomes differ between schools and which variables account for those differences. The variation between persons in repeated-measures data tells us how stable individual differences are. In many studies the dependence is itself the phenomenon of interest.

1.4How Multilevel Data Arise#

The structure of a data set follows from the way the data were collected. Two broad routes lead to multilevel data: observational studies with multistage sampling, and experiments in which groups are formed or assigned.

1.4.1Observational studies and multistage sampling#

Large surveys rarely draw a simple random sample of individuals. It is much cheaper to draw a random sample of schools, then a random sample of classes within each sampled school, and then pupils within each sampled class. Such multistage (or cluster) sampling is cost-effective, but it produces data in which the primary sampling units are clusters (Kish, 1965). The usual assumptions of a multilevel analysis mirror the sampling design: the level-2 units are a random (or at least representative) sample from a population of level-2 units, and within each sampled level-2 unit the level-1 units are a random (or representative) sample from that unit. When these assumptions hold approximately, the variation among sampled schools estimates the variation among schools in the population, and inferences can be made about schools in general and not only about the schools in the sample.

Many of the data sets used in this book come from designs of this kind. The National Education Longitudinal Study of 1988 (NELS:88) sampled schools and then eighth-grade students within schools; the High School and Beyond survey used in Chapter 5 has the same structure; international assessments such as PISA sample schools and then students within schools, and TIMSS samples schools and then intact classes within schools. In each case the analysis must respect the stages of the sample.

Taiwan case Taiwan Education Panel Survey

The Taiwan Education Panel Survey (TEPS), distributed by the Survey Research Data Archive (SRDA) of Academia Sinica, began in 2001 with two cohorts: seventh graders, and eleventh graders in senior high, vocational high, and five-year junior college programs. Its sampling design selected schools first, then classes within schools, and then students within classes, and separate questionnaires were given to the students’ parents, teachers, and school administrators. The TEPS design therefore has at least three levels (students, classes, schools), and because students were followed over several waves, a fourth level (measurement occasions within students) appears as soon as change is studied. Teacher reports add a further complication: a teacher describes a class, not an individual student, so teacher variables are class-level variables even when they are stored on every student’s row. The Taiwan Case boxes of later chapters return to TEPS as a design example; no TEPS data are analyzed in this book.

1.4.2Experiments#

Hierarchies are also created by experiments. When patients are randomly assigned to clinics, or pupils to classes, the random assignment ensures that the groups do not differ systematically at the outset, but it does not make the members of a group independent afterward. Patients treated in the same clinic share the same staff and procedures; pupils in the same class share the same teacher and influence one another. Grouping may be random at the start of a study, and yet over the course of the study the groups become differentiated, because the group affects its members and the members affect the group.

When whole groups (schools, clinics, communities) are randomly assigned to conditions, the design is a cluster-randomized or group-randomized trial (Murray, 1998; Raudenbush, 1997). The treatment then varies only between clusters, and, as Section 1.6 shows, the number of clusters, more than the number of individuals, determines how much information the study has about the treatment effect.

1.4.3Longitudinal, intensive, and crossed designs#

In longitudinal studies the same persons are measured repeatedly, so occasions are nested within persons. Traditional repeated-measures analysis of variance requires a balanced design: the same occasions for all persons and no missing data. A multilevel analysis allows persons to be measured at different times and on different numbers of occasions, can include explanatory variables that change over time, and can model how the residuals of one person are correlated over time (Hedeker & Gibbons, 2006; Singer & Willett, 2003). Diary and experience-sampling studies, in which persons report their mood or behavior several times a day over several weeks, produce dozens or hundreds of occasions per person and are analyzed with the same models (Bolger & Laurenceau, 2013). In experimental psychology, designs in which each participant responds to many stimuli produce trials classified by both participant and stimulus, and treating stimuli as a random factor is necessary whenever the conclusions are meant to generalize to other stimuli (Judd et al., 2012).

1.5A First Look at Multilevel Data#

1.5.1Two schools from NELS:88#

NELS:88 was conducted by the National Center for Education Statistics of the U.S. Department of Education. Its first wave, collected in the spring of 1988, measured about 25,000 eighth-grade students in about 1,000 schools. Kreft & de Leeuw (1998) made available a subsample of 23 schools with 519 students, which this book uses repeatedly (the data are described in Appendix C). The outcome is a mathematics test score, and one pupil-level explanatory variable is the time a pupil reports spending on mathematics homework each week, recorded in eight ordered categories from “none” to “10 or more hours.”

Figure 1.3 plots mathematics score against homework time for two of the 23 schools, and Table 1.4 summarizes them. School 24725 is a public school whose pupils have below-average socioeconomic status (SES); school 62821 is a private school whose pupils have well-above-average SES. The two schools differ in almost every respect. Their mean mathematics scores differ by about 19 points. Pupils in the private school report more homework (a mean category code of 3.3 versus 1.6; the codes are ordered categories, not hours). The spread of scores is larger in the public school. Most importantly for the question that opened the chapter, the relation between homework and mathematics differs: the least-squares slope is about 5.6 points per homework category in the public school and about 1.1 points in the private school.

Figure 1.3: Mathematics Score by Weekly Homework Time in Two NELS:88 Schools — Figure 1.3 Mathematics Score by Weekly Homework Time in Two NELS:88 Schools

Table 1.4 Two NELS:88 Schools: Descriptive Statistics and Within-School Regressions of Mathematics on Homework

		Mathematics				Within-school OLS
School	$n_{j}$	Mean	SD	Homework	SES	Intercept	Slope
24725 (public)	22	43.5	10.0	1.64	$-$ 0.70	34.4	5.59
62821 (private)	67	62.8	5.7	3.30	1.04	59.2	1.09
Pooled (both schools)	89	58.1	10.9	2.89	0.61	48.0	3.49

Note. $n_{j}$ = number of pupils; Homework = mean of the eight-category homework code (0–7); SES = mean standardized socioeconomic status; Intercept and Slope = ordinary least-squares (OLS) regression of mathematics on the homework code.

If the pupils of the two schools are pooled and a single regression line is fitted, the slope is about 3.5. This pooled slope describes neither school. It mixes the relation between homework and mathematics within schools with the difference between schools, which differ in both homework and mathematics for reasons that may have little to do with homework itself (sector, SES, and whatever else distinguishes the two schools). This small example already raises all three of the reviewer’s questions from the start of the chapter.

1.5.2Fifty schools from the Junior School Project#

The Junior School Project (JSP; Mortimore et al., 1988) followed nearly 2,000 pupils in 50 inner-London primary schools through the four junior-school years (1980–1984), and its data have been used to illustrate multilevel models since Goldstein (1987). The widely distributed version used here contains test scores from three of those years: 3,236 records on 1,192 pupils in 49 schools, of whom 1,093 have mathematics scores in both year 0 and year 1.

Figure 1.4 shows two views of the relation between mathematics in year 0 and mathematics in year 1. Panel (a) shows the pupils of one school (school 1). The relation is strong ( $R^{2} = .70$ for the 32 pupils) and approximately linear. Panel (b) shows the least-squares line for each of the 49 schools, drawn over the range of year-0 scores observed in that school. About 80% of the within-school $R^{2}$ values lie between .35 and .82 (the 10th and 90th percentiles), so in most schools prior achievement predicts later achievement well. The lines, however, are not identical: they differ in height (pupils with the same year-0 score tend to score higher in year 1 in some schools than in others) and in slope (the median slope is 0.85, but slopes range from near zero to 1.3). The first kind of difference is the concern of school-effectiveness research: which schools do better than expected given their intake (Aitkin & Longford, 1986)? The second kind of difference asks whether schools are equally effective for pupils with low and high prior achievement.

Figure 1.4: Mathematics in Year 1 Regressed on Mathematics in Year 0 in the Junior School Project — Figure 1.4 Mathematics in Year 1 Regressed on Mathematics in Year 0 in the Junior School Project

With two schools, or even with five, separate regressions for each school are a sensible way to describe the data, and they are the right analysis if the only interest is in those particular schools. With 49 schools, however, the separate regressions, which are noisiest in the smallest schools, do not answer the questions a researcher usually asks: How much do schools differ in general? Which school characteristics explain the differences? How confident can we be about the effect of a particular school, given how few pupils it contributed? Once the schools are regarded as a sample from a larger population of schools, these questions call for a model in which the school-specific intercepts and slopes are themselves random quantities with a distribution. That model is the multilevel model, and the random intercept and random slope versions of it are the subjects of Chapters 5 and 6.

1.6Dependence and Its Consequences#

1.6.1The intraclass correlation#

The simplest way to describe dependence in nested data is to divide the variance of the outcome into two parts. Write the score of pupil $i$ in school $j$ as the sum of an overall mean $γ_{00}$ , a school effect $u_{0 j}$ , and a pupil-level residual $e_{ij}$ :

y_{ij} = γ_{00} + u_{0 j} + e_{ij}, u_{0 j} \sim N (0, τ_{0}^{2}), e_{ij} \sim N (0, σ^{2}), (1.1)

with $u_{0 j}$ and $e_{ij}$ independent of each other. The between-school variance $τ_{0}^{2}$ measures how much the schools’ true means vary around the overall mean, and the within-school variance $σ^{2}$ measures how much pupils differ from the mean of their own school. Equation (1.1) is the empty (or null, or unconditional) multilevel model; it is the one-way random-effects analysis of variance model studied in Chapter 3 and the starting point of Chapter 5.

The intraclass correlation (ICC) is the proportion of the total variance that lies between clusters:

ρ_{I} = \frac{τ _{0}^{2}}{τ _{0}^{2} + σ ^{2}} . (1.2)

The ICC has two readings that are worth keeping apart. It is the proportion of variance in the outcome that is due to differences between clusters, and it is also the correlation between the outcomes of two different level-1 units drawn at random from the same cluster. Both readings follow from Equation (1.1): two pupils in school $j$ share the term $u_{0 j}$ , so their covariance is $τ_{0}^{2}$ ; because each score has variance $τ_{0}^{2} + σ^{2}$ , their correlation is $τ_{0}^{2} / (τ_{0}^{2} + σ^{2})$ .

For the 23 NELS:88 schools, the empty model gives estimates $\overset{τ}{^}_{0}^{2} = 26.1$ and $\overset{σ}{^}^{2} = 81.2$ , so the estimated ICC is $26.1/ (26.1 + 81.2) = .24$ . About a quarter of the variance in mathematics scores lies between schools; equivalently, the scores of two pupils from the same school correlate about .24. For achievement outcomes in national samples of U.S. schools, unconditional ICCs between about .15 and .25 are typical, and smaller values occur in more homogeneous samples of schools (Hedges & Hedberg, 2007). In repeated-measures data the ICC is often much larger, because stable differences between persons account for a large share of the variance.

1.6.2The design effect and the effective sample size#

Dependence reduces the amount of independent information in the data. If pupils within a school were perfect copies of each other ( $ρ_{I} = 1$ ), a school with 25 pupils would contribute no more information about the population mean than a single pupil. If pupils were completely independent ( $ρ_{I} = 0$ ), the school would contribute 25 independent observations. For values in between, the loss of information is summarized by the design effect (Kish, 1965). For clusters of equal size $n$ ,

DEFF = 1 + (n - 1) ρ_{I}, (1.3)

and the effective sample size is $N_{eff} = N / DEFF$ . The design effect is the factor by which the sampling variance of a mean is larger under cluster sampling than under simple random sampling of the same number of individuals (see the Technical Detail box below). For NELS:88, with the average cluster size of 22.6 pupils and the ICC to three decimals, .243, the design effect is about $1 + 21.6 \times .243 \approx 6.25$ (an approximation, because the schools differ in size), so the 519 pupils carry about as much information about the overall mean as $519/6.25 \approx 83$ independently sampled pupils would.

Technical detail Where the design effect comes from

Under Equation (1.1), the mean of the $n$ pupils in school $j$ is $\overset{y}{ˉ}_{j} = γ_{00} + u_{0 j} + \overset{e}{ˉ}_{j}$ , so $Var (\overset{y}{ˉ}_{j}) = τ_{0}^{2} + σ^{2} / n$ . The overall mean of $J$ schools of equal size is the average of $J$ independent school means, so

Var (\overset{y}{ˉ}) = \frac{τ _{0}^{2} + σ ^{2} / n}{J} = \frac{τ _{0}^{2} + σ ^{2}}{J n} \cdot \frac{n τ _{0}^{2} + σ ^{2}}{τ _{0}^{2} + σ ^{2}} = \frac{τ _{0}^{2} + σ ^{2}}{N} [1 + (n - 1) ρ_{I}] .

The first factor, $(τ_{0}^{2} + σ^{2}) / N$ , is the variance of a mean of $N$ independent observations with the same total variance. The second factor is the design effect of Equation (1.3). For a regression coefficient the inflation depends also on how strongly the predictor is clustered: with equal cluster sizes it is approximately $1 + (n - 1) ρ_{x} ρ_{e}$ , where $ρ_{x}$ is the ICC of the predictor and $ρ_{e}$ that of the residual (Moulton, 1990). A predictor measured at the cluster level has $ρ_{x} = 1$ and suffers the full design effect. For a predictor that varies only within clusters, the correlation between the values of two members of a cluster is $- 1/ (n - 1)$ , so the factor is $1 - ρ_{e}$ : its OLS standard error is slightly too large, not too small.

1.6.3What happens to significance tests#

Ordinary least-squares (OLS) regression, the $t$ test, and analysis of variance all assume that observations are independent. When the data are clustered and this assumption is ignored, the formulas treat the $N$ observations as if they were $N$ independent pieces of information. For predictors that are themselves clustered, and above all for cluster-level predictors, the estimated standard errors are then too small, confidence intervals are too narrow, and null hypotheses are rejected too often (Barcikowski, 1981; Kenny & Judd, 1986).

Figure 1.5 and Table 1.5 show the size of the problem in a small simulation (the script, with its seeds and all replication-level results, is in the book’s repository; see the Lab section at the end of the chapter). Forty clusters were generated, half assigned to a “treatment” and half to a “control” condition, as in a cluster-randomized trial. The true treatment effect was zero, so every significant result is a Type I error. The outcome ICC was set to 0, .05, .10, or .20, and the cluster size to 5, 10, 25, or 50. Each condition was replicated 5,000 times, and each data set was analyzed in two ways: by an OLS $t$ test on all individuals, which ignores clustering, and by a $t$ test on the 40 cluster means, which respects it (with equal cluster sizes this cluster-level test, with $J - 2 = 38$ degrees of freedom, is identical to the REML test of the treatment effect with the same degrees of freedom in a random intercept model, Chapter 5, whenever the estimated between-cluster variance is positive).

Figure 1.5: Type I Error Rate of a Test of a Cluster-Level Treatment When Clustering Is Ignored — Figure 1.5 Type I Error Rate of a Test of a Cluster-Level Treatment When Clustering Is Ignored

Table 1.5 Type I Error Rates in the Simulation of Figure 1.5

			Type I error rate ( $α = .05$ )
ICC	$n_{j}$	Design effect	Predicted (OLS)	OLS	Cluster-level
.00	5	1.00	.050	.051	.052
.00	10	1.00	.050	.046	.046
.00	25	1.00	.050	.053	.050
.00	50	1.00	.050	.048	.049
.05	5	1.20	.074	.080	.056
.05	10	1.45	.104	.098	.046
.05	25	2.20	.186	.192	.050
.05	50	3.45	.291	.286	.049
.10	5	1.40	.098	.099	.049
.10	10	1.90	.155	.156	.054
.10	25	3.40	.288	.295	.052
.10	50	5.90	.420	.422	.049
.20	5	1.80	.144	.140	.050
.20	10	2.80	.241	.234	.046
.20	25	5.80	.416	.420	.060
.20	50	10.80	.551	.566	.053

Note. ICC = intraclass correlation of the outcome; $n_{j}$ = cluster size (40 clusters in every condition); Design effect $= 1 + (n_{j} - 1) ICC$ ; Predicted (OLS) = two-sided rejection rate implied by the design effect, $2 Φ (- 1.96/ DEFF)$ ; OLS = test on individuals ignoring clustering; Cluster-level = test on the 40 cluster means. Monte Carlo standard errors are at most .007.

When the ICC is zero, both tests hold the nominal .05 level, as they should. As soon as the ICC is positive, the OLS test rejects too often, and the excess grows with both the ICC and the cluster size. With an ICC of only .05 and 50 individuals per cluster, the OLS test rejects a true null hypothesis in 29% of the samples; with an ICC of .20 and clusters of 50, it rejects in 57%. The rejection rates agree closely with the values predicted from the design effect, which shows that the design effect accounts for nearly all of the distortion here: the OLS standard error is too small by the factor $DEFF$ . The cluster-level test stays between .046 and .060 in all 16 conditions; only one value (.060, for ICC = .20 and $n_{j} = 25$ ) lies more than three Monte Carlo standard errors (about .003 each) from .05. The lesson is general: for effects of cluster-level variables, the number of clusters, more than the number of individuals, limits the effective sample size.

1.6.4A classic example#

The problem is not only theoretical. Bennett (1976) studied the progress in reading, mathematics, and English of primary-school pupils taught by teachers with different teaching styles and reported, on the basis of analyses with pupils as the units of analysis, that pupils taught in a “formal” style made significantly more progress than pupils taught in other styles. Aitkin et al. (1981) reanalyzed the data with a variance-components model in which pupils were nested within teachers. Once the dependence among pupils taught by the same teacher was represented, the differences between teaching styles were no longer statistically significant. The explanation is the one given above. Children within a classroom tend to be similar in their performance, so each child provides less information than an independently sampled child would. The teacher, not the pupil, was the unit at which teaching style varied, and the comparison of styles is a comparison of teachers. Adding pupils per teacher makes each teacher’s mean more precise, but the precision of the comparison between teachers is limited mainly by the number of teachers. The reanalysis, and the discussion that followed it, helped to establish multilevel modeling in educational research (Aitkin & Longford, 1986; Goldstein, 2011).

1.7The Unit-of-Analysis Problem#

Before multilevel models were available, researchers with nested data faced a choice between two imperfect analyses (Burstein, 1980; Kreft & de Leeuw, 1998). In a disaggregated analysis, cluster-level variables are copied onto every individual and the data are analyzed at the individual level, as in the teaching-styles study. The dependence is ignored, and tests of cluster-level variables are too liberal, as Figure 1.5 showed. In an aggregated analysis, individual scores are averaged within clusters and the cluster means are analyzed. The dependence problem disappears, but the within-cluster information is discarded, and, more seriously, the results describe clusters, not individuals.

1.7.1Within, between, and total relations#

The central difficulty is that a relation between two individual-level variables can be decomposed into two relations that need not be the same: a within-cluster relation, which describes how individuals who deviate from their cluster mean on $x$ deviate from their cluster mean on $y$ , and a between-cluster relation, which describes how clusters with different means on $x$ differ in their means on $y$ (Cronbach, 1976). The total relation, estimated by pooling all individuals and ignoring clusters, is a weighted mixture of the two.

Figure 1.6 shows an extreme but instructive case with simulated data. Within each of five groups, individuals with higher $x$ have lower $y$ (the within-group slope is $- 0.89$ ). Across groups, groups with a higher mean on $x$ have a higher mean on $y$ (the between-group slope is $+ 0.89$ ). A single regression that ignores groups (panel a) finds a positive slope of $0.70$ , which is dominated by the between-group relation and misrepresents the relation between $x$ and $y$ among individuals in the same group. Only when the within and between components are separated (panel b) does the structure become visible. Situations in which a relation reverses its sign when a grouping variable is taken into account are instances of Simpson’s paradox (Simpson, 1951).

Figure 1.6: Within-Group, Between-Group, and Total Relations Can Differ in Sign — Figure 1.6 Within-Group, Between-Group, and Total Relations Can Differ in Sign

1.7.2Ecological and atomistic fallacies#

Two errors of inference follow from confusing levels. The ecological fallacy is the error of drawing conclusions about individuals from relations observed among aggregates. Robinson (1950) showed that correlations computed from the averages of U.S. regions or states could differ greatly, even in sign, from the corresponding correlations among individuals. In Figure 1.6, concluding from the positive relation among group means that individuals with higher $x$ tend to have higher $y$ would be an ecological fallacy. The reverse error, the atomistic (or individualistic) fallacy, draws conclusions about groups from relations observed among individuals (Alker, 1969; Hox et al., 2018). Concluding from the negative within-group slope that groups with a higher mean on $x$ must have a lower mean on $y$ would be an atomistic fallacy.

Neither aggregation nor disaggregation solves the problem, because the problem is not the choice of a unit of analysis. Both relations exist, both can be of substantive interest, and both can be estimated. Multilevel models allow within-cluster and between-cluster relations to be estimated simultaneously and to differ, and they make it possible to ask whether they differ. How to do this, and how the answer depends on how the predictor is centered, is the subject of Chapter 7.

1.8Multilevel Theories and Propositions#

A multilevel model is a statistical representation of a substantive theory that involves more than one level. It is therefore worth stating the theory before writing any equation. Snijders & Bosker (2012, Ch. 2) give a compact way to do this, in the spirit of the macro–micro diagrams of Coleman (1990). Propositions are classified by the levels of the variables they connect; in this notation, macro-level (level-2) variables are written in capital letters and micro-level (level-1) variables in lower case. Figure 1.7 redraws the main types with UMG plates. The boundary of the inner (pupil) plate plays the role of the dotted line that separates the micro level from the macro level in Snijders and Bosker’s diagrams: variables inside the pupil plate are micro-level, and variables in the school plate but outside the pupil plate are macro-level.

Figure 1.7: Types of Multilevel Propositions Drawn With Plates — Figure 1.7 Types of Multilevel Propositions Drawn With Plates

A micro-level proposition (panel a) relates two micro-level (pupil) variables, for example a pupil’s homework time and the same pupil’s mathematics score in a study that sampled pupils within schools. The macro level plays no role in the theory, and dependence among pupils in the same school is a nuisance to be accounted for, not a subject of study. A macro-level proposition (panel b) relates two cluster-level variables, for example the average wealth of a school’s catchment area ( $Z$ ) and the school’s mean achievement ( $Y$ ). When $Y$ is not observed directly but is computed from pupil measurements, as a school mean is, the data are multilevel even though the proposition concerns only schools.

The more interesting propositions connect levels. In a macro-to-micro proposition (panel c), a cluster variable affects an individual outcome: pupils in schools with higher average SES have higher mathematics achievement. In a macro-and-micro-to-micro proposition (panel d), both a cluster variable and an individual variable affect the individual outcome: given the time a pupil spends on homework, pupils in schools with higher average SES achieve more; and given the school’s average SES, pupils who spend more time on homework achieve more. In a cross-level interaction (panel e), a cluster variable changes the relation between two individual variables. Snijders & Bosker (2012) give the example of ability grouping ( $Z$ ): the effect of a pupil’s aptitude ( $x$ ) on achievement ( $y$ ) may be small in classes where pupils are grouped by ability and large in classes where they are not. In the UMG this is drawn by making the slope of $y$ on $x$ a school-specific quantity $β_{j}$ that lives in the school plate, and letting $Z_{j}$ affect it. Chapter 6 estimates exactly this structure as a random slope with a cross-level interaction.

Propositions can also run upward. In a micro-to-macro (or emergent) proposition (panel f), individual characteristics affect a group-level outcome, for example pupils’ achievement affecting their teacher’s experience of stress; Kozlowski & Klein (2000) discuss such emergent processes in organizations. Finally, propositions can be chained (panel g): a teacher’s attitude toward learning ( $W$ ) affects pupils’ attitudes ( $x$ ), which affect their achievement ( $y$ ), which in turn affects the teacher’s professional prestige ( $Z$ ). Such chains explain through which micro-level variables two macro-level variables are associated. Models for upward and chained propositions go beyond the basic multilevel regression model; they are taken up in Chapter 29 on multilevel structural equation modeling.

1.9What Multilevel Models Offer, and What They Cost#

Multilevel models have become the standard tool for nested data because they offer several advantages over both disaggregated and aggregated analyses (Hoffman & Walters, 2022; Raudenbush & Bryk, 2002; Snijders & Bosker, 2012).

Correct inference. They represent the dependence among observations, so standard errors, confidence intervals, and significance tests are valid (under the model’s assumptions and with enough clusters), including tests of cluster-level variables.
Efficient estimation. When the model is correctly specified, they combine within-cluster and between-cluster information efficiently (Chapter 7 shows when the two kinds of information must be kept apart), and they borrow strength across clusters, so that estimates for small clusters are stabilized by information from the other clusters (Chapter 12).
Variables at every level. Explanatory variables measured at any level of the hierarchy can enter one model, and their effects can be tested at the appropriate level.
Within and between effects. Within-cluster and between-cluster relations can be estimated separately, which protects against the ecological and atomistic fallacies (Chapter 7).
Heterogeneity as a finding. Relations that vary across clusters (random slopes) and the explanation of that variation by cluster characteristics (cross-level interactions) become part of the model (Chapter 6).
Unbalanced and incomplete data. Clusters of different sizes, and in longitudinal data persons measured at different times and on different numbers of occasions, are handled without discarding data (Chapters 17 and 19).
Rich theories. Many of the propositions in Figure 1.7 can be translated directly into a model, and the model can be drawn and read as the theory it represents.

These advantages come at a cost. The tools of ordinary regression (estimation, tests, measures of fit, diagnostics) must be modified, and the modifications are not always obvious (Part III). Models with several random effects and cross-level interactions can become complex, and the number of possible models grows quickly with the number of variables and levels (Chapter 14). Estimation can fail or give implausible results, especially with few clusters or with variance components near zero (Chapter 9). And the models’ assumptions about the distribution of the random effects and about their independence from the explanatory variables must be examined rather than taken for granted (Chapters 7 and 15). Much of this book is about using multilevel models well in spite of these costs.

1.10The Range of Applications and a Map of This Book#

The applications of multilevel models extend well beyond pupils in schools. The same models are used for sample surveys with multistage designs, for school and teacher effectiveness research, for longitudinal and repeated-measures data, for discrete outcomes such as pass rates, graduation, or symptom presence, for multivariate outcomes in which several variables are measured on each individual (each response becoming a level-1 unit nested within the individual), for nonlinear growth patterns, for cross-classified designs such as participants crossed with stimuli, for measurement models in which item responses are nested within or crossed with persons (De Boeck & Wilson, 2004), for structural equation models with latent variables at several levels, for meta-analysis, in which effect sizes (or participants) are nested within studies, and for models in which the variability of an outcome, not only its mean, differs across persons or groups. Most of these applications have a chapter of their own in this book.

Table 1.6 gives the plan of the book. Part I lays the foundations: this chapter, a chapter on notation and on the three ways of writing (and the graphical way of drawing) a multilevel model, a chapter that derives the multilevel model from the analysis of variance and regression models that readers already know, and a chapter on exploring and managing multilevel data. Part II develops the core models. Part III opens the machinery of estimation and inference. Part IV turns the machinery into a workflow for building, checking, designing, and reporting analyses. Parts V to VII extend the models to longitudinal and intensive data, to generalized outcomes and non-nested structures, and to latent variables and causal questions.

Table 1.6 Plan of the Book

Part	Chapters	Main questions
I. Foundations	1–4	Why multilevel models? How are they written, drawn, and explored?
II. The core models	5–8	Random intercepts, random slopes, cross-level interactions, centering, three levels
III. Estimation and inference	9–13	How are the models estimated and tested? What are predicted random effects? Bayesian estimation
IV. Building, checking, designing, reporting	14–18	Model building, diagnostics, power, missing data, weights, reporting
V. Longitudinal and intensive data	19–22	Growth, covariance structures, within-person dynamics, latent curve models
VI. Generalized and extended models	23–28	Binary, ordinal, count, and survival outcomes; crossed and multiple-membership structures; multivariate and meta-analytic models; variability as an outcome
VII. Latent variables, causality, frontiers	29–32	Multilevel SEM and measurement, causal inference, alternatives and extensions

Note. Parts I–IV together with Chapters 19, 20, and 23 form the core of a one-semester graduate course; the remaining chapters serve a second course or independent study.

Each later chapter follows the same pattern. It opens with a substantive question and a data set, develops the concepts, writes each new model in three ways, draws it as a Unified Model Graph, shows the data displays the model implies, works through an example with annotated R code, and closes with advice on interpretation and reporting, common pitfalls, a software crosswalk for SAS, SPSS, Stata, Mplus, HLM, and JASP, a summary, and exercises. Longer derivations are placed in Technical Detail boxes, which can be skipped on a first reading.

Software crosswalk Estimating the intraclass correlation with the empty model

The empty model of Equation (1.1) is fitted to the NELS:88 data (file nels88_school23.csv, variables math and school). All programs below except Mplus use restricted maximum likelihood (REML) and should reproduce $\overset{τ}{^}_{0}^{2} = 26.1$ , $\overset{σ}{^}^{2} = 81.2$ , and $\overset{ρ}{^}_{I} = .24$ ; Mplus uses maximum likelihood (with the MLR estimator by default for multilevel models), so its estimates differ slightly. Estimation methods are explained in Chapter 9.

(lme4)

library(lme4)
nels <- read.csv("nels88_school23.csv")
m0 <- lmer(math ~ 1 + (1 | school), data = nels, REML = TRUE)
vc <- as.data.frame(VarCorr(m0))       # tau0^2 (school) and sigma^2 (residual)
vc$vcov[1] / sum(vc$vcov)              # intraclass correlation

(PROC MIXED)

SAS

proc mixed data=nels method=reml covtest;
  class school;
  model math = / solution;
  random intercept / subject=school;
run;

(MIXED)

SPSS

MIXED math
  /FIXED = | SSTYPE(3)
  /METHOD = REML
  /PRINT = SOLUTION TESTCOV
  /RANDOM = INTERCEPT | SUBJECT(school) COVTYPE(VC).

Stata

import delimited nels88_school23.csv, clear
mixed math || school: , reml
estat icc

(a copy of the CSV file without the header row; the section “Estimated Intraclass Correlations for the Y Variables” in the summary of the data reports the ICC)

Mplus

DATA:     FILE = nels88_school23.dat;
VARIABLE: NAMES = school student sex race homework schtype ses pared math
                  classstr schsize urban geo minority ratio female white public;
          USEVARIABLES = math;
          CLUSTER = school;
ANALYSIS: TYPE = TWOLEVEL BASIC;

HLM 8. Build a two-level MDM file from a pupil file and a school file linked by school; specify the level-1 model MATH $= β_{0} + r$ with no predictors and the level-2 model $β_{0} = γ_{00} + u_{0}$ . The “Final estimation of variance components” table gives $τ_{00}$ (HLM’s name for $τ_{0}^{2}$ ) and $σ^{2}$ .

JASP (Mixed Models module). Choose Linear Mixed Models, enter math as the dependent variable and school as a random-effects grouping factor, and add no fixed-effect predictors (the intercept is included by default). With the default test method (Satterthwaite) the model is fitted by REML. The variance estimates for the school intercept and the residual appear in the random-effects output; the ICC is computed from them by hand. Menu labels differ slightly between versions.

Chapter summary

Data in psychology and the social sciences are often nested (pupils in schools, occasions in persons) or cross-classified (trials by participants and stimuli). Every variable has a level: the level of the units on which it is measured. Plates show the level of each variable at a glance.
Members of a group resemble each other because of shared context, selection, interaction, and common measurement. The intraclass correlation $ρ_{I} = τ_{0}^{2} / (τ_{0}^{2} + σ^{2})$ measures this resemblance; in the NELS:88 schools it is about .24.
Ignoring dependence makes standard errors too small for predictors that are themselves clustered, especially cluster-level predictors. The design effect $1 + (n - 1) ρ_{I}$ gives the factor by which the variance of a mean is inflated; for cluster-level predictors the number of clusters, more than the number of individuals, governs precision. In the simulation, an ICC of .05 with clusters of 50 raised the Type I error rate of an OLS test from .05 to about .29.
Within-cluster, between-cluster, and total relations can differ, even in sign. Inferring individual relations from aggregates is the ecological fallacy; inferring group relations from individual data is the atomistic fallacy.
Research questions can be stated as micro, macro, macro-to-micro, cross-level, or micro-to-macro propositions. The type of proposition guides the choice of model.
Multilevel models give valid inference, efficient estimates, variables at every level, separate within and between effects, and a way to study heterogeneity across clusters, at the cost of more complex specification, estimation, and checking.

Key Terms#

aggregation and disaggregation
atomistic fallacy
between-cluster relation
cluster; cluster-randomized trial
cross-classified structure
cross-level interaction
design effect
ecological fallacy
effective sample size
empty (null) model
intraclass correlation
level; level-1 and level-2 units
macro-level and micro-level propositions
micro-to-macro (emergent) proposition
multiple-membership structure
multistage sampling
plate
Simpson’s paradox
total relation
unit-of-analysis problem
within-cluster relation

Exercises#

A. Concept checks#

For each study, name the levels, the units at each level, and whether the structure is nested, cross-classified, or multiple-membership: (a) 60 therapists each treat 8 clients, and each client completes a symptom checklist weekly for 12 weeks; (b) 150 participants each judge the attractiveness of the same 40 faces; (c) employees rate their team climate, and some employees belong to two teams; (d) 2,000 adolescents from 90 neighborhoods attend 30 high schools.
In a study of pupils in schools, state the level of each variable: mathematics score, school size, pupil’s gender, proportion of girls in the school, teacher’s years of experience (one teacher per class, several classes per school), and the average homework time of the pupils in a class.
An ICC of .08 is reported for an outcome in a sample of 60 schools with 120 pupils each. A colleague concludes that “the clustering is negligible, so ordinary regression is fine.” Explain in two or three sentences why this conclusion is unsafe.
Explain the difference between the ecological fallacy and the atomistic fallacy, and give an example of each from your own field.
Classify each proposition by the types in Figure 1.7: (a) “Classrooms with a more supportive climate have lower average anxiety”; (b) “Pupils’ anxiety is lower in classrooms with a more supportive climate”; (c) “The effect of pupils’ prior achievement on anxiety is weaker in classrooms with a more supportive climate”; (d) “Pupils’ disruptive behavior increases teachers’ burnout.”

B. Analytical exercises#

Using Equation (1.1), show that $Cov (y_{ij}, y_{i^{'} j}) = τ_{0}^{2}$ for two different pupils $i \neq = i^{'}$ in the same school, that $Cov (y_{ij}, y_{i^{'} j^{'}}) = 0$ for pupils in different schools, and hence that $Corr (y_{ij}, y_{i^{'} j}) = ρ_{I}$ . Solution →
Derive the design effect of Equation (1.3) by following the steps in the Technical Detail box. Then show that the effective sample size $N / DEFF$ lies between $J$ and $N$ , and state when each bound is reached. Solution →
A cluster-randomized trial will enroll 30 classes of 25 pupils. Pilot data suggest $ρ_{I} = .15$ . Compute the design effect and the effective sample size. By how much must the number of classes be increased to obtain the same effective sample size as a trial of 30 classes with no clustering? Solution →
Using the approximation $1 + (n - 1) ρ_{x} ρ_{e}$ in the Technical Detail box, explain why ignoring clustering hurts tests of cluster-level predictors much more than tests of predictors that vary mainly within clusters. Solution →

C. Computational exercises#

Using nels88_school23.csv, compute the number of pupils, the mean mathematics score, and the within-school OLS intercept and slope of mathematics on homework for each of the 23 schools. Plot the 23 regression lines in one panel, as in Figure 1.4(b). Describe how the schools differ.
Reproduce Table 1.4 for two other schools of your choice, and comment on whether the pooled slope describes either school.
Fit the empty model to the NELS:88 data in R and in one other program from the crosswalk. Confirm that the estimated ICC is about .24, and compute the design effect using the average cluster size.
Using jsp_wide.csv, compute for each school the correlation between mathematics in year 0 and year 1. Plot the correlations against the number of pupils in the school. What does the plot suggest about the reliability of school-specific estimates from small schools? Solution →

D. Simulation exercise#

Open sim_ch01_typeI_V01.R. (a) Change the number of clusters from 40 to 10 and rerun the simulation. Does the cluster-level test still hold the .05 level? Does the OLS test behave better or worse than with 40 clusters? (b) Add a pupil-level predictor $x_{ij}$ that varies only within clusters (for example, $x_{ij} = z_{ij} - \overset{z}{ˉ}_{j}$ for standard normal $z_{ij}$ ) and whose true effect is zero. Record how often the OLS test of this predictor rejects. Explain the result with the approximation $1 + (n - 1) ρ_{x} ρ_{e}$ . Save the new scripts with the suffix _V02.

E. Applied mini-project#

Choose a published empirical article in your area whose data are nested. (a) Draw the data structure as a unit diagram and with plates, placing every variable used in the analysis in the correct plate. (b) State the main research question as one or more of the proposition types in Figure 1.7. (c) Report the number of clusters and the cluster sizes, and state whether the analysis respected the design. (d) Write a short reviewer comment (about 150 words) on the unit of analysis.

Lab#

References

Aitkin, M., Anderson, D., & Hinde, J. (1981). Statistical modelling of data on teaching styles. Journal of the Royal Statistical Society: Series A (General), 144(4), 419–461. https://doi.org/10.2307/2981826

Aitkin, M., & Longford, N. (1986). Statistical modelling issues in school effectiveness studies. Journal of the Royal Statistical Society: Series A (General), 149(1), 1–43. https://doi.org/10.2307/2981882

Alker, H. R. (1969). A typology of ecological fallacies. In M. Dogan & S. Rokkan (Eds.), Quantitative ecological analysis in the social sciences (pp. 69–86). MIT Press.

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. https://doi.org/10.1016/j.jml.2007.12.005

Barcikowski, R. S. (1981). Statistical power with group mean as the unit of analysis. Journal of Educational Statistics, 6(3), 267–285. https://doi.org/10.3102/10769986006003267

Bennett, N. (1976). Teaching styles and pupil progress. Open Books.

Bolger, N., & Laurenceau, J.-P. (2013). Intensive longitudinal methods: An introduction to diary and experience sampling research. Guilford Press.

Browne, W. J., Goldstein, H., & Rasbash, J. (2001). Multiple membership multiple classification (MMMC) models. Statistical Modelling, 1(2), 103–124. https://doi.org/10.1177/1471082X0100100202

Buntine, W. L. (1994). Operations for learning with graphical models. Journal of Artificial Intelligence Research, 2, 159–225. https://doi.org/10.1613/jair.62

Burstein, L. (1980). The analysis of multilevel data in educational research and evaluation. Review of Research in Education, 8, 158–233. https://doi.org/10.3102/0091732X008001158

Coleman, J. S. (1990). Foundations of social theory. Belknap Press of Harvard University Press.

Cronbach, L. J. (1976). Research on classrooms and schools: Formulation of questions, design, and analysis [Occasional paper]. Stanford Evaluation Consortium, Stanford University.

De Boeck, P., & Wilson, M. (Eds.). (2004). Explanatory item response models: A generalized linear and nonlinear approach. Springer. https://doi.org/10.1007/978-1-4757-3990-9

Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press. https://doi.org/10.1017/CBO9780511790942

Goldstein, H. (1987). Multilevel models in educational and social research. Griffin.

Goldstein, H. (2011). Multilevel statistical models (4th ed.). Wiley. https://doi.org/10.1002/9780470973394

Hedeker, D., & Gibbons, R. D. (2006). Longitudinal data analysis. Wiley. https://doi.org/10.1002/0470036486

Hedges, L. V., & Hedberg, E. C. (2007). Intraclass correlation values for planning group-randomized trials in education. Educational Evaluation and Policy Analysis, 29(1), 60–87. https://doi.org/10.3102/0162373707299706

Hoffman, L. (2015). Longitudinal analysis: Modeling within-person fluctuation and change. Routledge. https://doi.org/10.4324/9781315744094

Hoffman, L., & Walters, R. W. (2022). Catching up on multilevel modeling. Annual Review of Psychology, 73, 659–689. https://doi.org/10.1146/annurev-psych-020821-103525

Hox, J. J., Moerbeek, M., & van de Schoot, R. (2018). Multilevel analysis: Techniques and applications (3rd ed.). Routledge. https://doi.org/10.4324/9781315650982

Judd, C. M., Westfall, J., & Kenny, D. A. (2012). Treating stimuli as a random factor in social psychology: A new and comprehensive solution to a pervasive but largely ignored problem. Journal of Personality and Social Psychology, 103(1), 54–69. https://doi.org/10.1037/a0028347

Kenny, D. A., & Judd, C. M. (1986). Consequences of violating the independence assumption in analysis of variance. Psychological Bulletin, 99(3), 422–431. https://doi.org/10.1037/0033-2909.99.3.422

Kish, L. (1965). Survey sampling. Wiley.

Kozlowski, S. W. J., & Klein, K. J. (2000). A multilevel approach to theory and research in organizations: Contextual, temporal, and emergent processes. In K. J. Klein & S. W. J. Kozlowski (Eds.), Multilevel theory, research, and methods in organizations: Foundations, extensions, and new directions (pp. 3–90). Jossey-Bass.

Kreft, I. G. G., & de Leeuw, J. (1998). Introducing multilevel modeling. Sage. https://doi.org/10.4135/9781849209366

Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963–974. https://doi.org/10.2307/2529876

Longford, N. T. (1993). Random coefficient models. Clarendon Press.

Mortimore, P., Sammons, P., Stoll, L., Lewis, D., & Ecob, R. (1988). School matters: The junior years. Open Books.

Moulton, B. R. (1990). An illustration of a pitfall in estimating the effects of aggregate variables on micro units. The Review of Economics and Statistics, 72(2), 334–338. https://doi.org/10.2307/2109724

Murray, D. M. (1998). Design and analysis of group-randomized trials. Oxford University Press.

Rasbash, J., & Goldstein, H. (1994). Efficient analysis of mixed hierarchical and cross-classified random structures using a multilevel model. Journal of Educational and Behavioral Statistics, 19(4), 337–350. https://doi.org/10.3102/10769986019004337

Raudenbush, S. W. (1997). Statistical analysis and optimal design for cluster randomized trials. Psychological Methods, 2(2), 173–185. https://doi.org/10.1037/1082-989X.2.2.173

Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical linear models: Applications and data analysis methods (2nd ed.). Sage.

Robinson, W. S. (1950). Ecological correlations and the behavior of individuals. American Sociological Review, 15(3), 351–357. https://doi.org/10.2307/2087176

Searle, S. R., Casella, G., & McCulloch, C. E. (1992). Variance components. Wiley. https://doi.org/10.1002/9780470316856

Simpson, E. H. (1951). The interpretation of interaction in contingency tables. Journal of the Royal Statistical Society: Series B (Methodological), 13(2), 238–241. https://doi.org/10.1111/j.2517-6161.1951.tb00088.x

Singer, J. D., & Willett, J. B. (2003). Applied longitudinal data analysis: Modeling change and event occurrence. Oxford University Press. https://doi.org/10.1093/acprof:oso/9780195152968.001.0001

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

Verbeke, G., & Molenberghs, G. (2000). Linear mixed models for longitudinal data. Springer. https://doi.org/10.1007/978-1-4419-0300-6

Yu, H.-T. (2026a). The Unified Model Graph: Drawing the statistical models of psychology [Unpublished manuscript].

Yu, H.-T. (2026b). umg: Unified Model Graphs for statistical models in psychology (Version 0.6.1) [Computer software]. https://github.com/hsiutingyu/umg