The problem with the multivariate procedure outlined in the above is that it makes no assumptions regarding the temporal correlation structure of the data, and hence, may be overparameterized leading to poor parameter estimates. The mixed model procedure allows us to look at temporal correlation functions involving a limited number of parameters. The mixed model procedure falls beyond the scope of this class. The following brief outline is intended to be just an overview.
Approach 3 - Mixed Model Analysis
The mixed model initially looks identical to the split-plot model considered earlier.
\(Y_{ijk} = \mu + \alpha_i + \beta_{j(i)}+ \tau_k + (a\tau)_{ik} + \epsilon_{ijk}\)
where
- \(\mu\) = overall mean
- \(\alpha _ { i }\) = effect of treatment i
- \(\beta j ( i )\) = random effect of dog j receiving treatment i
- \(\tau_k\) = effect of time k
- \((a\tau)_{ik}\) = treatment by time interaction
- \(\epsilon_{ijk}\) = experimental error
Assumptions
- The dog effects \(\beta_{j ( i )}\) are independently sampled from a normal distribution with mean 0 and variance \(\sigma^2_\beta\) .
- The errors \(\epsilon_{ijk}\) from different dogs are independently sampled from a normal distribution with mean 0 and variance \(\sigma^2_\epsilon\) .
- The correlation between the errors for the same dog depends only on the difference in observation times:\(|k-k'|\)
Several covariance and correlation functions are listed below.
- Compound Symmetry: \(cov(\epsilon_{ijk}, \epsilon_{ijk'}) = \sigma^2_\epsilon+\sigma^2_\beta\) if \(k=k'\) and \(\sigma^2_\beta\) otherwise. This is the default structure for split plots.
- Autoregressive: \(corr(\epsilon_{ijk}, \epsilon_{ijk'}) = \rho^{|k-k'|}\)
- Autoregressive Moving Average:
\(corr(\epsilon_{ijk}, \epsilon_{ijk'}) = \left\{\begin{array}{cl}\gamma; & \text{if } |k-k'| = 1 \\ \gamma\rho^{|k-k'|-1}; & \text{if } |k-k'| \ge 2\end{array}\right.\)
- Toeplitz: \(corr(\epsilon_{ijk}, \epsilon_{ijk'}) = \rho(|k-k'|)\)
- The autoregressive model is a special case of a autoregressive moving average model with γ = 1.
- The autoregressive moving average model is a special case of a toeplitz model with
\(\rho(|k-k'|) = \left\{\begin{array}{cl}\gamma; & \text{if } |k-k'| = 1 \\ \gamma\rho^{|k-k'|-1}; & \text{if } |k-k'| \ge 2\end{array}\right.\)
Analysis
Approach 1: If one model is a special case of another, they can be compared using the -2 log likelihood values output. The difference is approximately chi-squared with degrees of freedom equal to the difference between the numbers of estimated parameters. For example, when comparing the AR(1) model with the ARMA(1,1) model, the difference between their -2 log likelihood values is:
\(237.426 - 237.329 = .097\)
which is less than the chi-square critical value of \(3.85 = \chi^2_{1, 0.05}\) (the df is 1 because there is one additional parameter estimated with the ARMA(1,1) model). This would not be significant evidence to claim the ARMA(1,1) fits better and the AR(1) model would be preferred.
Approach 2: Models that are not special cases of each other can be compared using AICC or BIC values from the output. Smaller values are better. For example, based on the AICC values for the CS, AR(1), ARMA(1,1), and Toeplitz models below, the AR(1) would be preferred.
Compound Symmetry | 243.9 |
AR(1) | 243.6 |
ARMA(1,1) | 245.7 |
Toeplitz | 247.8 |
Using SAS
The syntax for using SAS's Mixed Model procedure can be seen in the program below.
Download the SAS Program here: dog3.sas
options ls=78;
title "Mixed Model Analysis - Dog Data";
data dogs;
infile "D:\Statistics\STAT 505\data\dog1.txt";
input treat dog p1 p2 p3 p4;
time=1; k=p1; output;
time=5; k=p2; output;
time=9; k=p3; output;
time=13; k=p4; output;
drop p1 p2 p3 p4;
run;
proc mixed;
class treat dog time;
model k=treat|time;
random dog(treat);
lsmeans treat|time;
run;
proc mixed;
class treat dog time;
model k=treat|time;
random dog(treat);
repeated / subject=dog(treat) type=ar(1);
run;
proc mixed;
class treat dog time;
model k=treat|time;
random dog(treat);
repeated / subject=dog(treat) type=arma(1,1);
run;
proc mixed;
class treat dog time;
model k=treat|time;
random dog(treat);
repeated / subject=dog(treat) type=toep;
run;
The general format here for the mixed model procedure requires that the data are on separate lines for separate points in time.
Except for the first model, each of the various models will have repeated statements. The second, third and fourth models contain the repeated statements where subject is specified to be the dog within treatments, indicating within which units we have our repeated measures, in this case within each of the dogs.
This is followed by the type option which specifies what model you want. Here we set ar(1) for an autoregressive model., arma(1,1) for a 1,1 autoregressive moving average model and toep, short for a Toeplitz model.
Based on the AR(1) model above, the following hypothesis tests are obtained:
Effect | F | d.f. | p-value |
Treatments | 6.09 | 3,32 | 0.0021 |
Time | 9.84 | 3,96 | < 0.0001 |
Treatment by Time | 1.89 | 9,96 | 0.0631 |