9.7 - Approach 3: Mixed Model Analysis

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

  1. The dog effects \(\beta_{j ( i )}\) are independently sampled from a normal distribution with mean 0 and variance \(\sigma^2_\beta\) .
  2. The errors \(\epsilon_{ijk}\) from different dogs are independently sampled from a normal distribution with mean 0 and variance \(\sigma^2_\epsilon\) .
  3. 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'|)\) 
Note!
  • 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.

Note! The first instance of the mixed procedure (without the repeated statement) is using the compound symmetry structure.

Download the SAS Program here: dog3.sas

SAS program for dog3 data - mixed model analysis

 

 

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