It is common to fit a model where a variable (or variables) has an effect on the
expected mean. You may also want to fit a model where a variable has an effect
on the variance, that is a model with heteroskedastic errors.
One example is a hierarchical model where the error variance is not constant
over some categorical variable (e.g., gender). Another example is a longitudinal model where you might want to allow the errors
to vary across
time. This page describes how
to fit such a model and demonstrates with an example.

Example 1: Heteroskedastic errors by group

It may be useful to note that a multi-level model (or a
single level model) that allows different error terms between groups is similar
to a multi-group structural equation model in which all parameters in the model
are constrained to equality except for the error terms.

The dataset for this example includes data on 7185 students in 160
school. This dataset and model come from the HLM manual. The variable mathach is a measure of each student’s
math achievement, the variable female is a binary variable equal to one if the student is
female and zero otherwise, and the variable id contains the school
identifier.

Below we run a model that uses female to predict mathach. This
model assumes equal error variances for males and females.

use https://stats.idre.ucla.edu/stat/stata/faq/hsb, clear

mixed mathach female || id:, var mle nolog

Mixed-effects ML regression                     Number of obs      =      7185
Group variable: id                              Number of groups   =       160

                                                Obs per group: min =        14
                                                               avg =      44.9
                                                               max =        67


                                                Wald chi2(1)       =     62.89
Log likelihood =  -23526.66                     Prob > chi2        =    0.0000

------------------------------------------------------------------------------
     mathach |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      female |   -1.35939   .1714111    -7.93   0.000     -1.69535    -1.02343
       _cons |   13.34526    .253935    52.55   0.000     12.84756    13.84297
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   8.108982   1.018274      6.339834    10.37181
-----------------------------+------------------------------------------------
               var(Residual) |   38.84481   .6555316      37.58101    40.15112
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) =   936.66 Prob >= chibar2 = 0.0000

In order to fit a model with heteroskedastic errors, we must modify the above model. The model shown above can be written as:

mathach_ij = b0 + b1*female_ij + u_i + e_ij                (1)

This model assumes:

e_ij ~ N(0,s2)
u_i ~ N(0,t2)

Where e_ij is the level 1 errors (i.e., residuals), and u_i is the random intercept across classrooms (i.e., the level 2
variance). We want to allow
the variance of e_ij (i.e., s2) to be different
for male and female students. One way to think about this would be to replace the single
e_ij
term, with
separate terms for males and females, in equation form:

e_ij = e(m)_ij*male + e(f)_ij*female                          (2)

Where male is a dummy variable equal to one if the student is male and
equal to zero otherwise, and e(m)_ij is the error term for males. The
variable female and the term e(f)_ij are the same except that they
are for females. Another way to write this is:

e_ij ~  N(0,s2_1)  males
        N(0,s2_2)  females

Another way to think about this is to use one
group as the reference group, and estimate the difference between the variance
of the two groups. In equation form this would be:

e_ij = e(m)_ij + e(f)_ij*female                                (3)

Where e(m)_ij is the error term for males, and e(f)_ij is the
difference between the errors for males and females. This can also be written:

e_ij ~ N(0,s2 + s2_2*female)

This second way of looking
at the error structure is part of how we will restructure our model to allow for
heteroskedastic error terms.

Equation 3 says that the residual variance is a function of gender. So we can
write it as r_ijk, where i is the group, j  is the subject, and k is the
gender.

mathach_ij = b0 + b1*female + u_i + r_ijk               (4)

To model this, we add a level to our model. Now
the third level will be classrooms (previously level 2), the second level will be
students (previously level 1), and level 1 will be a single case within each
student.
The only random effect at level 1 is gender (even the intercept is fixed). The new model can be written as:

mathach_ij = b0 + b1*female_ij + u_i + e_ij*female + r_ij0
r_ijk =   r_ij0   male 
          r_ij1   female=e_ij+r_ij0

Where the level one error
is represented by r_ij0. Because males are the omitted category in the
random portion of the model, the variance of r_ij0 is the variance of the
errors males (i.e., female=0). This is similar to the interpretation of
the intercept in the fixed portion of the mode, where b0 represents the
intercept for males, since they are the omitted category. The difference between
the error terms for males and females is e_ij so the error term for females is
r_ij + e1_ij.

