Germán Rodríguez

Multilevel Models
Princeton University
We now consider 3-level models. We will use time-series data on 1721 students nested within 60 urban public primary schools. The outcome of interest is mathematics achievement. The data were collected at the end of first grade and annually thereafter up to sixth grade, but not all students have six observations.

This dataset is used in Chapter 4 of the manual for HLM 6. It came in three files, one for each level, but I have merged them into a single file available in the course website.

A correspondent asked about the equations for the models fitted here. I have included links to the equations for the growth curve model, and for the growth curve model with predictors, skipping the simpler variance-components model.

We read the data and fit a simple variance-components model. Note that we have two specifications for the random part, one at the school level (level 3), and one at the child-level (level 2). They are both null (intercept-only) models, so we are simply adding an error term at each level. We specify these going *down* the hierarchy. We also specify maximum likelihood estimation. (The R default is REML).

. use https://data.princeton.edu/pop510/egm, clear (Data from HLM 6 manual on math achievement in 60 urban schools) . mixed math || schoolid: || childid: , mle Performing EM optimization: Performing gradient-based optimization: Iteration 0: log likelihood = -12652.991 Iteration 1: log likelihood = -12652.99 Computing standard errors: Mixed-effects ML regression Number of obs = 7,230 ------------------------------------------------------------- ¦ No. of Observations per Group Group Variable ¦ Groups Minimum Average Maximum ----------------+-------------------------------------------- schoolid ¦ 60 18 120.5 387 childid ¦ 1,721 2 4.2 6 ------------------------------------------------------------- Wald chi2(0) = . Log likelihood = -12652.99 Prob > chi2 = . ------------------------------------------------------------------------------ math ¦ Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _cons ¦ -.510021 .0779639 -6.54 0.000 -.6628275 -.3572144 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters ¦ Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ schoolid: Identity ¦ var(_cons) ¦ .317646 .0687236 .207866 .4854037 -----------------------------+------------------------------------------------ childid: Identity ¦ var(_cons) ¦ .5704854 .0337543 .5080202 .6406312 -----------------------------+------------------------------------------------ var(Residual) ¦ 1.523894 .0290756 1.467959 1.58196 ------------------------------------------------------------------------------ LR test vs. linear model: chi2(2) = 1404.55 Prob > chi2 = 0.0000 Note: LR test is conservative and provided only for reference.

> library(foreign) > egm <- read.dta("https://data.princeton.edu/pop510/egm.dta") > library(lme4) > vc <- lmer(math ~ 1 + (1|childid) + (1 | schoolid), data=egm, REML=FALSE) > summary(vc, corr = FALSE) Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: math ~ 1 + (1 | childid) + (1 | schoolid) Data: egm AIC BIC logLik deviance df.resid 25314.0 25341.5 -12653.0 25306.0 7226 Scaled residuals: Min 1Q Median 3Q Max -2.6726 -0.6576 0.0265 0.6735 3.3344 Random effects: Groups Name Variance Std.Dev. childid (Intercept) 0.5705 0.7553 schoolid (Intercept) 0.3176 0.5636 Residual 1.5239 1.2345 Number of obs: 7230, groups: childid, 1721; schoolid, 60 Fixed effects: Estimate Std. Error t value (Intercept) -0.51002 0.07796 -6.542

We have 60 schools with an average of 120.5 students for a total of 1721 children, who have between 2 and 6 observations each for a total of 7230 observations.

Our estimate of the mean score is -0.51 and the standard deviations are 0.56, 0.76 and 1.23 at the school, child and observation levels. We can turns these into intraclass correlations from first principles.

. mat list e(b) // to get names e(b)[1,4] math: lns1_1_1: lns2_1_1: lnsig_e: _cons _cons _cons _cons y1 -.51002095 -.57340893 -.28063386 .21063435 . scalar v3 = exp(2 * _b[lns1_1_1:_cons]) . scalar v2 = exp(2 * _b[lns2_1_1:_cons]) . scalar v1 = exp(2 * _b[lnsig_e:_cons]) . display v3/(v3 + v2 + v1), (v2 + v3)/(v3 + v2 + v1) .13169264 .36820984

> v <- c(sigma(vc)^2, unlist(VarCorr(vc))) > v[3]/(v[1] + v[2] + v[3]) schoolid 0.1316919 > (v[2] + v[3])/(v[1] + v[2] + v[3]) childid 0.3682094

