
ANTE(1)

specifies a firstorder antedependence structure
(Kenward, 1987; Patel, 1991) parameterized in terms of variances and correlation parameters. If t ordered random variables have a firstorder antedependence structure, then each , , is independent of all other , given . This Markovian structure is characterized by its inverse variance matrix, which is tridiagonal. Parameterizing an ANTE(1)
structure for a random vector of size t requires 2t – 1 parameters: variances and t – 1 correlation parameters . The covariances among random variables and are then constructed as
PROC GLIMMIX constrains the correlation parameters to satisfy . For variableorder antedependence models see Macchiavelli and Arnold (1994).

AR(1)

specifies a firstorder autoregressive structure,
The values and are derived for the ith and jth observations, respectively, and are not necessarily the observation numbers. For example, in the following statements the
values correspond to the class levels for the time
effect of the ith and jth observation within a particular subject:
proc glimmix;
class time patient;
model y = x x*x;
random time / sub=patient type=ar(1);
run;
PROC GLIMMIX imposes the constraint for stationarity.

ARH(1)

specifies a heterogeneous firstorder autoregressive structure,
with . This covariance structure has the same correlation pattern as the TYPE=AR(1) structure, but the variances are allowed to
differ.

ARMA(1,1)

specifies the firstorder autoregressive movingaverage
structure,
Here, is the autoregressive parameter, models a movingaverage component, and is a scale parameter. In the notation of Fuller (1976, p. 68), and
The example in Table 43.19 and imply that
where and . PROC GLIMMIX imposes the constraints and for stationarity, although for some values of and in this region the resulting covariance matrix is not positive definite. When the estimated value of becomes negative, the computed covariance is multiplied by to account for the negativity.

CHOL<(q)>

specifies an unstructured variancecovariance matrix
parameterized through its Cholesky root. This parameterization ensures that the resulting variancecovariance matrix is at
least positive semidefinite. If all diagonal values are nonzero, it is positive definite. For example, a unstructured covariance matrix can be written as
Without imposing constraints on the three parameters, there is no guarantee that the estimated variance matrix is positive
definite. Even if and are nonzero, a large value for can lead to a negative eigenvalue of . The Cholesky root of a positive definite matrix is a lower triangular matrix such that . The Cholesky root of the above matrix can be written as
The elements of the unstructured variance matrix are then simply , , and . Similar operations yield the generalization to covariance matrices of higher orders.
For example, the following statements model the covariance matrix of each subject as an unstructured matrix:
proc glimmix;
class sub;
model y = x;
random _residual_ / subject=sub type=un;
run;
The next set of statements accomplishes the same, but the estimated matrix is guaranteed to be nonnegative definite:
proc glimmix;
class sub;
model y = x;
random _residual_ / subject=sub type=chol;
run;
The GLIMMIX procedure constrains the diagonal elements of the Cholesky root to be positive. This guarantees a unique solution
when the matrix is positive definite.
The optional order parameter determines how many bands below the diagonal are modeled. Elements in the lower triangular portion of in bands higher than q are set to zero. If you consider the resulting covariance matrix , then the order parameter has the effect of zeroing all offdiagonal elements that are at least q positions away from the diagonal.
Because of its good computational and statistical properties, the Cholesky root parameterization is generally recommended
over a completely unstructured covariance matrix (TYPE=UN). However, it is computationally slightly more involved.

CS

specifies the compoundsymmetry structure, which
has constant variance and constant covariance
The compound symmetry structure arises naturally with nested random effects, such as when subsampling error is nested within
experimental error. The models constructed with the following two sets of GLIMMIX statements have the same marginal variance
matrix, provided is positive:
proc glimmix;
class block A;
model y = block A;
random block*A / type=vc;
run;
proc glimmix;
class block A;
model y = block A;
random _residual_ / subject=block*A
type=cs;
run;
In the first case, the block*A
random effect models the Gside experimental error. Because the distribution defaults to the normal, the matrix is of form (see Table 43.20), and is the subsampling error variance. The marginal variance for the data from a particular experimental unit is thus . This matrix is of compound symmetric form.
Hierarchical random assignments or selections, such as subsampling or splitplot designs, give rise to compound symmetric
covariance structures. This implies exchangeability of the observations on the subunit, leading to constant correlations between
the observations. Compound symmetric structures are thus usually not appropriate for processes where correlations decline
according to some metric, such as spatial and temporal processes.
Note that Rside compoundsymmetry structures do not impose any constraint on . You can thus use an Rside TYPE=CS structure to emulate a variancecomponent model with unbounded estimate of the variance
component.