A final consideration before we can actually fit the model is
that Stata restricts the estimates of the error variances to be greater than zero. As a result, we must use the
group with the smallest variance as the omitted category, since this will result in the values of the error variances for the
group(s) being positive. To establish which group has the lowest variance
you can examine the residuals from a model without heteroskedastic errors. In this example,
examination of the residuals shows that the residual variance is smaller for females than
males.

The code below creates two new variables, male, and nid. The variable nid
will be the identifier for the students.

recode female (0=1) (1=0), gen(male)
gen nid = _n

The resulting data set should look something like the one shown below. Where each classroom
(identified by the variable id) has multiple rows in the dataset, but each case
(i.e., each student, identified by the variable nid) has only one row.

      id   nid  
    1224     1  
    1224     2  
    1224     3  
     ...
    1288    48  
    1288    49  
    1288    50  
     ...
    1296    73  
    1296    74  
    1296    75

Once the necessary variables are created, we can run the model as shown
below, which allows for a difference in the variance of the errors for males and
females. It is necessary to specify the nocons
option suppresses the random intercept at level 2, so that the only random
effect at level 2 is gender (i.e., male).

mixed mathach female || id: || nid: male, nocons var mle nolog

Mixed-effects ML regression                     Number of obs      =      7185

-----------------------------------------------------------
                |   No. of       Observations per Group
 Group Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
             id |      160         14       44.9         67
            nid |     7185          1        1.0          1
-----------------------------------------------------------

                                                Wald chi2(1)       =     63.03
Log likelihood = -23522.932                     Prob > chi2        =    0.0000

------------------------------------------------------------------------------
     mathach |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      female |  -1.363969   .1718025    -7.94   0.000    -1.700695   -1.027242
       _cons |   13.34707   .2548619    52.37   0.000     12.84755    13.84659
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   8.090527   1.016457      6.324639    10.34946
-----------------------------+------------------------------------------------
nid: Identity                |
                   var(male) |   3.622413   1.333432       1.76062    7.452988
-----------------------------+------------------------------------------------
               var(Residual) |   37.13827   .8657943      35.47953    38.87456
------------------------------------------------------------------------------
LR test vs. linear regression:       chi2(2) =   944.12   Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

The second table in the output shows the estimated random effects. The residual variance for females is equal to var(Residuals) = 37.138,
while the variance for males is var(Residuals) + var(male) =
37.1383 + 3.622 = 40.7607. Since the 95% confidence interval for var(male)
does not include zero, we can say that the difference between the variances is
statistically significant at the p<0.05 level.

Traditionally, and in some other packages (e.g., HLM), the variances in models with heteroskedastic errors have been parameterized as:

s2_ij = exp(a_0 + a_1*x_ij)

Where s2_ij is the variance of the residuals, x_ij is a
dichotomous variable equal to one for one group and zero for the other, and a_0 and a_1 are estimated parameters for the
variance. a_0 is the variance in the group where x_ij
= 0, and a_1 is the difference in the variances between the two groups.
If you want to compare the results of your analysis in Stata to those from another
package, it is possible to translate between the two. As a first step, let’s
think about what Stata is estimating a little more. One way to write what Stata models
is:

s2_ij = var(Residuals) + var(x_ij)*x_ij

Where var(Residuals) is the variance of the level 1 errors, and
var(x_ij) is the random effect of the dummy variable x_ij. Because Stata models the natural log of the standard deviation of the error term,
the above is visually clear, but not quite correct. The actual model is:

s2_ij = exp(2*lnsd_0 ) + exp(2*lnsd_1)*x_ij

Where lnsd_0 is the natural log of the standard deviation of the level
1 errors, and lnsd_1 is the natural log of the standard deviation of
the level 2 random effect. We can convert the estimates Stata gives us to a_0
and a_1 using the formulas below (each formula is written twice, the
first for visual clarity, the second to reflect what is actually modeled by
Stata):

a_0 = ln(var(Residuals))
a_0 = ln(exp(2*lnsd_0))

a_1 = ln(var(Residuals)+var(x_ij)) - ln(var(Residuals))
a_1 = ln[exp(2*lnsd_0) + exp(2*lnsd_1)] - ln[exp(2*lnsd_0)]