We find intra-class correlations of 0.13 at the school level and 0.37 at the child level. These numbers represent the correlation in math achievement between two children in the same school, and between two measurements on the same child (in the same school). We can also conclude that 13% of the variation in math achievement can be attributed to the schools and 37% to the children (which includes the school).

We can also obtain posterior Bayes estimates of the school and child random effects (these come up in the same order as the equations) and compare with school and child means

. predict r3 r2, reffects . // 3 . egen m3 = mean(math), by(schoolid) . gen eb3 = _b[_cons] + r3 . egen lev3 = tag(schoolid) . twoway (scatter eb3 m3 if lev3) (function y=x, range(-2.5 .8)), /// > legend(off) title(School Means) name(school, replace) . // 2 . egen m2 = mean(math), by (childid) . egen n2 = count(math), by(childid) //? . gen eb2 = eb3 + r2 . egen lev2 = tag(childid) . twoway (scatter eb2 m2 if lev2 & n2 > 2) (function y=x, range(-3.6 3.1)), /// > legend(off) title(Child Means) name(child, replace) . // graph . graph combine school child, xsize(6) ysize(3) . graph export egmfig1.png, width(700) replace (file egmfig1.png written in PNG format)

> library(dplyr) > library(ggplot2) > library(gridExtra) > re <- ranef(vc) > # schools > schools <- group_by(egm, schoolid) %>% + summarize(ml3 = mean(math)) %>% + mutate(eb3 = fixef(vc) + re$schoolid[schoolid, ] ) > r <- data.frame(x=range(schools$ml3)) > g3 <- ggplot(schools, aes(ml3, eb3)) + geom_point() + + geom_line(data=r, aes(x,x)) > # students > students <- group_by(egm, childid) %>% + summarize(ml2 = mean(math), schoolid = first(schoolid) ) %>% + mutate(eb2 = fixef(vc) + re$schoolid[schoolid,1] + re$childid[childid,1] ) > r <- data.frame(x=range(students$ml2)) > g2 <- ggplot(students, aes(ml2, eb2)) + geom_point() + + geom_line(data=r, aes(x,x)) > g <- arrangeGrob(g3, g2, ncol=2) > ggsave("egmfig1r.png", plot = g, width=10, height=5, dpi=72)

We see the usual shrinkage, particularly for the child means.

The obvious predictor at the observation level is study `year`

, which is coded -2.5, -1.5, -0.5, 0.5, 1.5 and 2.5 for the six years of data collection. There’s also a variable called `grade`

reflecting the grade minus one, which is almost the same thing:

. tab year grade ¦ grade year ¦ 0 1 2 3 4 5 ¦ Total -----------+------------------------------------------------------------------+---------- -2.5 ¦ 131 0 0 0 0 0 ¦ 131 -1.5 ¦ 1,346 0 0 0 0 0 ¦ 1,346 -.5 ¦ 104 1,412 4 0 0 0 ¦ 1,520 .5 ¦ 1 194 1,470 7 0 0 ¦ 1,672 1.5 ¦ 0 5 174 1,200 8 0 ¦ 1,387 2.5 ¦ 0 1 5 149 1,010 9 ¦ 1,174 -----------+------------------------------------------------------------------+---------- Total ¦ 1,582 1,612 1,653 1,356 1,018 9 ¦ 7,230

> table(egm$year, egm$grade) 0 1 2 3 4 5 -2.5 131 0 0 0 0 0 -1.5 1346 0 0 0 0 0 -0.5 104 1412 4 0 0 0 0.5 1 194 1470 7 0 0 1.5 0 5 174 1200 8 0 2.5 0 1 5 149 1010 9

We will follow the HLM manual in using years to track progress in math achievement over time. Year zero is near the start of the third grade, and I will refer to it as ‘third grade’ for simplicity.

Our first model will allow for random intercepts and slopes, with no predictors at levels 2 or 3. I specify `covariance(unstructured)`

at each level, so the intercept and slope can be correlated at level 2 as well as level 3.

The equations for this model are here.

