The new meintreg command fits models in which the outcome is interval measured (interval-censored) and the observations are clustered.

Interval measured means that rather than the outcome (y) being observed precisely, it is known only that yl ≤ y ≤ yu in some or all observations. Observations can also be left-censored (y ≤ yl) or right-censored (y ≥ yu). An interval-measured outcome might be income, recorded in income brackets, or weekly minutes of exercise, recorded as less than 30 minutes, 31-59 minutes, 60-89 minutes, etc.

Multilevel mixed effects means that the fitted model accounts for clustering, such as when people live near each other or students attend the same school or students are tested repeatedly.

Let’s see it work

We have fictional data on 8,424 people living in 50 cities. Different cities have different policies. One aspect of this is that some cities are more walker-friendly than others. Our data include a binary variable for being walker-friendly.

In these data, respondents were asked how many minutes per week they engage in physical activity outside of their jobs and were offered the following response categories: 0–30, 30–59, 60–89, 90–119, 120–239, and more than 240 minutes. This variable will be our dependent variable. It is both interval- and right-censored. The dataset contains variables exerlo and exerup. The values recorded in the variables are

 Walking per week exerlo exerup Meaning 0 30 0-30 minutes 30 59 30-59 minutes 60 89 60-89 minutes 90 119 90-119 minutes 120 239 120-239 minutes 240 . more than 240 minutes

Our independent variables will be

 walk walker-friendly city age respondent’s age in years work hours/week worked kids children under age 12 in household

Even so, we expect that there are unmeasured characteristics about cities that matter—like weather, for example. Unmeasured characteristics would result in people who live in the same city behaving in overly similar ways, and so we will allow for city-level random effects. Variable cid records city ID number.

To fit a random-intercept model using exercise in the range exerlo and exerup, we type

. meintreg exerlo exerup age work kids walk || cid:

Fitting fixed-effects model:

Iteration 0:   log likelihood = -16228.837
Iteration 1:   log likelihood = -16211.002
Iteration 2:   log likelihood = -16209.922
Iteration 3:   log likelihood =  -16209.92
Iteration 4:   log likelihood =  -16209.92

Refining starting values:

Grid node 0:   log likelihood =  -16183.23

Fitting full model:

Iteration 0:   log likelihood =  -16183.23
Iteration 1:   log likelihood = -16021.782
Iteration 2:   log likelihood = -15902.566
Iteration 3:   log likelihood = -15856.315
Iteration 4:   log likelihood = -15836.069
Iteration 5:   log likelihood =  -15830.93
Iteration 6:   log likelihood =  -15830.08
Iteration 7:   log likelihood = -15830.052
Iteration 8:   log likelihood = -15830.052

Mixed-effects interval regression               Number of obs     =      8,424
Uncensored     =          0
Left-censored  =          0
Right-censored =         54
Interval-cens. =      8,370

Group variable:             cid                 Number of groups  =         50
Obs per group:
min =         81
avg =      168.5
max =        249

Integration method: mvaghermite                 Integration pts.  =          7

Wald chi2(4)      =     125.15
Log likelihood = -15830.052                     Prob > chi2       =     0.0000
 Coef. Std. Err. z P>|z| [95% Conf. Interval] age -.1745526 .0695285 -2.51 0.012 -.3108259 -.0382792 work -.3635766 .0794423 -4.58 0.000 -.5192807 -.2078726 kids -10.80659 1.102927 -9.80 0.000 -12.96829 -8.644894 walk 13.7689 6.514321 2.11 0.035 1.001061 26.53673 _cons 78.62265 5.01203 15.69 0.000 68.79925 88.44605 cid var(_cons) 269.5222 56.84988 178.2592 407.5089 var(e.exerlo) 2440.81 40.2559 2363.172 2520.999 LR test vs. interval model: chibar2(01) = 759.74 Prob >= chibar2 = 0.0000

Just above the coefficient table is a header. It reports the number of censored and uncensored observations, reports the number of observations per city, and provides information about the fitted model.

The coefficient table itself is divided into three or four parts. The first part reports the fixed part of the model, that is, the effects of variables age, work, etc.

The second part reports the variance of the random intercept. We fit a model with a random intercept (_cons) whose estimated mean is 79 and standard deviation is 16 (a variance of 270).

The third part, var(e.exerlo), reports the variance of the model’s error. Stata labels errors as e. followed by the dependent variable’s name. Even so, e.exerlo is a bit misleading. A better name would have been e.exer, with exer being the imaginary name of the unobserved variable that is bounded by exerlo and exerup.

The fourth part—the part below the table—reports a test of whether the intercept is random. Anyone would reject the null hypothesis that it has variance equal to zero.

To obtain standard deviations instead of variances for the random portion of the model, we can type estat sd:

. estat sd

 Coef. Std. Err. z P>|z| [95% Conf. Interval] cid sd(_cons) 16.41713 1.731419 13.35137 20.18685 sd(e.exerlo) 49.40456 .4074108 48.61246 50.20955

We can also obtain the residual intraclass correlation.

. estat icc

Residual intraclass correlation

 Level ICC Std. Err. [95% Conf. Interval] cid .0994425 .0189461 .0679834 .1432221

Tell me more

You can also fit Bayesian multilevel interval regression using the bayes prefix.