We can use the display command to calculate these values. Stata stores the ln(sd) for the
level 1 residuals as [lnsig_e]_cons, and the ln(sd) of the random effect
of male in [lns2_1_1]_cons. Below we use these values to calculate a_0 and a_1.

di "a_0=" ln(exp(2 * [lnsig_e]_cons))
a_0=3.6147746

di "a_1=" ln(exp(2 * [lnsig_e]_cons)+exp(2 * [lns2_1_1]_cons)) - ln(exp(2 * [lnsig_e]_cons))
a_1=.09308814

Example 2: A growth curve model

When you run a growth curve model in the structural equation modeling
framework, it is not uncommon to allow the error variance to be different across
time points, or at least to test whether such differences exist. A longitudinal
model that allows different error variances across time points is similar to a
growth model in the structural equation modeling setting, where all parameters
except for the error terms across time points are constrained to equality. Below we show
how to allow for different error variances across time points, and how to test
whether the error variances are significantly different from each other.

The example dataset contains data on 239 subjects. The dataset includes
observations at up to five time points per subject, for a total of 1079 observations.
The variable id uniquely identifies subjects. The time at which an
observation was taken is indicated by the variable time, which takes on
the values 0 to 4.

Below we use mixed to fit a model where
the variables x and time
predict the variable attit. This model assumes the variance is the same
across time points.

use https://stats.idre.ucla.edu/stat/stata/faq/nys2, clear

mixed attit x time ||id:, var mle nolog

Mixed-effects ML regression                     Number of obs      =      1079
Group variable: id                              Number of groups   =       239

                                                Obs per group: min =         1
                                                               avg =       4.5
                                                               max =         5


                                                Wald chi2(2)       =    324.65
Log likelihood =  140.31043                     Prob > chi2        =    0.0000

------------------------------------------------------------------------------
       attit |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |   .0241245   .0033107     7.29   0.000     .0176357    .0306134
        time |   .0600884   .0039557    15.19   0.000     .0523353    .0678415
       _cons |   .1191882   .0184893     6.45   0.000     .0829498    .1554267
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   .0308536    .003533      .0246511    .0386167
-----------------------------+------------------------------------------------
               var(Residual) |   .0311131   .0015204      .0282714    .0342405
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) =   324.66 Prob >= chibar2 = 0.0000

The model above can be written as:

attit_it = b0 + b1*x_it + b2*time_it + u_i + e_it               (4)

It is assumed that:

e_it ~ N(0,s2)
u_i ~ N(0,t2)

Where u_i is the random intercept across individuals (i.e., the level 2
variance), and e_it represents the level 1 errors (i.e., residuals). We want to allow
the variance of e_it to be different
at each time point. One way to think about this would be to replace the single
e_ij
term, with
separate terms for each time point, in equation form:

e_it = e_i0*t0 + e_i1*t1 + e_i2*t2 + e_i3*t3 + e_i4*t4

Where t0 is a dummy variable equal to one if time = 0 (i.e. the
first time point) and equal to zero otherwise, and e_i0 is the error for
the first measurement occasion (time=0). The variables t1 to t4, are dummy variables for time
points 1 to 4. The terms
e_i1 to e_i4 are the error terms for time points 1 to 4.

In order to
model the heteroskedastic errors, we add a third level to our model. In this new model,
the third level will be individuals (previously level 2), the second level will be time
points (previously level 1), and level 1 will be a single case within each time point.
Since the effect of time is in the level at model 2, only random effects for
time are included at level 1. The new model can be written as:

attit_itk = b0 + b1*x_it + b2*time_it + e_i1*t1 + e_i2*t2 + e_i3*t3 + e_i4*t4 + r_ij0

Where the level one error
is represented by r_it0. Because the time=0 is the omitted category, the variance of
r_it0 is the variance of the
errors at time=0. The
random effects e_i1 to e_i4 represent the difference between the variance
of the errors and the variance of the errors at time=0 for each time
point.

A final consideration before we can actually fit the model is
that Stata restricts the estimates of the error variances to be greater than zero. As a result, we must use the time point (or other
category/group) with the smallest variance as the omitted category, since this will result in the values of the error variances for the other time points
being positive. To establish which time point or group has the lowest variance
you can examine the residuals from a model without heteroskedastic errors. In
this example, the residual variance is smallest when time=0.

In order to fit the model we will need dummy variables for each of the other
time points.