. mixed math year /// > || schoolid: year, covariance(unstructured) /// > || childid: year, covariance(unstructured) mle Performing EM optimization: Performing gradient-based optimization: Iteration 0: log likelihood = -8169.1304 Iteration 1: log likelihood = -8163.1176 Iteration 2: log likelihood = -8163.1156 Iteration 3: log likelihood = -8163.1156 Computing standard errors: Mixed-effects ML regression Number of obs = 7,230 ------------------------------------------------------------- ¦ No. of Observations per Group Group Variable ¦ Groups Minimum Average Maximum ----------------+-------------------------------------------- schoolid ¦ 60 18 120.5 387 childid ¦ 1,721 2 4.2 6 ------------------------------------------------------------- Wald chi2(1) = 2499.87 Log likelihood = -8163.1156 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ math ¦ Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- year ¦ .7630273 .0152609 50.00 0.000 .7331164 .7929382 _cons ¦ -.7793053 .0578295 -13.48 0.000 -.8926489 -.6659616 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters ¦ Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ schoolid: Unstructured ¦ var(year) ¦ .0110171 .0025621 .0069842 .0173788 var(_cons) ¦ .1653165 .0359122 .1079959 .2530611 cov(year,_cons) ¦ .0170462 .0071472 .003038 .0310545 -----------------------------+------------------------------------------------ childid: Unstructured ¦ var(year) ¦ .0112561 .001961 .0080001 .0158373 var(_cons) ¦ .6404597 .025183 .5929559 .6917691 cov(year,_cons) ¦ .0467854 .005063 .0368621 .0567086 -----------------------------+------------------------------------------------ var(Residual) ¦ .3014384 .0066403 .2887005 .3147382 ------------------------------------------------------------------------------ LR test vs. linear model: chi2(6) = 5595.58 Prob > chi2 = 0.0000 Note: LR test is conservative and provided only for reference. . estat recov, corr Random-effects correlation matrix for level schoolid ¦ year _cons -------------+---------------------- year ¦ 1 _cons ¦ .3994262 1 Random-effects correlation matrix for level childid ¦ year _cons -------------+---------------------- year ¦ 1 _cons ¦ .5510235 1

> gc <- lmer(math ~ year + (1 + year|childid) + (1 + year|schoolid), + data = egm, REML=FALSE) > summary(gc) Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: math ~ year + (1 + year | childid) + (1 + year | schoolid) Data: egm AIC BIC logLik deviance df.resid 16344.2 16406.2 -8163.1 16326.2 7221 Scaled residuals: Min 1Q Median 3Q Max -3.2047 -0.5640 -0.0372 0.5309 5.2629 Random effects: Groups Name Variance Std.Dev. Corr childid (Intercept) 0.64046 0.8003 year 0.01126 0.1061 0.55 schoolid (Intercept) 0.16532 0.4066 year 0.01102 0.1050 0.40 Residual 0.30144 0.5490 Number of obs: 7230, groups: childid, 1721; schoolid, 60 Fixed effects: Estimate Std. Error t value (Intercept) -0.77931 0.05783 -13.48 year 0.76303 0.01526 50.00 Correlation of Fixed Effects: (Intr) year 0.357

The first thing you will notice is that as models become more complex they take a bit longer to run. The results can be compared with the output from HLM 6 starting on page 80 of the manual, and are in very close agreement.

The fixed-effect estimates yield an intercept of -0.78 and a slope of 0.76. These values can be interpreted as the average achievement in third grade, and the average rate of change per year. Clearly average math achievement increases significantly over time.

The output also shows the variances and covariances of the random effects. The correlations can be obtained with `estat recovariance`

, which has an option `corr`

.

We see that the intercepts and slopes have correlations of 0.40 at the school level and 0.55 at the child level. (In a hierarchical model the correlation at level 2 is never lower than at level 3.) The child estimate (0.55) can be interpreted as the correlation between a child’s achievement around third grade and her rate of growth. The school estimate (0.40) represents the correlation between a school’s achievement around third grade and the rate of change per year.

Most of the variance in the intercept comes from the child level, whereas the variance in the slope is about the same at each level. These results indicate fairly significant variation among schools in mean achievement around third grade, as well as in the rates of change over time.

We can now estimate growth curves for each school and for each child, combining the estimated fixed and random effects. I’ll just predict the random effects, and plot the slope versus the intercept at each level. Note that at the child level I include both the school and child residuals. In Stata we use `predict`

with the `reffects`

option. With four random effects we need four variables to store the Bayes estimates. I’ll use numbers for the levels, `a`

