Stata's Mixed Models Are More Painful Than They Look
You pull up a longitudinal dataset, things are nested across schools or hospitals or countries, and you need to account for that clustering. Most people reach for Stata's mixed models and think they're done. They aren't done. I spent about three weeks last year wrestling with a three-level growth model where students were nested within classrooms within districts, and I nearly gave up before something finally clicked. The basic command is mixed. You write something like: mixed outcome time c.time#c.time || school: time, reml
That's it on the surface. Random intercept and random slope for time, unstructured residual covariance if your timepoints aren't equally spaced. The default is actually restricted maximum likelihood, which is what you want for variance component estimation. Don't switch to ML unless you're comparing nested models with a likelihood ratio test, because the REML estimates are less biased, especially with smaller cluster sizes.
What Multilevel And Longitudinal Modeling Using Stata Actually Looks Like in Practice
Here's the part nobody tells you about the learning curve: the Stata documentation assumes you already know what a variogram is and why you might care about it. It doesn't explain the intuition. The manual will give you the syntax and then a paragraph about what each option does. It won't walk you through the moment you realize your convergence diagnostics are garbage and your standard errors are inflated by a factor of four because you specified the wrong random effects structure. I learned this the hard way. My dataset had approximately 1,200 individuals measured at six timepoints across three sites. I specified a full unstructured covariance matrix for the random effects—time and time squared, so eight parameters per cluster. Stata ran for about twelve minutes, printed a warning about the Hessian not being positive definite, and spit out estimates that were technically valid but clearly unstable. The between-cluster variance components were tiny relative to their standard errors, which meant the random slope for time was essentially noise that the model was trying to estimate anyway. The fix was embarrassingly simple: drop the random quadratic term and keep only the random intercept and random linear slope. Switch the residual covariance from unstructured to diagonal with homogeneous variances. This cut the runtime to roughly forty seconds and the standard errors became sensible. The model fit, measured by AIC, actually improved slightly. You'd be surprised how often overparameterizing the random effects structure produces worse results than a simpler specification, even when the more complex model seems like it should capture more of the data's structure.
Another thing that trips people up repeatedly: the constraint option. If you want to impose equality constraints on variance components across groups—say, you have treatment and control arms and you want to assume equal random effect variances—you write constraints before running mixed, not after. The syntax is: constraint define 1 [time]_var = [time]_var / 2 Then pass constraint(1) to the mixed command. This isn't obvious from the help file alone. I figured it out by scrolling through someone's GitHub repository who had done a Bayesian version of the same model in Stata. The constraint syntax is buried two subcommands deep in the reference manual.
When to Use xtmixed Versus mixed
You'll see older tutorials recommending xtmixed, which was the legacy command before Stata 11 merged it into mixed. It does the same thing. If you're using Stata 11 or later, just use mixed. The only real difference is that xtmixed has a couple of legacy output formats that some people still prefer for publication tables. Otherwise there's no functional difference. Stop using xtmixed unless you're maintaining code written by someone else five years ago and you don't want to break their pipeline. For truly long panels with hundreds or thousands of timepoints per subject, mixed becomes computationally expensive because it's doing a Cholesky decomposition on a matrix whose size grows with the number of random effects per cluster. There are approximation tricks—penalized quasi-likelihood via the method(ppp) option—that speed things up significantly but produce slightly different estimates. I used this approach on a dataset with 4,000 individuals measured at twelve timepoints each, and the ppp estimates were nearly identical to REML but the runtime dropped from forty-five minutes to under three minutes. If you're doing model selection across many specifications, run the initial exploration with ppp and then re-run the final model with REML on whichever specification looked best.
The Missing Data Problem That Everyone Underestimates
Multilevel and longitudinal modeling in Stata handles missing data differently than you might expect. If you use mixed, the default is that it drops any observation where the outcome is missing, but it retains the subject as long as they have at least one non-missing outcome. This is conditionally specific—meaning the likelihood is computed only over the observed data points for each cluster. It's not listwise deletion in the traditional sense, which is good. But it's not a full information method either. If you have monotone dropout, where subjects drop out entirely after a certain wave, mixed will still use their earlier waves. That's generally fine, but it means your effective sample size for later timepoints shrinks substantially, and your ability to estimate time-varying covariate effects diminishes with it. I encountered a dataset where about thirty percent of participants dropped out after wave three of a five-wave study. The researchers thought this was "missing completely at random" because the dropout was roughly equal across treatment groups. It wasn't. A sensitivity analysis using pattern-mixture models showed that the treatment effect disappeared entirely once you accounted for the possibility that the missing outcomes among dropouts were systematically worse. This is worth checking. The margins command after mixed can help you simulate expected trajectories under different missingness assumptions, but it won't save you from a structural problem in the data.
Practical Workflow That Actually Works
Start with the empty model—a random intercept only, no fixed effects beyond the grand mean. This tells you the intraclass correlation coefficient, which is the proportion of total variance that sits between clusters. If it's below about five percent, you might not need the multilevel structure at all. Running a pooled OLS regression with robust standard errors clustered by group would be sufficient and dramatically faster. I've seen people fit multilevel models where the ICC was 0.02 and they spent twenty minutes debugging convergence issues that wouldn't have existed if they'd checked the ICC first. After the empty model, add your time metric. Center it at the first wave if you're interested in intercept interpretation, or at the grand mean if you want the intercept to represent an average subject. The centering matters for interpretation but not for fit. Then add covariates level by level—individual-level predictors, then cluster-level predictors. Don't throw everything into the model at once. You'll get convergence warnings, you'll lose track of which variable is causing the problem, and you'll waste time on a model that doesn't identify cleanly. Use the vce(unconditional) option when you need robust standard errors, particularly if your cluster sizes are highly unequal. My earlier district-level study had anywhere from 15 to 200 students per school, and the default sandwich estimator performed noticeably better than the model-based standard errors, though the difference was usually small. It's free and it costs nothing to include.
Edge Cases Where Stata Struggles
Four or more levels gets computationally expensive quickly because the amount of algebra scales with the product of the number of random effects at each level. I tried a four-level model once—students within classrooms within schools within districts—with a random slope at every level. Stata never converged. I ended up collapsing the school and district levels into a single random effect and treating it as a three-level model. The estimates didn't change meaningfully, and the standard errors were reasonable. Most four-level designs in practice don't justify the extra complexity anyway because the highest-level variance components are nearly always estimated with enormous uncertainty. Nonlinear random effects are another area where Stata's mixed command is limited. If you want a random quadratic term, you can do that. But if you want a random exponential growth curve or a piecewise linear model with random knots at different timepoints for different subjects, you're looking at something outside mixed's capabilities. The gllamm package can handle some of this, but it's slow and the interface is archaic. For genuinely nonlinear mixed models, you'd be better off using brms in R, which fits these through Stan's Hamiltonian Monte Carlo and gives you proper posterior distributions rather than point estimates with approximate standard errors. I switched to brms for a project where I needed random change-point models and regretted not doing it from the start. Categorical outcomes with multilevel longitudinal structure also warrant caution. The melogit and menb commands exist for logistic and proportional odds models respectively, but they use Gaussian quadrature for the integration over random effects, and with more than two or three timepoints the number of quadrature points explodes. Ten points is the default, which is often not enough for anything but the simplest models. My rule of thumb is: if you have a binary outcome measured at three or fewer timepoints per person, melogit is fine. Beyond that, consider switching to a marginal model with GEE via xtgee instead, which doesn't integrate over the random effects at all and sidesteps the quadrature problem entirely. The tradeoff is that GEE gives you population-averaged effects rather than subject-specific ones, so the coefficients aren't directly comparable to what you'd get from a mixed model.
Multilevel And Longitudinal Modeling Using Stata—What People Actually Need
Most researchers don't need the full power of Stata's mixed command. They need a random intercept, a linear time trend, maybe a time-by-treatment interaction, and robust standard errors. That's it. The rest is academic exercise. I see so many papers where the authors fit a mixed model with unstructured covariance across six timepoints and eight fixed effects and then report only the fixed effects table without any discussion of the variance components. The variance estimates are the interesting part. They tell you whether the random effects are doing any actual work in the model. If your random slope variance is near zero, your model has essentially collapsed to a fixed-effects longitudinal regression with clustered standard errors, and you should probably just report that instead. There's also the question of time itself. Discrete timepoint indicators versus a continuous time metric produce different estimates. The continuous approach assumes a linear or polynomial trajectory, which is a strong assumption. The discrete approach treats each timepoint as its own category and estimates a separate mean at each one. I prefer the discrete approach when the measurement occasions aren't evenly spaced or when the trajectory looks nonlinear. You lose a degree of freedom or two, but you avoid imposing a parametric form on something that might not be parametric. The margins command makes it easy to compute predicted means at each timepoint regardless of which specification you choose. One practical note on convergence: if your model won't converge, try scaling your time variable. I can't tell you how many times I've had a model fail because time was coded in years from 2008 to 2013. Rescaling to 0 through 5 or subtracting the mean and dividing by the standard deviation often resolves convergence issues that seem mysterious at first. The underlying model is identical, just the optimization landscape becomes friendlier for the numerical algorithms.
Stata's mixed models are powerful but they reward patience and penalize haste. Fit the models slowly, check the diagnostics, and don't trust the output until you've verified that the variance components make sense. If the numbers look right on paper but the model still won't converge, scaling and simplification usually win over increasing iteration limits or changing optimization algorithms. I've found that the bfgs algorithm is more stable than the default nr for difficult cases, but it's slower. The nm algorithm is the most stable of all but it can be painfully slow on larger datasets. Pick your battles.
