What’s this about?
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.
Learn more about Stata’s multilevel mixed-effects models features.
Read more about multilevel interval regression in the Stata Multilevel Mixed-Effects Reference Manual.
