Germán Rodríguez
Multilevel Models Princeton University

3-Level Models

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.

Variance Components

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 http://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("http://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.

Growth Curves

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.

.  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.

Growth Curve Predictors

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 a school with no poor students!

. 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.

Plotting Illustrative Growth Curves

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.