Stata's mixed command now supports spatial error covariance structures, letting you model correlation between errors as a function of distance computed from coordinate variables. This feature is a part of StataNow™.

Spatial data arise in many scientific and applied settings: soil samples across a floodplain, crop yields across field plots, or sensor measurements across a region. In these settings, observations that are close together tend to have more similar errors than observations that are far apart. Modeling that dependence can reduce bias and produce more accurate inference.

In a multilevel mixed-effects model, random effects capture systematic differences across groups (such as farms, sites, or soil types). Spatial dependence often remains within groups because nearby locations share unobserved conditions. Spatial error covariance structures let you represent that within-group dependence directly through a distance-based correlation function.

You specify a spatial structure inside the residuals() option, along with the coordinate variables used to compute distances via the coordinates() suboption.

Four new spatial structures

mixed adds four widely used spatial correlation kernels:

spexponential		spatial exponential—correlation decays exponentially with distance
spgaussian		spatial Gaussian—correlation decays smoothly and is flat near the origin
spspherical		spatial spherical—correlation becomes exactly 0 beyond a finite range
sppower		spatial power—correlation follows a power law

Each structure can be fit using Euclidean distances by default, or using alternative metrics (Manhattan or maximum) via the metric() suboption.

Real spatial measurements often include microscale variation and measurement error. A nugget effect captures this by allowing correlation at zero distance to be less than one. All four spatial structures support the nugget suboption, and Stata reports the nugget proportion \(eta\), which summarizes how much of the total error variance is attributable to the nugget component.

These spatial error covariance structures are also useful for modeling intensive longitudinal data with repeated measurements taken at irregular time points.

Let's see it work

Modeling spatially correlated errors in a field trial

To illustrate spatial error covariance structures in mixed, we analyze data from a simulated agricultural field trial measuring wheat yield across multiple farms. The fictional study design includes 12 farms, with 60 plots measured per farm, totaling 720 observations. At each plot, wheat yield (yield) was recorded along with nitrogen fertilizer application (nitrogen) and rainfall (rainfall). The spatial coordinates of each plot within its farm are stored in variables x and y, measured in meters.

. webuse wheatyield
(Fictional dataset of wheat yield from multi-farm agricultural trial)

. describe

Contains data from wheatyield.dta
 Observations:           720                  
    Variables:             9                  26 Nov 2025 21:03


Variable      Storage   Display    Value
    name         type    format    label      Variable label

plotid          int     %12.0g                Plot ID
farmid          byte    %12.0g                Farm ID
yield           double  %10.0g                Wheat yield (tons/ha)
nitrogen        double  %10.0g                Nitrogen fertilizer application
                                                (kg/ha)
rainfall        double  %10.0g                Rainfall during growing season
                                                (mm)
nitrogen_c      double  %10.0g                Nitrogen centered at 100 kg/ha
rainfall_c      double  %10.0g                Rainfall centered at 500 mm
x               double  %10.0g                X coordinate within farm (meters)
y               double  %10.0g                Y coordinate within farm (meters)

Sorted by:

We want to quantify the effects of nitrogen and rainfall on yield while accounting for two sources of dependence: between-farm variation, captured by random intercepts, and within-farm spatial correlation, modeled through a spatial error covariance structure. Because yield measurements are subject to measurement error and because microscale spatial variation exists below the sampling resolution, we also allow for a nugget effect.

Spatial spherical covariance structure with a nugget effect

We fit a linear mixed-effects model in which the predictors nitrogen_c and rainfall_c are centered at their sample means. Spatial correlation is modeled using a spherical spatial (spspherical) error covariance structure with a nugget effect (via the nugget suboption). The spatial coordinates x and y define plot locations within each farm on a 300-meter-by-300-meter grid.

. mixed yield nitrogen_c rainfall_c || farmid:, 
     residuals(spspherical, coordinates(x y) nugget) reml

Obtaining starting values by EM ...

Performing gradient-based optimization:
Iteration 0:  Log restricted-likelihood = -466.70434  (not concave)
Iteration 1:  Log restricted-likelihood = -453.72989
Iteration 2:  Log restricted-likelihood = -452.87084
Iteration 3:  Log restricted-likelihood = -452.52054
Iteration 4:  Log restricted-likelihood = -452.51159
Iteration 5:  Log restricted-likelihood = -452.51046
Iteration 6:  Log restricted-likelihood = -452.51046

Computing standard errors ...

Mixed-effects REML regression                        Number of obs    =    720
Group variable: farmid                               Number of groups =     12
                                                     Obs per group:
                                                                  min =     60
                                                                  avg =   60.0
                                                                  max =     60
                                                     Wald chi2(2)     = 892.70
Log restricted-likelihood = -452.51046               Prob > chi2      = 0.0000



       yield 
 
 Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
 
 
 
  nitrogen_c 
 
   .0145382   .0005215    27.88   0.000     .0135161    .0155602
  rainfall_c 
 
   .0082663   .0007632    10.83   0.000     .0067704    .0097622
       _cons 
 
   4.384632   .1695019    25.87   0.000     4.052414    4.716849




  Random-effects parameters  
 
   Estimate   Std. err.     [95% conf. interval]
 
 
 
farmid: Identity             
 
                                                
                  var(_cons) 
 
   .3347489   .1472153      .1413771    .7926093
 
 
 