gen t1 = (time==1)
gen t2 = (time==2)
gen t3 = (time==3)
gen t4 = (time==4)

Now we are ready to fit our model using mixed. The nocons option
suppresses the random intercept at level 2.

mixed attit x time ||id: ||time: t1 t2 t3 t4, nocons var mle nolog

Mixed-effects ML regression                     Number of obs      =      1079

-----------------------------------------------------------
                |   No. of       Observations per Group
 Group Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
             id |      239          1        4.5          5
           time |     1079          1        1.0          1
-----------------------------------------------------------

                                                Wald chi2(2)       =    319.37
Log likelihood =  143.67787                     Prob > chi2        =    0.0000

------------------------------------------------------------------------------
       attit |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |   .0238029   .0033059     7.20   0.000     .0173235    .0302823
        time |   .0597567   .0039638    15.08   0.000     .0519878    .0675256
       _cons |      .1208   .0178047     6.78   0.000     .0859035    .1556965
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   .0280981   .0034551      .0220805    .0357557
-----------------------------+------------------------------------------------
time: Independent            |
                     var(t1) |   .0027798   .0046901      .0001018    .0758888
                     var(t2) |   .0058393   .0053178      .0009799    .0347972
                     var(t3) |   .0125841   .0060606      .0048964    .0323421
                     var(t4) |     .01354   .0059784      .0056988    .0321702
-----------------------------+------------------------------------------------
               var(Residual) |   .0247609   .0033388      .0190103     .032251
------------------------------------------------------------------------------
LR test vs. linear regression:       chi2(5) =   331.40   Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

Once we have run the model, we may want to test whether the estimated
variances are different either from zero or different from each other. Below nlcom is used to
get the estimated variance at each time point.

nlcom   (time0: exp(2 * [lnsig_e]_cons)) ///
        (time1: exp(2 * [lnsig_e]_cons) + exp(2 * [lns2_1_1]_cons)) ///
        (time2: exp(2 * [lnsig_e]_cons) + exp(2 * [lns2_1_2]_cons)) ///
        (time3: exp(2 * [lnsig_e]_cons) + exp(2 * [lns2_1_3]_cons)) ///
        (time4: exp(2 * [lnsig_e]_cons) + exp(2 * [lns2_1_4]_cons))

       time0:  exp(2 * [lnsig_e]_cons)
       time1:  exp(2 * [lnsig_e]_cons) + exp(2 * [lns2_1_1]_cons)
       time2:  exp(2 * [lnsig_e]_cons) + exp(2 * [lns2_1_2]_cons)
       time3:  exp(2 * [lnsig_e]_cons) + exp(2 * [lns2_1_3]_cons)
       time4:  exp(2 * [lnsig_e]_cons) + exp(2 * [lns2_1_4]_cons)

------------------------------------------------------------------------------
       attit |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       time0 |   .0247609   .0033388     7.42   0.000      .018217    .0313048
       time1 |   .0275406   .0034705     7.94   0.000     .0207385    .0343427
       time2 |   .0306002   .0035976     8.51   0.000     .0235491    .0376513
       time3 |    .037345   .0044345     8.42   0.000     .0286536    .0460364
       time4 |   .0383009   .0044824     8.54   0.000     .0295155    .0470862
------------------------------------------------------------------------------

We can also use nlcom to estimate the difference in the variance between time points,
in this case, the difference between the variance at time = 1 and time = 4.

nlcom (t4_t1: exp(2 * [lns2_1_4]_cons) - exp(2 * [lns2_1_1]_cons)) 

       t4_t1:  exp(2 * [lns2_1_4]_cons) - exp(2 * [lns2_1_1]_cons)

------------------------------------------------------------------------------
       attit |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       t4_t1 |   .0107602   .0060541     1.78   0.076    -.0011055    .0226259
------------------------------------------------------------------------------

References

Gutierrez, Roberto G. (2008) Tricks of the Trade: Getting the most out of
xtmixed. 2008 Fall North American Stata Users Group Meeting, San Francisco, CA.

Raudenbush, Stephen, Anthony Bryk, Yuk Fai Cheong, and Richard Congdon
(2001) HLM 5: Hierarchical Linear and Nonlinear Modeling. Scientific
Software International. Lincolnwood, IL.