CSH

specifies the heterogeneous compoundsymmetry structure, which is an
equicorrelation structure but allows for different variances

FA(q)

specifies the factoranalytic structure with q factors
(Jennrich and Schluchter, 1986). This structure is of the form , where is a rectangular matrix and is a diagonal matrix with t different parameters. When , the elements of in its upperright corner (that is, the elements in the ith row and jth column for ) are set to zero to fix the rotation of the structure.

FA0(q)

specifies a factoranalytic structure with q factors of the
form , where is a rectangular matrix and t is the dimension of . When , is a lower triangular matrix. When —that is, when the number of factors is less than the dimension of the matrix—this structure is nonnegative definite but not
of full rank. In this situation, you can use it to approximate an unstructured covariance matrix.

HF

specifies a covariance structure that satisfies the general
HuynhFeldt condition (Huynh and Feldt, 1970). For a random vector with t elements, this structure has positive parameters and covariances
A covariance matrix generally satisfies the HuynhFeldt condition if it can be written as . The preceding parameterization chooses . Several simpler covariance structures give rise to covariance matrices that also satisfy the HuynhFeldt condition. For
example, TYPE=CS, TYPE=VC, and TYPE=UN(1) are nested within TYPE=HF. You can use the COVTEST statement to test the HF structure against one of these simpler structures. Note also that the HF structure is nested within
an unstructured covariance matrix.
The TYPE=HF covariance structure can be sensitive to the choice of starting values and the default MIVQUE(0) starting values
can be poor for this structure; you can supply your own starting values with the PARMS statement.

LIN(q)

specifies a general linear covariance structure with q parameters.
This structure consists of a linear combination of known matrices that you input with the LDATA= option. Suppose that you want to model the covariance of a random vector of length t, and further suppose that are symmetric ) matrices constructed from the information in the LDATA= data set. Then,
where denotes the element in row i, column j of matrix .
Linear structures are very flexible and general. You need to exercise caution to ensure that the variance matrix is positive
definite. Note that PROC GLIMMIX does not impose boundary constraints on the parameters of a general linear covariance structure. For example, if classification variable A
has 6 levels, the following statements fit a variance component structure for the random effect without boundary constraints:
data ldata;
retain parm 1 value 1;
do row=1 to 6; col=row; output; end;
run;
proc glimmix data=MyData;
class A B;
model Y = B;
random A / type=lin(1) ldata=ldata;
run;

PSPLINE<(options)>

requests that PROC GLIMMIX form a Bspline basis and fits a
penalized Bspline (Pspline, Eilers and Marx 1996) with random spline coefficients. This covariance structure is available only for Gside random effects and only a single
continuous random effect can be specified with TYPE=PSPLINE. As for TYPE=RSMOOTH, PROC GLIMMIX forms a modified matrix and fits a mixed model in which the random variables associated with the columns of are independent with a common variance. The matrix is constructed as follows.
Denote as the matrix of Bsplines of degree d and denote as the matrix of rthorder differences. For example, for K = 5,
Then, the matrix used in fitting the mixed model is the matrix
The construction of the Bspline knots is controlled with the KNOTMETHOD= EQUAL(m) option and the DEGREE=d suboption of TYPE=PSPLINE. The total number of knots equals the number m of equally spaced interior knots plus d knots at the low end and knots at the high end. The number of columns in the Bspline basis equals K = m + d + 1. By default, the interior knots exclude the minimum and maximum of the randomeffect values and are based on m – 1 equally spaced intervals. Suppose and are the smallest and largest randomeffect values; then interior knots are placed at
In addition, d evenly spaced exterior knots are placed below and exterior knots are placed above . The exterior knots are evenly spaced and start at times the machine epsilon. For example, based on the defaults d = 3, r = 3, the following statements lead to 26 total knots and 21 columns in , m = 20, K = m + d + 1 = 24, K – r = 21:
proc glimmix;
model y = x;
random x / type=pspline knotmethod=equal(20);
run;
Details about the computation and properties of Bsplines can be found in de Boor (2001). You can extend or limit the range of the knots with the KNOTMIN= and KNOTMAX= options. Table 43.18 lists some of the parameters that control this covariance type and their relationships.
Table 43.18: PSpline Parameters
Parameter