for intercepts and b` for slopes.

. predict r3b r3a r2b r2a, reffects . local labels xtitle(intercept) ytitle(slope) . gen a3 = _b[_cons] + r3a . gen b3 = _b[year] + r3b . twoway scatter b3 a3 if lev3, `labels' title(School) /// > name(school, replace) . gen a2 = a3 + r2a // or _b[_cons] + r3a + r2a . gen b2 = b3 + r2b // or _b[year] + r3b + r3b . twoway (scatter b2 a2 if lev2, `labels') /// > (scatter b2 a2 if lev2 & b3 < .55), legend(off) /// > title(Child) name(child, replace) . graph combine school child, xsize(6) ysize(3) /// > title(School and Child Random Intercepts and Slopes) . graph export egmfig2.png, width(700) replace (file egmfig2.png written in PNG format)

> b = fixef(gc) > gcre <- ranef(gc) > names(gcre$schoolid) <- names(gcre$childid) <- c("intercept","slope") > school <- mutate(gcre$schoolid, intercept = b[1] + intercept, slope = b[2] + slope) > g3 <- ggplot(school, aes(intercept,slope)) + geom_point() > child <- select(egm, schoolid, childid) %>% + group_by(childid) %>% filter(row_number()==1) %>% mutate( + intercept = b[1] + gcre$schoolid[schoolid,1] + gcre$childid[childid,1], + slope = b[2] + gcre$schoolid[schoolid,2] + gcre$childid[childid,2]) > g2 <- ggplot(child, aes(intercept,slope)) + geom_point() + + geom_point(data=filter(child, slope < .55), aes(intercept,slope), color="blue") > g <- arrangeGrob(g3, g2, ncol=2) > ggsave("egmfig2r.png", plot = g, width=10, height=5, dpi=72)

We see substantial variation in school levels and rates of change in math achievement. One school is notable for having a reasonable third-grade mean, but a low rate of change per year (just above 0.5). I highlight the children from this school in the graph on the right, which shows their expected third-grade scores and rates of change combining school and child effects.

Note the positive correlation between intercepts and slopes at both levels, so children and schools with higher third-grade math achievement also tend to show greater gains over time.

We can also plot fitted lines at the school and child level. It might not be a bad idea to test for linearity of the growth curves.

One way to look for school or individual characteristics that may be associated with levels and trends in mathematics achievement is to plot or regress the Bayes estimates of intercepts and slopes on potential predictors at each level.

Here I’ll follow the HLM manual in considering a model where a child’s growth curve depends on ethnicity, with different intercepts and slopes for whites, blacks and hispanics, and where the school average curve depends on the percent of students with low income.

One way to develop a model of this type is to write down the equation for each level and then substitute down the hierarchy to obtain the fixed and random components in the complete model. This process shows the need to create cross-level interactions between `year`

and the predictors at the school and child level.

Before creating cross-products we consider centering. This is not needed (not, in my view, advisable) for dummy variables, where zero is meaningful because it identifies the reference category. One might consider centering the percent poor, but for consistency with the HLM manual I will not. Bear in mind that the main effects then refer to

The equations for this model are here.

. gen yearBlack = year * black . gen yearHispanic = year * hispanic . gen yearLowinc = year * lowinc . mixed math year black hispanic lowinc /// > yearBlack yearHispanic yearLowinc /// > || schoolid: year, covariance(unstructured) /// > || childid: year, covariance(unstructured) mle Performing EM optimization: Performing gradient-based optimization: Iteration 0: log likelihood = -8125.8349 Iteration 1: log likelihood = -8119.6062 Iteration 2: log likelihood = -8119.6035 Iteration 3: log likelihood = -8119.6035 Computing standard errors: Mixed-effects ML regression Number of obs = 7,230 ------------------------------------------------------------- ¦ No. of Observations per Group Group Variable ¦ Groups Minimum Average Maximum ----------------+-------------------------------------------- schoolid ¦ 60 18 120.5 387 childid ¦ 1,721 2 4.2 6 ------------------------------------------------------------- Wald chi2(7) = 3324.78 Log likelihood = -8119.6035 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ math ¦ Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- year ¦ .8745122 .0391404 22.34 0.000 .7977984 .9512259 black ¦ -.5021085 .0778753 -6.45 0.000 -.6547414 -.3494756 hispanic ¦ -.3193816 .0860935 -3.71 0.000 -.4881218 -.1506414 lowinc ¦ -.0075778 .0016908 -4.48 0.000 -.0108918 -.0042638 yearBlack ¦ -.0309253 .0224586 -1.38 0.169 -.0749432 .0130927 yearHispanic ¦ .0430865 .024659 1.75 0.081 -.0052442 .0914172 yearLowinc ¦ -.0013689 .0005226 -2.62 0.009 -.0023933 -.0003446 _cons ¦ .1406379 .1274908 1.10 0.270 -.1092396 .3905153 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters ¦ Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ schoolid: Unstructured ¦ var(year) ¦ .0079801 .0020562 .004816 .0132231 var(_cons) ¦ .0780901 .019642 .0476973 .1278494 cov(year,_cons) ¦ .0008172 .0044663 -.0079366 .009571 -----------------------------+------------------------------------------------ childid: Unstructured ¦ var(year) ¦ .0110938 .0019517 .0078583 .0156615 var(_cons) ¦ .6222512 .0245399 .5759658 .6722562 cov(year,_cons) ¦ .0466258 .0049886 .0368483 .0564033 -----------------------------+------------------------------------------------ var(Residual) ¦ .3015912 .0066415 .2888511 .3148932 ------------------------------------------------------------------------------ LR test vs. linear model: chi2(6) = 4797.28 Prob > chi2 = 0.0000 Note: LR test is conservative and provided only for reference.

