Defines a variance structure for random effects in a `REML`

model.

### Options

`TERMS` = formula |
Model terms for which the covariance structure is to be defined |
---|---|

`FORMATION` = string token |
Whether the structure is formed by direct product construction or by definition of the whole matrix (`direct` , `whole` ); default `dire` |

`CORRELATE` = string token |
Whether to impose correlation across the model terms if several are specified (`none` , `positivedefinite` , `unrestricted` ); default `none` |

`CINITIAL` = scalars |
Initial values for covariance matrix across terms |

`COORDINATES` = matrix or variates |
Coordinates of the data points to be used in calculating distance-based models |

### Parameters

`MODELTYPE` = string tokens |
Type of covariance model associated with the term(s), or with individual factors in the term(s) if `FORMATION=direct` (`identity` , `fixed` , `AR` , `MA` , `ARMA` , `power` , `boundedlinear` , `circular` , `spherical` , `linearvariance` , `banded` , `correlation` , `antedependence` , `unstructured` , `diagonal` , `uniform` , `FA` , `FAequal` ) default `iden` |
---|---|

`ORDER` = scalar |
Order of model |

`HETEROGENEITY` = string token |
Heterogeneity for correlation matrices (`none` , `outside` ); default `none` |

`METRIC` = string token |
How to calculate distances when `MODELTYPE=power` (`cityblock` , `squared` , `euclidean` ); default `city` |

`FACTOR` = factors |
Factors over which to form direct products |

`MATRIX` = symmetric matrices, diagonal matrices or pointers |
Defines matrix values for a term or the factors when `MODELTYPE=fixed` |

`INVERSE` = symmetric matrices, diagonal matrices or pointers |
Define values for matrix inverses (instead of the fixed matrices themselves) when `MODELTYPE=fixed` |

`DISTANCES` = symmetric matrices |
Symmetric matrix of pre-formed distances to be used in distance-based models of order one |

`COORDINATES` = matrices, variates or pointers |
Specifies coordinates of each factor level to be used in calculating distance-based models |

`INITIAL` = scalars, variates, matrices, symmetric matrices or pointers |
Initial parameter values for each correlation matrix (supplied in the structures appropriate for the model concerned) |

`CONSTRAINTS` = texts |
Texts containing strings `none` , `fix` or `positive` to define constraints for the parameters in each model |

`EQUALITYCONSTRAINTS` = variates |
Non-zero values in the variate indicate groups of parameters whose values are to be constrained to be equal |

### Description

`VSTRUCTURE`

can be used to define the form of covariance structure for any term in the random model defined for `REML`

by `VCOMPONENTS`

. By default, the effects for each random term are assumed to be independent with common variance σ_{j}^{2} for term *j*, that is, the random term has covariance matrix σ_{j}^{2}*I*. `VSTRUCTURE`

is used to define correlation between random effects within terms, to allow a changing variance within a term, and to define correlations between different random terms. These models are particularly useful when fitting linear models to repeated measurements or spatial data and for random coefficient regression.

`VSTRUCTURE`

can only be used after `VCOMPONENTS`

has been used to define the fixed and random models. It can be used more than once to define different structures for different random terms. The information is accumulated within Genstat, and it will all be used by subsequent `REML`

commands. You can check on the model and covariance structures defined at any time by using the `VSTATUS`

directive. To cancel a covariance structure for a term you simply need to use `VSTRUCTURE`

to change the model back to the scaled identity matrix σ_{j}^{2}*I*. To cancel all covariance structures you can give a new `VCOMPONENTS`

command.

For a random term constructed from more than one factor, the covariance matrix can be formed either as a single matrix for the whole term, or as the direct product of several matrices corresponding to the factors. Consider an analysis of repeated measurements where data has been taken weekly from each subject, and one of several different treatments has been applied to each subject. It is likely that data taken from the same subject will be correlated, with correlation decreasing over time, but that subjects will be independent. This corresponds to an *I* ⊗ *C* covariance structure, where the identity matrix *I* corresponds to the independent subjects, and the covariance matrix *C* corresponds to the correlated measurements over time within subjects. If we take *C* to be an auto-regressive process of order 1, this can be defined and fitted as follows:

`VCOMPONENTS [FIXED=Tmt] RANDOM=Subject.Week`

`VSTRUCTURE [TERM=Subject.Week] MODELTYPE=I,AR; ORDER=1;\`

` FACTOR=Subject,Week`

`REML Y`

The `TERM`

option is used to specify the term to which the covariance structure is to be applied. For each factor in the term you can then specify the covariance model to be applied (see below for list of available models). However, it is not necessary to specify factors for which the default identity model is required, so the following is an equivalent specification:

`VCOMPONENTS [FIXED=Tmt] RANDOM=Subject.Week`

`VSTRUCTURE [TERM=Subject.Week] MODELTYPE=AR; ORDER=1; FACTOR=Week`

To cancel the covariance structure for the term, a null setting is sufficient:

`VSTRUCTURE [TERM=Subject.Week]`

Although the covariance structure for each term here is of the form *G _{j}* =

*I*, the variance matrix for the data is of the form

*V* = σ^{2 }( Σ_{j} γ_{j}*Z _{j}G_{j}Z_{j}*′ +

*I*)

In this case the random subject term generates correlations that are equal across all the times within subjects. It is important to remember that including a random term in the model will generate uniform correlations between units with the same values of the random factor(s). Retaining these terms in the model as well as specifying a correlated structure may be appropriate for some data sets, but can sometimes lead to difficulties in parameter estimation.

The possible settings for the `MODELTYPE`

parameter, generating symmetric covariance matrices *C* (*C _{i}*

_{, j}=

*C*

_{j}_{, i}for all

*i*,

*j*), are listed below. Where more than one model order can be used, the default is shown in bold and can be changed by using the

`ORDER`

option. For the `AR`

, `MA`

, `ARMA`

, `power`

and `banded`

models, the order is the same as the number of parameters to be fitted. For the `banded`

, `correlation`

, `ante-dependence`

and `unstructured`

models, the order is the number of non-zero off-diagonal bands in the matrix. For the FAequal and FA models, the order is the number of columns in the matrix Λ.`identity` |
identity matrix | C_{i}_{, i} = 1, C_{i}_{,j} = 0, for i≠j |

`fixed` |
fixed matrix | C_{i}_{, j} specified |

`AR` |
auto-regressive order |
_{, i} = 1C_{i}_{+1, i} = φ_{1} / (1-φ_{2})C_{i}_{, j} = φ_{1} C–_{i}_{1, j} + φ_{2} C–_{i}_{2, j},i > j+1, -1 < φ_{1}, φ_{2} < 1,|φ _{1}+φ_{2}|<1, φ_{2}-φ_{1}<1, φ_{2}>-1 |

`MA` |
moving average order |
_{, i} = 1C_{i}_{+1, i} = -θ_{1}(1-θ_{2})/(1+θ_{1}^{2}+θ_{2}^{2})C_{i+2, i} = -θ_{2} / (1+θ_{1}^{2}+θ_{2}^{2})C_{i}_{, j} = 0, i>j+2-1 < θ _{1}, θ_{2} < 1, θ_{2}+θ_{1} < 1 |

`ARMA` |
auto-regressive moving-average order 1 |
_{, i} = 1C_{i}_{+1, i} = (θ-φ)(1-φθ)/(1+θ^{2}-2φθ)C_{i}_{, j} = φC–_{i}_{1, j }, i>j+1-1 < φ, θ < 1 |

`power` |
based on distance order |
_{, i} = 1C_{i}_{, j} = φ_{1}^{d1}φ_{2}^{d2
}d_{1}, d_{2} = distance in 1st and 2nd dimensions0 < φ _{1}, φ_{2} < 1 |

`boundedlinear` |
based on distance order 1 |
_{, j} = 1 – d/φ for d ≤ φ,C_{i}_{, j} = 0 for d > φ0 < φ |

`circular` |
based on distance order 1 |
_{, j} = 1 – (2/π) {(d/φ)√(1-(d/φ)^{2}) + sin^{-1}(d/φ)} for d ≤ φ,C_{i}_{, j} = 0 for d > φ0 < φ |

`spherical` |
based on distance order 1 |
_{, j} = 1 – 1.5 (d/φ) + 0.5 (d/φ)^{3} for d ≤ φ,C_{i}_{, j} = 0 for d > φ0 < φ |

`linearvariance` |
based on distance order 1 |
_{, j} = 1 – 2φ d / max(d)0 < φ < 1 |

`banded` |
equal bands |
_{, i} = 1C_{i}_{+k, i} = θ_{k }, 1 < k < order-1 < θ _{k} < 1C_{i}_{+k, i} = 0, otherwise |

`correlation` |
general correlation matrix |
_{, i} = 1C_{i}_{, j} = θ_{ij },1 < | i–j| ≤ orderC_{i}_{, j} = 0, |i–j| > order-1 < θ _{ij} < 1 |

`uniform` |
uniform matrix | C_{i}_{, j} = θ for all i,j |

`diagonal` |
diagonal matrix |
_{, i} = θ_{i
}C_{i}_{, j} = 0, i≠j |

`antedependence` |
ante-dependence model |
_{, i}^{-1} = d_{i}^{-1},D_{i}_{, j} = 0 for i≠jU_{i}_{, i} = 1,U_{i}_{, j} = u ,_{ij}1 ≤ j–i ≤ orderU_{i}_{, j} = 0, for i>j |

`unstructured` |
general covariance matrix |
_{, j} = θ_{ij },0 < | i–j| ≤ orderC_{i}_{, j} = 0, |i–j| > order |

`FA` |
factor analytic |
= ψ_{i} for _{i}i=1…nrows |

`FAequal` |
factor analytic with common variance |
= ψ for _{i}i=1…nrows |

Initial parameter values can be specified using the `INITIAL`

parameter. For most models, the number of initial values required is the number of parameters, and default values will be generated. However, for `unstructured`

models, a full covariance matrix of initial values must be given, and for the `correlation`

model a full correlation matrix must be provided. For the `ante-dependence`

model, either a full covariance matrix can be provided, or a pointer to a *U* and a *D*^{-1} matrix of the correct forms. For the `FA`

and `FAequal`

models, a pointer must be used to give the initial Λ and Ψ matrices, otherwise default initial values are generated. The `FAequal`

model can be used to get initial values for the `FA`

model. Initial values are required for these models because the algorithm may not converge when many parameters are fitted if the starting values are not realistic. Initial values might be generated from covariance matrices estimated by fitting simpler models, or from residuals from a null variance model.

A missing value in the initial values is taken to mean that the value is inestimable and it will be fixed at a small value for the analysis. Alternatively, a parameter can be fixed at its initial value using the `CONSTRAINTS`

parameter. The codes (not case sensitive and able to be abbreviated) may take value `fix`

to indicate the parameter is to be fixed at its initial value, `positive`

to indicate it is to remain positive or `none`

to indicate no constraints. The default is a positive constraint or no constraint depending on context; for example scaling parameters are always constrained to remain positive. The `EQUALITYCONSTRAINTS`

parameter allows you to constrain some of the parameters to have the same value. The variate that it specifies contains a zero value if there is no constraint, and an identical integer value for any set of parameters whose values are to be equal. So, a variate containing the values (0,1,2,1,2) would constrain the second parameter to be equal to the fourth parameter, and the third parameter to be equal to the fifth parameter.

It may sometimes be desirable to allow for unequal variances for the models defined in terms of correlation matrices: that is, for the `AR`

, `MA`

, `ARMA`

, `uniform`

, `power`

, `boundedlinear`

, `circular`

, `spherical`

, `linearvariance`

, `banded`

and `correlation`

models. This can be done by setting option `HETEROGENEITY=outside`

. This means a diagonal matrix *D *of standard errors will be applied to the correlation matrix *C* to generate a matrix *D*^{½}*CD*^{½}. In this case, a number of extra parameters (equal to the number of effects in the factor or term) should be added to the vector of initial values. These models allow investigation of a structured correlation pattern for changing variances and are particularly useful in the analysis of repeated measurements data when variance increases over time. For example, to allow for changing variance over time in our example above, we can specify

`VCOMPONENTS [FIXED=Tmt] RANDOM=Subject.Week`

`VSTRUCTURE [TERM=Subject.Week] MODELTYPE=AR; ORDER=1;\`

` FACTOR=Week; HETEROGENEITY=outside`

`REML Y`

In some circumstances, you may wish to define a single model to apply to the whole term, instead of using the direct product form illustrated above. In this case, you should set option `FORM=whole`

. Note that, when a term consists of a single factor, it is not necessary to set the `FACTOR`

option.

With `MODELTYPE=fixed`

, you must either use the `MATRIX`

option to specify the values of the covariance matrix *C*, or the `INVERSE`

option to specify the inverse matrix. The values of the matrix or its inverse can be supplied as diagonal matrices or symmetric matrices. In addition, values for the inverse matrix can be supplied in sparse form as a pointer. The output from `VPEDIGREE`

is designed for input here, but you can also define the inverse matrix explicitly. The second element of the pointer should then be a variate containing the non-zero values of the inverse in lower triangular order. The first element should be a factor, with number of levels equal to the number of rows *n*(*n*+1)/2 of the matrix. This has firstly a block of *n* values giving the position in the variate of the first value stored for each row. There is then a block of values for each row in turn, giving the columns in which each non-zero value appears.

When `MODELTYPE=power`

is used to define a distance-based model, the model can be of order 1 (isotropic) or 2 (anisotropic). For models with `ORDER=1`

, a single set of distances must be formed. The necessary information can be supplied using either the `COORDINATES`

option, or the `COORDINATES`

parameter, or the `DISTANCES`

parameter. With the `COORDINATES`

option you can specify either a matrix, or a list of variates, to define multi-dimensional coordinates for each unit of the data. The length of the variates, or the number of rows of the matrix, must be equal to the number of data values. The number of variates, or the number of columns of the matrix, is equal to the number of dimensions. The coordinates for the levels of each `FACTOR`

are then calculated as the mean values of the coordinates of the units included in the analysis with those levels. Alternatively, you can use the `COORDINATES`

parameter to specify a single variate, a pointer to several variates or a matrix to define multi-dimensional coordinates for each level of the `FACTOR`

. This parameter takes precedence over the `COORDINATES`

option. The length of the variates, or the number of rows of the matrix, must be equal to the number of levels of the `FACTOR`

. The number of variates, or the number of columns of the matrix, is again equal to the number of dimensions.

The distance calculation is defined by the `METRIC`

option. For levels *i* and *j* with *n*-dimensional coordinates {*c _{ik}*:

*k*=1…

*n*} and {

*c*:

_{jk}*k*=1…

*n*} the distance

*d*is defined as

_{ij} d = Σ_{ij}_{k} |c – _{ik}c|_{jk} |
for `METRIC=cityblock` (the default); |
---|---|

d = Σ_{ij}_{k} (c – _{ik}c)_{jk}^{2} |
for `METRIC=squared` ; and |

d = {Σ_{ij}_{k} (c – _{ik}c)_{jk}^{2}}^{1/2} |
for `METRIC=euclidean` . |

Finally, you can supply a symmetric matrix of pre-calculated distances, using the `DISTANCES`

parameter, and this takes precedence over the `COORDINATES`

parameter and option. The number of rows of the `DISTANCES`

matrix must be equal to the number of levels of the `FACTOR`

.

When `MODELTYPE=power`

and `ORDER=2`

, the `DISTANCES`

parameter cannot be used, and only two-dimensional coordinates are allowed. The coordinates must be specified using either the `COORDINATES`

option or parameter, as described above. The distances are calculated within each dimension separately, according to the setting of the `METRIC`

option. In this case the Euclidean and city-block distances are equivalent.

The spherical family of geostatistical models correspond to the `MODELTYPE`

settings `boundedlinear`

(for one-dimensional distances), `circular`

(for one or two dimensions) and `spherical`

(for one or two dimensions). For further details, see Webster & Oliver (2007). These models are based on distances, and require coordinates to be supplied using either the `COORDINATES`

option (to give coordinates for each data value), or the `COORDINATES`

parameter (to give coordinates for each factor level), as described for `MODELTYPE=power`

above. The parameter φ is interpreted as the range at which the correlation is considered to have decayed to zero. A small value therefore indicates weak correlation, and a large value indicates stronger correlation. These models do not have continuous second derivatives, and their log-likelihood may be multi-modal. To detect this potential problem, it is therefore important to start their estimation from several different initial values; this can be done using the `INITIAL`

parameter as described above. To ensure that the estimated correlation matrix differs from the identity matrix, it is necessary for the range parameter to be larger than the minimum distance specified by the coordinates; any initial value smaller than this will be adjusted.

The setting `MODELTYPE=linearvariance`

specifies the linear variance model of Williams (1986), extended by Piepho & Williams (2010). This model is parameterized so that the parameter φ lies in the range [0,1], which allows correlations in the range [-1,1]. Values of φ close to one indicate weak correlation and values close to zero indicate strong correlation between neighbouring observations.

The `CORRELATE`

option allows you to specify correlations between model terms which have equal numbers of effects. A common correlation will then be fitted between parallel effects. For example, consider a random coefficient regression model where the fixed model contains common response to covariate `X`

and the random model allows for deviations in the intercept and slope about this line for each subject. The random intercept and slope for each subject may be correlated, but subjects are independent. This correlation across terms is defined using the `CORRELATE`

option as follows:

`VCOMPONENTS [FIXED=X] RANDOM=SUBJECT+SUBJECT.X`

`VSTRUCTURE [SUBJECT+SUBJECT.X; CORRELATE=positivedefinite;\`

` CINITIAL=!(1,0.1,0.3); FORM=whole]`

The `CORRELATE`

option setting `positivedefinite`

is used to ensure that the correlation matrix between the terms remains positive definite. This constraint can be relaxed using the setting `unrestricted`

(an unstructured covariance matrix is then used to describe covariance across the terms). The model fitting is done here in terms of a covariance matrix, where the diagonal elements are the gammas for the correlated terms. The `CINITIAL`

option is used to give initial values for this matrix. If no initial values are given, the initial values are taken from initial gamma values given in `VCOMPONENTS`

when the model is declared. When correlations are declared between terms, you must set `FORMATION=whole`

. In the random coefficient regression model above, no correlation structure is declared within terms since the subjects are independent. However, it is possible to declare correlation/covariance models within terms as usual. For example, an animal breeding model might use `VPEDIGREE`

to set up covariances within terms as follows:

`VPEDIGREE animal; FEMALE=dam; MALE=sire; INVERSE=Ainv`

`VCOMPONENTS [FIXED=Trt] RANDOM=animal+dam+env`

`VSTRUCTURE [animal+dam; CORRELATE=unr; FORM=whole]\`

` MODELTYPE=fixed; INVERSE=Ainv`

These declarations set up random terms with covariance structures of the form:

cov(`animal`

) = σ_{a}^{2} *A*, cov(`dam`

) = σ_{d}^{2} *A*, cov(`animal`

, `dam`

) = σ_{ad} *A*.

**Direct Products**

Although the direct product construction used to build the covariance structures does not generally constrain the models that can be fitted to any data set, you should be aware of the implications that arise when defining covariance structures for the residual term. The REML algorithm used by Genstat detects the presence of the residual term in the model by searching for terms with number of levels equal to the number of data values, *n*. When no covariance structures are specified, the first term with number of levels > *n* is used as the residual. However, when covariance structures are defined, the form of the variance model is

*V* = σ^{2 }( Σ_{j} γ_{j}*Z _{j}G_{j}Z_{j}*′ +

*R*)

where matrix *R* corresponds to the residual term and has *n* rows. For this reason, any term found with > *n* rows will not be used as the residual if it has a covariance matrix. If no valid residual term is found, a residual term will automatically be added to the model. This may result in an extra error term being fitted unintentionally. An example where this may happen is in repeated measurements data where unequal numbers of measurements have been taken on subjects. If direct product construction is used, the matrix generated will have more rows than the data and cannot be used as *R*. A work-around is to put missing values in the data set to give equal replication and use `REML`

option `MVINCLUDE=yvariate`

to retain the missing values in the analysis. Alternatively, you could fix the residual component at a small value.

Note that in the repeated measurements example above, if measurements are taken at different times for each subject, the direct product structure is not appropriate. In this case, a power model may be fitted over the whole term, constraining the between subject correlation to zero:

`VSTRUCTURE [TERM=Subject.Week; FORM=whole;\`

` COORD=subject,week] MODELTYPE=power; ORDER=2;\`

` INITIAL=!(0,0.1); CONSTRAIN=!T(Fix,None)`

Note that the parameters run in the order of the coordinates vectors (which are variate forms of the model factors).

Options: `TERMS`

, `FORMATION`

, `CORRELATE`

, `CINITIAL`

, `COORDINATES`

.

Parameters: `MODELTYPE`

, `ORDER`

, `HETEROGENEITY`

, `METRIC`

, `FACTOR`

, `MATRIX`

, `INVERSE`

, `DISTANCES`

, `COORDINATES`

, `INITIAL`

, `CONSTRAINTS`

, `EQUALITYCONSTRAINTS`

.

### References

Piepho, H.P. & Williams, E.R. (2010). Linear variance models for plant breeding trials. *Plant Breeding*, 129, 1-8.

Webster, R. & Oliver, M.A. (2007). *Geostatistics for Environmental Scientists, 2nd edition*. Wiley, Chichester.

Williams, E.R. (1986). A neighbour model for field experiments. *Biometrika*, 73, 279-87.

### See also

Directives: `REML`

, `VCOMPONENTS`

, `VRESIDUAL`

, `VSTATUS`

.

Procedure: `VFSTRUCTURE`

, `VNEARESTNEIGHBOUR`

.

Commands for: REML analysis of linear mixed models, Repeated measurements.

### Example

" Examples 2:5.4.3a-e, 2:5.4.5a-b, 2:5.7a-c " " Repeated measurements analysis: growth of 14 plants measured after 1,3,5,7,10 weeks." FACTOR Plant & [LABELS=!T(HC,MAV)] Treatment OPEN '%GENDIR%/Examples/GuidePart2/Plant.dat'; CHANNEL=2 READ [CHANNEL=2] Plant,Treatment,Time,Height;\ FREPRESENTATION=levels,labels,levels,levels CLOSE 2 GROUPS Time; FACTOR=Week DREPMEASURES [GROUPS=Plant; TIMEPOINTS=Week] Height " Anova split-plot analysis." VCOMPONENTS [FIXED=Treatment*Week] Plant/Week REML Height " Equivalent analysis using uniform correlation structure" VCOMPONENTS [FIXED=Treatment*Week] Plant.Week VSTRUCTURE [TERM=Plant.Week] MODELTYPE=uniform; FACTOR=Week REML [PRINT=model,components] Height " Power model: - correlation decreases as time between measurements increase - takes account of unequally spaced measurements - co-ordinates must be specified as a list of variates or a matrix." VCOMPONENTS [FIXED=Treatment*Week] Plant.Week VSTRUCTURE [TERM=Plant.Week; COORDINATES=Time] MODELTYPE=power; FACTOR=Week REML [PRINT=model,components,wald,deviance] Height " Heterogeneous power model - correlations follow power model, variance allowed to change over time." VSTRUCTURE [TERM=Plant.Week; COORDINATES=Time] MODELTYPE=power;\ ORDER=1; FACTOR=Week; HETEROGENEITY=outside REML [PRINT=model,components,wald,deviance] Height " Save fitted covariance model as initial values for unstructured model." VKEEP TERM=Plant.Week; COVARIANCEMODEL=Covpower PRINT Covpower['Week'] " Unstructured model." VSTRUCTURE [TERM=Plant.Week] MODELTYPE=unstructured; FACTOR=Week;\ INITIAL=Covpower['Week'] REML [PRINT=model,components,wald,deviance] Height " Alternatively, generate initial values using residuals generated after fitting with no variance model." VCOMPONENTS [FIXED=Treatment*Week] REML [PRINT=*] Height; RESIDUALS=r " Residuals are in order plants within weeks, so matrix rows correspond to weeks and columns to plants." MATRIX [ROWS=5; COLUMNS=14; VALUES=#r] mres SYMMETRIC [ROWS=5] vcov CALCULATE vcov = mres *+ TRANSPOSE(mres) " Dividing by number of replicates within each week gives an easy but conservative estimate as it takes no account of treatment d.f." CALCULATE vcov = vcov / 14 " Unstructured model from new initial values." VCOMPONENTS [FIXED=Treatment*Week] Plant.Week VSTRUCTURE [TERM=Plant.Week] MODELTYPE=unstructured; FACTOR=Week;\ INITIAL=!(#vcov) REML [PRINT=deviance] Height " Ante-dependence model order 1 - also requires initial values." VSTRUCTURE [TERM=Plant.Week] antedependence; FACTOR=Week;\ INITIAL=Covpower['Week']; ORDER=1 REML [PRINT=model,components,deviance] Height " Ante-dependence model order 2 - use initial values from power model again ." VSTRUCTURE [TERM=Plant.Week] antedependence; FACTOR=Week; INITIAL=Covpower['Week']; ORDER=2 REML [PRINT=deviance] Height " Random coefficient regression Growth of 14 plants measured after 1,3,5,7,10 weeks. Profiles suggest use of quadratic functions: fit random intercept and slope for plants - no correlation." CALCULATE Timesqrd = Time * Time VCOMPONENTS [FIXED=Treatment*(Time+Timesqrd)] Plant+Plant.Time REML [PRINT=model,components,deviance] Height " Fit random intercept and slope (with correlation) for plants, using previous estimates as initial values." VCOMPONENTS [FIXED=Treatment*(Time+Timesqrd)] Plant+Plant.Time VSTRUCTURE [TERMS=Plant+Plant.Time; FORMATION=whole;\ CORRELATE=unrestricted; CINITIAL=!(3,0.1,0.1)] REML [PRINT=#,deviance] Height " Save random coefficient regression matrix for use as initial values " VKEEP TERM=Plant; COVARIANCEMODEL=CovRCR PRINT CovRCR['Across terms'] " Random cubic spline models: growth of 14 plants measured after 1,3,5,7,10 weeks. Baseline model without spline terms." VCOMPONENTS [FIXED=Treatment*Time] Plant+Plant.Time VSTRUCTURE [TERMS=Plant+Plant.Time; FORMATION=whole;\ CORRELATE=unrestricted; CINITIAL=CovRCR['Across terms']] REML [PRINT=components,deviance] Height " Include a random cubic spline term over time." VCOMPONENTS [FIXED=Treatment*Time; SPLINE=Time] Plant+Plant.Time VSTRUCTURE [TERMS=Plant+Plant.Time; FORMATION=whole;\ CORRELATE=unrestricted; CINITIAL=CovRCR['Across terms']] REML [PRINT=model,components,deviance] Height " Plot the spline term." VKEEP Time; SPLBLUP=Tblup; SPLDESIGN=Tdes; SPLX=Tknot CALCULATE Tspline = Tdes[1] *+ Tblup[1] YAXIS 1; LOWER=-50; UPPER=50 PEN 1,2; METHOD=line DGRAPH [TITLE='Common spline effect over time'] Tspline; Tknot[1] " Fit separate splines for each treatment." VCOMPONENTS [FIXED=Treatment*Time; SPLINE=Treatment.Time]\ Plant+Plant.Time VSTRUCTURE [TERMS=Plant+Plant.Time; FORMATION=whole;\ CORRELATE=unrestricted; CINITIAL=CovRCR['Across terms']] REML [PRINT=components,deviance] Height " Plot the splines." VKEEP [RMETHOD=all; RESIDUALS=R1] VKEEP [RMETHOD=notspline; RESIDUALS=R2] CALCULATE TTspline = R1 - R2 DGRAPH [TITLE='Separate treatment splines'] TTspline; Time;\ PEN=Treatment