Description

d

Degree of Bspline, default d = 3

r

Order of differencing in construction of , default r = 3

m

Number of interior knots, default


Total number of knots


Number of columns in Bspline basis


Number of columns in

You can specify the following options for TYPE=PSPLINE:
 DEGREE=d

specifies the degree of the Bspline. The default is d = 3.
 DIFFORDER=r

specifies the order of the differencing matrix . The default and maximum is r = 3.

RSMOOTH<(m  NOLOG)>

specifies a radial smoother covariance structure for Gside random
effects. This results in an approximate lowrank thinplate spline where the smoothing parameter is obtained by the estimation
method selected with the METHOD= option of the PROC GLIMMIX statement. The smoother is based on the automatic smoother in Ruppert, Wand, and Carroll (2003, Chapter 13.4–13.5), but with a different method of selecting the spline knots. See the section Radial Smoothing Based on Mixed Models for further details about the construction of the smoother and the knot selection.
Radial smoothing is possible in one or more dimensions. A univariate smoother is obtained with a single random effect, while
multiple random effects in a RANDOM statement yield a multivariate smoother. Only continuous random effects are permitted
with this covariance structure. If denotes the number of continuous random effects in the RANDOM statement, then the covariance structure of the random effects
is determined as follows. Suppose that denotes the vector of random effects for the ith observation. Let denote the vector of knot coordinates, , and K is the total number of knots. The Euclidean distance between the knots is computed as
and the distance between knots and effects is computed as
The matrix for the GLMM is constructed as
where the matrix has typical element
and the matrix has typical element
The exponent in these expressions equals , where the optional value m corresponds to the derivative penalized in the thinplate spline. A larger value of m will yield a smoother fit. The GLIMMIX procedure requires p > 0 and chooses by default m = 2 if and otherwise. The NOLOG option removes the and terms from the computation of the and matrices when is even; this yields invariance under rescaling of the coordinates.
Finally, the components of are assumed to have equal variance . The “smoothing parameter” of the lowrank spline is related to the variance components in the model, . See Ruppert, Wand, and Carroll (2003) for details. If the conditional distribution does not provide a scale parameter , you can add a single Rside residual parameter.
The knot selection is controlled with the KNOTMETHOD= option. The GLIMMIX procedure selects knots automatically based on the vertices of a kd tree or reads knots from a data set that you supply. See the section Radial Smoothing Based on Mixed Models for further details on radial smoothing in the GLIMMIX procedure and its connection to a mixed model formulation.

SIMPLE

is an alias for TYPE=VC.

SP(EXP)(clist)

models an exponential spatial or temporal covariance structure,
where the covariance between two observations depends on a distance metric . The clist contains the names of the numeric variables used as coordinates to determine distance. For a stochastic process in , there are k elements in clist. If the vectors of coordinates for observations i and j are and , then PROC GLIMMIX computes the Euclidean distance
The covariance between two observations is then
The parameter is not what is commonly referred to as the range parameter in geostatistical applications. The practical range of a (secondorder
stationary) spatial process is the distance at which the correlations fall below 0.05. For the SP(EXP) structure, this distance is . PROC GLIMMIX constrains to be positive.

SP(GAU)(clist)

models a Gaussian covariance structure,
See TYPE=SP(EXP) for the computation of the distance . The parameter is related to the range of the process as follows. If the practical range is defined as the distance at which the correlations fall below 0.05, then . PROC GLIMMIX constrains to be positive. See TYPE=SP(EXP) for the computation of the distance from the variables specified in clist.