> egm <- mutate(egm, + eth = factor(1 + black + 2 * hispanic, labels=c("white","black","hispanic"))) > gcp <- lmer(math ~ year * (eth + lowinc) + + (year | childid) + (year | schoolid), data=egm, REML=FALSE) > summary(gcp, corr=FALSE) Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: math ~ year * (eth + lowinc) + (year | childid) + (year | schoolid) Data: egm AIC BIC logLik deviance df.resid 16269.2 16372.5 -8119.6 16239.2 7215 Scaled residuals: Min 1Q Median 3Q Max -3.2482 -0.5645 -0.0328 0.5348 5.2765 Random effects: Groups Name Variance Std.Dev. Corr childid (Intercept) 0.62225 0.78883 year 0.01109 0.10533 0.56 schoolid (Intercept) 0.07809 0.27945 year 0.00798 0.08933 0.03 Residual 0.30159 0.54917 Number of obs: 7230, groups: childid, 1721; schoolid, 60 Fixed effects: Estimate Std. Error t value (Intercept) 0.1406379 0.1274907 1.103 year 0.8745123 0.0391403 22.343 ethblack -0.5021084 0.0778753 -6.448 ethhispanic -0.3193816 0.0860935 -3.710 lowinc -0.0075778 0.0016908 -4.482 year:ethblack -0.0309253 0.0224586 -1.377 year:ethhispanic 0.0430865 0.0246590 1.747 year:lowinc -0.0013689 0.0005226 -2.619

Results can be compared with the HLM output starting on page 88 of the manual. In this model each child has her own growth curve, but the intercept and slope depend on the school’s percent poor, school unobserved characteristics, the child’s owh ethnicity, and unobserved characteristics of the child.

One way to think about the estimates is as follows. We start with a growth curve that has 0.141 points of math achievement by the middle of third grade and increases by 0.875 points per year, provided the school has no poor students. Each percentage point of low income students is associated with a third grade mean that’s 0.008 lower and increases 0.001 points less per year, so a school with 100% poor students would have a third-grade mean of -0.617 and an increase of 0.738 points per year. The intercept and slope in any particular school will vary around the average line for its percent poor.

Next come characteristics of the individual, in this case ethnicity. Blacks have a third-grade mean 0.502 points below whites in the same school, and the increase per year is 0.031 points less than for whites in the same school. The results for hispanics indicate a third-grade mean 0.319 points below whites, but a growth rate per year 0.043 points higher than whites in the same school. Each child’s growth curve, however, will vary around the average line for kids of the same ethnicity in the same school.

The random effects indicate that the average third-grade score varies by school, with a standard deviation of 0.28, while the rate of change has a standard deviation of 0.089. Interestingly, the correlation between intercept and slope has practically vanished, suggesting that it can be attributed to the school’s socio-economic composition.

Variation at the child level is even more sustantial; expected third-grade scores vary with a standard deviation of 0.789, and the rates of change vary with a standard deviation of 0.105. These two random effects remain highly correlated, so kids who have high levels of math proficiency in third grade also tend to improve faster over time.