Residual: Spatial spherical  
 
                                                
                       Range 
 
   109.2307   10.80181      89.98486    132.5928
                      Nugget 
 
   .3861334   .0704936      .2599008    .5297884
                      var(e) 
 
   .2187123   .0138238      .1932292    .2475562

LR test vs. linear model: chi2(3) = 673.00                Prob > chi2 = 0.0000

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

The estimated range parameter indicates that plots within roughly 109 meters exhibit spatial correlation, while plots farther apart are effectively uncorrelated. The nugget accounts for about 39% of the total error variance, suggesting substantial microscale variation and measurement error.

We store this model for further analysis:

. estimates store spher_nug

Assessing stability using alternative initial values

Spatial covariance parameters enter the likelihood nonlinearly, so it is often good practice to assess sensitivity to starting values. For the spherical structure, this can be done by varying the initial value of the range parameter using the range0() suboption within residuals(). Below, we specify four initial values (50, 100, 200, and 400) for the range parameter, and we store the corresponding fitted models as sph_50, sph_100, sph_200, and sph_400.

. local init_vals 50 100 200 400

. foreach i of local init_vals {
     quietly mixed yield nitrogen_c rainfall_c || farmid:,
     residuals(spspherical, coordinates(x y) range0(`i') nugget) reml
     estimates store sph_`i'
  }

We summarize the results from the above four fitted models by using the etable command. We also report the log restricted-likelihood (via the mstat(ll) option).

. etable, estimates(sph_50 sph_100 sph_200 sph_400) 
     column(estimates) mstat(ll)



                                sph_50   sph_100  sph_200  sph_400

Nitrogen centered at 100 kg/ha    0.015    0.015    0.015    0.015
                                (0.001)  (0.001)  (0.001)  (0.001)
Rainfall centered at 500 mm       0.008    0.008    0.008    0.008
                                (0.001)  (0.001)  (0.001)  (0.001)
Intercept                         4.385    4.385    4.385    4.385
                                (0.170)  (0.170)  (0.170)  (0.170)
var(_cons)                        0.335    0.335    0.335    0.335
                                (0.147)  (0.147)  (0.147)  (0.147)
Range                           109.231  109.231  109.231  109.230
                               (10.787) (10.805) (10.802) (10.804)
Nugget                            0.386    0.386    0.386    0.386
                                (0.070)  (0.070)  (0.070)  (0.070)
var(e)                            0.219    0.219    0.219    0.219
                                (0.014)  (0.014)  (0.014)  (0.014)
Log restricted-likelihood       -452.51  -452.51  -452.51  -452.51

All four models yield identical results, indicating that the fitted model is robust to the choice of initial range. In applications where different initial values lead to different solutions, one would typically retain the model with the highest log restricted-likelihood.

Comparing alternative spatial-error structures

To explore alternative specifications, we also fit two additional models: one using the same spherical structure without a nugget effect and one using a Gaussian spatial-error covariance structure.

. quietly mixed yield nitrogen_c rainfall_c || farmid:, residuals(spspherical, coordinates(x y)) reml nolog

. estimates store spher

. mixed yield nitrogen_c rainfall_c || farmid:, 
     residuals(spgaussian, coordinates(x y)) reml nolog

Mixed-effects REML regression                        Number of obs    =    720
Group variable: farmid                               Number of groups =     12
                                                     Obs per group:
                                                                  min =     60
                                                                  avg =   60.0
                                                                  max =     60
                                                     Wald chi2(2)     = 798.80
Log restricted-likelihood = -483.63586               Prob > chi2      = 0.0000



       yield 
 
 Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
 
 
 
  nitrogen_c 
 
   .0134571   .0005136    26.20   0.000     .0124505    .0144636
  rainfall_c 
 
   .0090491   .0007753    11.67   0.000     .0075295    .0105687
       _cons 
 
   4.374288   .1672699    26.15   0.000     4.046445    4.702131




  Random-effects parameters  
 
   Estimate   Std. err.     [95% conf. interval]
 
 
 
farmid: Identity             
 
                                                
                  var(_cons) 
 
   .3296813    .143176      .1407438    .7722525
 
 
 
Residual: Spatial Gaussian   
 
                                                
                         phi 
 
   26.45716   1.511448       23.6546    29.59176
                      var(e) 
 
   .2232835    .012767      .1996118    .2497623

LR test vs. linear model: chi2(2) = 610.75                Prob > chi2 = 0.0000

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

. estimates store gauss

Compared with the mixed model with the spherical structure, the Gaussian structure replaces the finite-range parameter (range) with a scale parameter phi, which controls how quickly spatial correlation decays with distance. In this example, observations become approximately uncorrelated (correlation <= 0.05) at distances greater than about \(3\phi\).

We compare the three candidate models using information criteria:

. estimates stats spher_nug spher gauss

Akaike's information criterion and Bayesian information criterion



       Model 
 
          N   ll(null)  ll(model)      df        AIC        BIC
 
 
 
   spher_nug 
 
        720          .  -452.5105       7   919.0209   951.0757
       spher 
 
        720          .  -465.4181       6   942.8361   970.3117
       gauss 
 
        720          .  -483.6359       6   979.2717   1006.747

Note: BIC uses N = number of observations. See [R] IC note.

The spherical spatial model with a nugget effect provides the best fit for this dataset, as indicated by having the lowest AIC and BIC values.