SP(MAT)(clist)

models a covariance structure in the Matérn class of covariance
functions (Matérn, 1986). The covariance is expressed in the parameterization of Handcock and Stein (1993); Handcock and Wallis (1994); it can be written as
The function is the modified Bessel function of the second kind of (real) order . The smoothness (continuity) of a stochastic process with covariance function in the Matérn class increases with . This class thus enables datadriven estimation of the smoothness properties of the process. The covariance is identical
to the exponential model for (TYPE=SP(EXP)(clist)), while for the model advocated by Whittle (1954) results. As , the model approaches the Gaussian covariance structure (TYPE=SP(GAU)(clist)).
Note that the MIXED procedure offers covariance structures in the Matérn class in two parameterizations, TYPE=SP(MATERN) and
TYPE=SP(MATHSW). The TYPE=SP(MAT) in the GLIMMIX procedure is equivalent to TYPE=SP(MATHSW) in the MIXED procedure.
Computation of the function and its derivatives is numerically demanding; fitting models with Matérn covariance structures can be timeconsuming. Good
starting values are essential.

SP(POW)(clist)

models a power covariance structure,
where . This is a reparameterization of the exponential structure, TYPE=SP(EXP). Specifically, . See TYPE=SP(EXP) for the computation of the distance from the variables specified in clist. When the estimated value of becomes negative, the computed covariance is multiplied by to account for the negativity.

SP(POWA)(clist)

models an anisotropic power covariance structure in k dimensions,
provided that the coordinate list clist has k elements. If denotes the coordinate for the ith observation of the mth variable in clist, the covariance between two observations is given by
Note that for k = 1, TYPE=SP(POWA) is equivalent to TYPE=SP(POW), which is itself a reparameterization of TYPE=SP(EXP). When the estimated
value of becomes negative, the computed covariance is multiplied by to account for the negativity.

SP(SPH)(clist)

models a spherical covariance structure,
The spherical covariance structure has a true range parameter. The covariances between observations are exactly zero when
their distance exceeds . See TYPE=SP(EXP) for the computation of the distance from the variables specified in clist.

TOEP

models a Toeplitz covariance structure. This structure can be viewed
as an autoregressive structure with order equal to the dimension of the matrix,

TOEP(q)

specifies a banded Toeplitz structure,
This can be viewed as a movingaverage structure with order equal to q – 1. The specification TYPE=TOEP(1) is the same as , and it can be useful for specifying the same variance component for several effects.

TOEPH<(q)>

models a Toeplitz covariance structure. The correlations of this
structure are banded as the TOEP or TOEP(q) structures, but the variances are allowed to vary:
The correlation parameters satisfy . If you specify the optional value q, the correlation parameters with are set to zero, creating a banded correlation structure. The specification TYPE=TOEPH(1) results in a diagonal covariance
matrix with heterogeneous variances.

UN<(q)>

specifies a completely general (unstructured) covariance matrix
parameterized directly in terms of variances and covariances,
The variances are constrained to be nonnegative, and the covariances are unconstrained. This structure is not constrained
to be nonnegative definite in order to avoid nonlinear constraints; however, you can use the TYPE=CHOL structure if you want
this constraint to be imposed by a Cholesky factorization. If you specify the order parameter q, then PROC GLIMMIX estimates only the first q bands of the matrix, setting elements in all higher bands equal to 0.

UNR<(q)>

specifies a completely general (unstructured) covariance matrix
parameterized in terms of variances and correlations,
where denotes the standard deviation and the correlation is zero when and when , provided the order parameter q is given. This structure fits the same model as the TYPE=UN(q) option, but with a different parameterization. The ith variance parameter is . The parameter is the correlation between the ith and jth measurements; it satisfies . If you specify the order parameter q, then PROC GLIMMIX estimates only the first q bands of the matrix, setting all higher bands equal to zero.

VC

specifies standard variance components and is the default structure
for both Gside and Rside covariance structures. In a Gside covariance structure, a distinct variance component is assigned
to each effect. In an Rside structure TYPE=VC is usually used only to add overdispersion effects or with the GROUP= option to specify a heterogeneous variance model.