It’s useful to compare fixed and random effects on a similar scale. Here I will consider four types of schools, with 50 and 100 poor and random effects one standard deviation above and below the mean. In each school I will consider kids from the three ethnic groups whose intercepts and slopes fall one standard deviation above and below the mean of their ethnic group in their school.

To do this we need to create a prediction dataset with the combinations of interest. I will store these in the first few rows of the dataset.

. local k = 1 . foreach lowinc in 50 100 { 2. foreach eth in "white" "black" "hispanic" { 3. foreach year in -2.5 2.5 { 4. quietly replace lowinc = `lowinc' in `k' 5. quietly replace black = "`eth'" == "black" in `k' 6. quietly replace hispanic = "`eth'" == "hispanic" in `k' 7. quietly replace year = `year' in `k' 8. local k = `k' + 1 9. } 10. } 11. } . keep in 1/12 (7,218 observations deleted) . quietly replace yearBlack = year * black . quietly replace yearHispanic = year * hispanic . quietly replace yearLowinc = year * lowinc

> pdf <- data.frame( + lowinc = rep(c(50,100),c(24,24)), + schoolz = rep(rep(c(-1,1), c(12,12)),2), + black = rep(rep(c(0,1,0),c(2,2,2)),8), + hispanic = rep(rep(c(0,0,1),c(2,2,2)),8), + childz = rep(rep(c(-1,1),c(6,6)), 4), + year = rep(c(-2.5,2.5), 24) + ) %>% mutate( + eth = factor(1 + black + 2 * hispanic, labels=c("white","black","hispanic")) + )

We can then predict the fixed effects and add the random effects. In Stata we first expand the dataset.

. predict xb, xb . expand 4 (36 observations created) . bysort lowinc black hispanic year: gen rep=_n . gen schoolz = -1 + 2 * (rep > 2) . gen childz = -1 + 2 * (rep == 2 | rep ==4) . scalar sa3 = exp(_b[lns1_1_2:_cons]) // school cons . scalar sa2 = exp(_b[lns2_1_2:_cons]) // child cons . scalar sb3 = exp(_b[lns1_1_1:_cons]) // school year . scalar sb2 = exp(_b[lns2_1_1:_cons]) // child year . gen fv = xb + (schoolz * sa3 + childz * sa2) + /// > (schoolz * sb3 + childz * sb2) * year

> sdre <- rbind(attr(VarCorr(gcp)$schoolid,"stddev"), + attr(VarCorr(gcp)$childid,"stddev")) > xb <- predict(gcp, newdata = pdf, re.form = ~0) > pdf <- mutate(pdf, + xb = xb, + fv = xb + + (schoolz * sdre[1,1] + childz * sdre[2,1]) + + (schoolz * sdre[1,2] + childz * sdre[2,2]) * year + )

We are now ready to plot the results

. label define schoolz -1 "-1 school sigs" 1 "+1 school sigs" . label values schoolz schoolz . label define lowinc 50 "50% poor" 100 "100% poor" . label values lowinc lowinc . gen eth = 1 + black + 2*hispanic . sort schoolz lowinc childz eth year . twoway line fv year if eth == 1, c(asc) lpat(solid) /// > || line fv year if eth == 2, c(asc) lpat(shortdash) lc(red) /// > || line fv year if eth == 3, c(asc) lpat(dash) lc(blue) /// > , by(lowinc schoolz, legend(off) /// > note("Legend: solid=white, long-dash=hispanic, short-dash=black") /// > title("Growth Curves by School and Child Characteristics")) /// > xtitle(year) ytitle(math) . graph export egmfig3.png, width(700) replace (file egmfig3.png written in PNG format)

> ggplot(filter(pdf,childz==-1), aes(year, fv, linetype=eth)) + + geom_line() + + geom_line(data=filter(pdf,childz==1), aes(year, fv, linetype=eth)) + + facet_grid(lowinc ~ schoolz) > ggsave("egmfig3r.png", width=10, height=8, dpi=72)

The plot shows the relative importance of fixed and random effects at the school and individual levels. At the school level the difference between 50% and 100% poor is smaller than the two-sigma difference in unobserved school characteristics. At the individual level the differences between whites, blacks and hispanics are dwarfed by the two-sigma differences in unobserved child characteristics. Interestingly, within each school type the average line for hispanics starts at about the average level of blacks in first grade, but catches up with the average level of whites by sixth grade.

Continue with Logit Models