StataNow spotlight: mixed across space and time
StataNow™ introduces a new type of residual structure to the mixed command: spatial structures. Spatial error covariance structures let you model correlation between errors as a function of distance computed from coordinate variables (Cressie 2015). You can also use these structures to account for distance across time. This is particularly useful for intensive longitudinal data, ecological momentary assessment data, or any longitudinal data with unequally spaced time.
Let’s take a look.
The data
We begin by loading a dataset that tracks children’s weight using doctor-visit data from their first 2.5 years of life, taken from Goldstein (1986).
. use https://www.stata.com/exdata/childwt.dta
(Weight data on Asian children)
. describe
Contains data from https://www.stata.com/exdata/childwt.dta
Observations: 266 Weight data on Asian children
Variables: 4 3 Nov 2026 06:37
(_dta has notes)
| Variable Storage Display Value | ||
| name type format label Variable label | ||
| id int %8.0g Child identifier age float %4.2f Age (years) weight float %5.2f Weight (kg) girl byte %9.0g bg Female | ||
Listing the first 11 observations reveals the data structure.
. list in 1/11, sepby(id)
| id age weight girl | |
| 1. | 45 0.00 4.14 Boy |
| 2. | 45 0.14 5.17 Boy |
| 3. | 45 0.66 10.86 Boy |
| 4. | 45 1.22 13.15 Boy |
| 5. | 45 1.43 13.20 Boy |
| 6. | 45 2.27 15.88 Boy |
| 7. | 258 0.00 3.15 Girl |
| 8. | 258 0.19 5.30 Girl |
| 9. | 258 0.69 9.74 Girl |
| 10. | 258 1.13 9.98 Girl |
| 11. | 258 2.31 11.34 Girl |
Because the age variable is not measured in integers, we cannot xtset these data with age. However, we can still visualize weight changes over time with xtline by explicitly specifying the panel ID and time variables.
. xtline weight, i(id) t(age) overlay legend(off)
;)
Based on this graph, we suspect a quadratic growth rate. We also see more variability between children as they get older.
The model
To analyze change in weight over time, we use a multilevel mixed-effects model. We estimate fixed-effects coefficients for age, the square of age, and their interactions with sex to fit a quadratic growth trajectory, allowing average trajectories to differ by sex. In addition, we estimate random-effects parameters for intercept and linear slope to account for individual differences in growth trajectories.
. mixed weight c.age##c.age##girl || id: age, covariance(unstructured)
Performing EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log likelihood = -305.5123
Iteration 1: Log likelihood = -305.44659
Iteration 2: Log likelihood = -305.44647
Iteration 3: Log likelihood = -305.44647
Computing standard errors ...
Mixed-effects ML regression Number of obs = 266
Group variable: id Number of groups = 68
Obs per group:
min = 2
avg = 3.9
max = 6
Wald chi2(5) = 3820.38
Log likelihood = -305.44647 Prob > chi2 = 0.0000
| weight | Coefficient Std. err. z P>|z| [95% conf. interval] | |
| age | 8.696155 .2272288 38.27 0.000 8.250795 9.141515 | |
| c.age#c.age | -1.971854 .0915805 -21.53 0.000 -2.151349 -1.792359 | |
| girl | ||
| Girl | -.1705086 .1524447 -1.12 0.263 -.4692947 .1282775 | |
| girl#c.age | ||
| Girl | -.9934062 .3236593 -3.07 0.002 -1.627767 -.3590456 | |
| girl#c.age# | ||
| c.age | ||
| Girl | .289226 .1308179 2.21 0.027 .0328278 .5456243 | |
| _cons | 3.298761 .1074758 30.69 0.000 3.088112 3.509409 | |
| Random-effects parameters | Estimate Std. err. [95% conf. interval] | |
| id: Unstructured | ||
| var(age) | .240249 .0687079 .1371607 .4208173 | |
| var(_cons) | .1713331 .0610408 .0852283 .3444281 | |
| cov(age,_cons) | .1122623 .0438132 .02639 .1981346 | |
| var(Residual) | .3157933 .0383252 .2489428 .4005957 | |
We have strong evidence for a quadratic growth trajectory in which the linear and quadratic terms differ by sex but not the intercept.
After we account for age and sex, weight varies more within children than between children. We also see that there is more variability in slope than in intercept between children.
To see the within-child correlation matrix that this model implies, we use the estat wcorrelation postestimation command. By default, it shows the first child.
. estat wcorrelation Standard deviations and correlations for id = 45: Standard deviations:
| obs | 1 2 3 4 5 6 | |
| sd | 0.698 0.723 0.859 1.057 1.140 1.496 |
| obs | 1 2 3 4 5 6 | |
| 1 | 1.000 | |
| 2 | 0.370 1.000 | |
| 3 | 0.409 0.454 1.000 | |
| 4 | 0.418 0.476 0.632 1.000 | |
| 5 | 0.417 0.479 0.645 0.736 1.000 | |
| 6 | 0.408 0.478 0.668 0.777 0.802 1.000 |
Within-child standard deviations increase over time, and time points farther apart generally have lower correlations.
Have we fully captured the autocorrelation of the data? This is important because we know that we need to correctly specify the level-1 residual structure for our inferences to be valid (for example, Ferron, Dailey, and Yi [2002]). We cannot try adding an autoregressive structure because our time variable is not an integer; we would need to do some data management first.
A more flexible option is to add a spatial residual structure.
Spatial error covariance structures
In mixed, there are four options for spatial correlation structures:
- Exponential (spexponential)—exponential decay of spatial correlation
- Gaussian (spgaussian)—smooth Gaussian (squared-exponential) correlation
- Spherical (spspherical)—compact-support correlation with a finite range
- Power (sppower)—power-law decay
Each of these requires a coordinates() suboption in which you specify geographic coordinates used to compute distances. With longitudinal data, you specify the time variable as the single coordinate.
We can visualize each structure by graphing correlations across a two-year distance.
;)
These graphs show the correlation across time differences (d) at different values of their spatial-scale parameter, under each spatial structure. Notice that in the case of one coordinate, as in longitudinal data, the exponential and power-law decay structures are identical with different parameterizations.
We are not sure which of these is most appropriate for our dataset, so let’s fit a model with each, store the estimates, and compare their fit by using AIC and BIC.
We first rerun and store our model estimates with the default, independent structure. For the Gaussian structure, we specify a starting value in the phi0() suboption to help with convergence.
. quietly mixed weight c.age##c.age##girl || id: age, covariance(unstructured)
. estimates store independent
. quietly mixed weight c.age##c.age##girl || id: age, covariance(unstructured)
residual(sppower, coordinates(age))
. estimates store power
. quietly mixed weight c.age##c.age##girl || id: age, covariance(unstructured)
residual(spgaussian, coordinates(age) phi0(.4))
. estimates store gaussian
. quietly mixed weight c.age##c.age##girl || id: age, covariance(unstructured)
residual(spspherical, coordinates(age))
. estimates store spherical
To pick the best-fitting residual structure, we compare models using AIC and BIC. We use etable to additionally compare estimates across models.
. etable, estimates(*) column(estimates) mstat(aic) mstat(bic)
| independent power gaussian spherical | ||
| Age (years) 8.696 8.933 8.952 8.999 (0.227) (0.259) (0.262) (0.271) Age (years) # Age (years) -1.972 -2.063 -2.073 -2.091 (0.092) (0.105) (0.105) (0.110) Female Girl -0.171 -0.133 -0.116 -0.108 (0.152) (0.170) (0.173) (0.173) Female # Age (years) Girl -0.993 -1.152 -1.148 -1.222 (0.324) (0.368) (0.373) (0.385) Female # Age (years) # Age (years) Girl 0.289 0.351 0.343 0.377 (0.131) (0.150) (0.151) (0.158) Intercept 3.299 3.220 3.231 3.196 (0.107) (0.120) (0.123) (0.122) var(age) 0.240 0.217 0.186 0.204 (0.069) (0.068) (0.060) (0.065) var(_cons) 0.171 0.074 0.113 0.087 (0.061) (0.045) (0.051) (0.047) cov(age,_cons) 0.112 0.127 0.145 0.133 (0.044) (0.035) (0.032) (0.033) var(e) 0.316 0.424 0.411 0.425 (0.038) (0.055) (0.049) (0.047) rho 0.119 (0.054) phi 0.254 (0.032) Range 0.983 (0.109) AIC 630.89 596.22 607.71 591.59 BIC 666.73 635.64 647.13 631.01 | ||
We see that the spherical structure fits best, followed closely by the power structure. We also see differences in fixed-effects coefficient estimates, especially for the interaction between age and girl.
The rho parameter of the power structure indicates that the residual correlation between any two measurement occasions one year apart is 0.12 (SE = 0.05). Occasions six months apart imply a residual correlation of 0.12^0.5 = 0.35, and so on.
The phi parameter of the Gaussian structure indicates that when d = phi, here at 0.254 (SE = 0.03) years apart, about three months, the correlation between residuals is exp(-12) = 0.368.
The Range parameter of the spherical structure indicates the point at which the spatial error correlation reaches exactly 0, here at about one year (0.98, SE = 0.11).
Using the spherical-structure model, we use estat wcorrelation again to see the ultimate model-implied correlation matrix for the first child.
. estat wcorrelation, list Standard deviations and correlations for id = 45: Standard deviations:
| obs | 1 2 3 4 5 6 | |
| sd | 0.715 0.743 0.880 1.067 1.144 1.472 |
| obs | 1 2 3 4 5 6 | |
| 1 | 1.000 | |
| 2 | 0.831 1.000 | |
| 3 | 0.376 0.505 1.000 | |
| 4 | 0.326 0.380 0.639 1.000 | |
| 5 | 0.339 0.394 0.579 0.888 1.000 | |
| 6 | 0.370 0.430 0.603 0.710 0.744 1.000 |
| id weight age girl | |
| 1. | 45 4.14 0.00 Boy |
| 2. | 45 5.17 0.14 Boy |
| 3. | 45 10.86 0.66 Boy |
| 4. | 45 13.15 1.22 Boy |
| 5. | 45 13.20 1.43 Boy |
| 6. | 45 15.88 2.27 Boy |
This structure is similar to what we saw with the independent-structure model but with larger standard deviations in the first several occasions and higher correlations between closer occasions.
Conclusion
The new spatial residual structures for the mixed command in StataNow offer means to account for correlations across space and time in your model. For this dataset, an unbalanced longitudinal dataset with unequally spaced measurement occasions, we saw that the spherical structure fit best.
These spatial structures offer additional flexibility. You can specify a nugget effect to account for measurement error or microscale variation and choose among Euclidean, Manhattan, and maximum distance metrics. You can also specify initial values for the spatial correlation parameters by using phi0(), range0(), and rho0() or obtain them from the distribution of within-group distances by using distptile0(). Read more about spatial residual-error structures in the Stata Multilevel Mixed-Effects Reference Manual; see [ME] mixed.
References
Cressie, N. Statistics for Spatial Data. 2015. New York: John Wiley & Sons.
Ferron, J., R. Dailey, and Q. Yi. 2002. Effects of misspecifying the first-level error structure in two-level models of change. Multivariate Behavioral Research 37: 379–403. 10.1207/S15327906MBR3703_4.
Goldstein, H. 1986. Efficient statistical modeling of longitudinal data. Human Biology 13: 129–142.
— Meghan Cain
Director of Educational Services
