Help the Stat Consulting Group by
giving a gift
Loading
|
|
|
||||
|
Help the Stat Consulting Group by
giving a gift
| |||||
|
Loading
|
|||||
This page will examine how to analyze data from a repeated measures Analysis of Variance. We will use data based on a real consulting problem we received. (We are grateful to Dr. Holly Hazlett-Stevens for giving us permission to use and adapt her data.) This page is not a comprehensive examination of how to analyze repeated measures data, but is a summary of the advice we would give for this particular analysis. The most common statistical procedures we see people use for repeated measures data are SAS proc glm and SPSS glm, but also SAS proc mixed and SPSS manova can be useful for such analyses, and can offer analyses not possible in SAS proc glm and SPSS glm. We will use this case study to illustrate some of the advice we gave in how to choose a proper analysis, how to interpret the results, and how to choose among statistical procedures within a package, and between packages to get the desired analysis.
This example is based on an experiment on 42 subjects who were randomly assigned to one of three groups (conditions, or cond), control (1), treatment A (2) and treatment B (3). There were three dependent measures (ib ms and sr) and they were each measured at five time points (e.g., ib1 ms1 and sr1 were all measured at time 1 and ib5 ms5 and sr5 were all measured at time 5). So this experiment has one between subjects variable (cond) and one repeated measures variable (time). The data are shown in the file repeat.txt).
The study predicts that the responses of the control group will not differ from those who received treatment A. However, it is predicted that the responses of those who receive treatment B should differ from the control group and the treatment A group. In particular, it is predicted that the responses of the treatment B group will tend to oscillate (go up and down) over the five trials (time periods) more than the other two groups. In other words, it is predicted that there will be an interaction between cond and time.
You can see and/or download the file repeat.txt that is referred to in the SAS program below.
options nocenter; data repeat; infile "repeat.txt"; input id cond ib1-ib5 ms1-ms5 sr1-sr5 ; run; proc print data=repeat(obs=10); run;
The output of the proc print is shown below.
Obs id cond ib1 ib2 ib3 ib4 ib5 1 1 1 1415.94 1457.06 1484.38 1583.71 1603.07 2 2 3 1385.95 1550.62 1528.70 1541.53 1582.64 3 3 2 1211.38 1313.07 1326.06 1294.30 1411.97 4 4 2 1318.06 1350.23 1423.97 1386.87 1392.59 5 5 3 1648.54 1655.81 1672.03 1756.02 1758.25 6 6 3 1467.75 1506.05 1595.74 1675.13 1570.27 7 7 1 1627.09 1659.87 1734.72 1736.64 1756.33 8 8 3 1391.31 1362.23 1381.47 1585.56 1553.41 9 9 3 999.29 1127.63 1295.83 1400.72 1382.79 10 10 2 1600.34 1645.25 1665.01 1717.93 1722.37 Obs id cond ms1 ms2 ms3 ms4 ms5 sr1 sr2 sr3 sr4 sr5 1 1 1 64.71 51.04 54.40 86.96 90.23 8 7 5 5 6 2 2 3 36.51 50.29 47.32 43.08 46.45 14 15 17 17 16 3 3 2 61.08 71.80 76.45 72.78 79.32 18 18 18 16 17 4 4 2 28.75 29.64 39.95 31.05 33.89 16 15 11 13 14 5 5 3 88.97 90.67 87.73 91.87 104.56 9 11 10 12 11 6 6 3 71.59 66.88 78.89 100.50 76.31 13 15 14 11 9 7 7 1 109.85 97.48 122.51 122.40 139.64 10 11 11 9 9 8 8 3 . . . . . 15 16 13 13 11 9 9 3 22.20 48.71 57.71 61.36 76.34 13 10 4 11 10 10 10 2 78.60 69.38 74.36 79.24 82.49 13 11 13 13 11
Let's do a standard between/within repeated measures ANOVA using SAS proc glm, with cond as a between subjects variable, and time as a within subjects variable.
PROC GLM DATA=repeat; class cond; model sr1-sr5 = cond ; repeated time 5 ; Run;
You can see the entire output called an1.txt, and we have shown important excerpts below.
Manova Test Criteria and Exact F Statistics for the Hypothesis of no time Effect
H = Type III SSCP Matrix for time
E = Error SSCP Matrix
S=1 M=1 N=17
Statistic Value F Value Num DF Den DF Pr > F
Wilks' Lambda 0.64788274 4.89 4 36 0.0030
Pillai's Trace 0.35211726 4.89 4 36 0.0030
Hotelling-Lawley Trace 0.54348918 4.89 4 36 0.0030
Roy's Greatest Root 0.54348918 4.89 4 36 0.0030
Manova Test Criteria and F Approximations for the Hypothesis of no time*cond Effect
H = Type III SSCP Matrix for time*cond
E = Error SSCP Matrix
S=2 M=0.5 N=17
Statistic Value F Value Num DF Den DF Pr > F
Wilks' Lambda 0.62780347 2.36 8 72 0.0259
Pillai's Trace 0.39613669 2.28 8 74 0.0303
Hotelling-Lawley Trace 0.55472196 2.46 8 49.161 0.0253
Roy's Greatest Root 0.47432775 4.39 4 37 0.0053
NOTE: F Statistic for Roy's Greatest Root is an upper bound.
NOTE: F Statistic for Wilks' Lambda is exact.
Repeated Measures Analysis of Variance Tests of Hypotheses for Between Subjects Effects Source DF Type III SS Mean Square F Value Pr > F cond 2 573.895238 286.947619 5.70 0.0067 Error 39 1962.528571 50.321245
Repeated Measures Analysis of Variance
Univariate Tests of Hypotheses for Within Subject Effects
Adj Pr > F
Source DF Type III SS Mean Square F Value Pr > F G - G H - F
time 4 88.1142857 22.0285714 8.36 <.0001 0.0001 <.0001
time*cond 8 26.6285714 3.3285714 1.26 0.2669 0.2851 0.2807
Error(time) 156 411.2571429 2.6362637
Greenhouse-Geisser Epsilon 0.6496
Huynh-Feldt Epsilon 0.7358
Let's examine these results, especially the within subjects results. Notice that the output shows two separate tests of the within subjects effects: one set of tests labeled the MANOVA tests and the other the UNIVARIATE tests. Now look at the results for the time by cond interaction. The UNIVARIATE tests give a p value of 0.2669 (and adjusted values of .2851 and .2807, but we will not be so concerned with these). The MANOVA tests give 4 p-values, repeated below.
Wilks' Lambda 0.62780347 2.36 8 72 0.0259 Pillai's Trace 0.39613669 2.28 8 74 0.0303 Hotelling-Lawley Trace 0.55472196 2.46 8 49.161 0.0253 Roy's Greatest Root 0.47432775 4.39 4 37 0.0053 NOTE: F Statistic for Roy's Greatest Root is an upper bound. NOTE: F Statistic for Wilks' Lambda is exact.
These UNIVARIATE and MANOVA tests give quite different results, and when this occurs, the MANOVA tests should be used. The UNIVARIATE results assume homogeneity of variance across the trials (that the variance of the scores on each of the five trials are the same) and homogeneity of covariance (that the covariances among all of the pairs of trials are the same). This is a very dicey assumption, and can be easily and frequently violated. By contrast, the MANOVA tests are based on a much more realistic (but possibly overly complex) assumption about the variances and covariances. The MANOVA tests model a separate variance for each trial (in this case, a separate variance for each of the five trials, and a separate covariance for each of the pairs of trials, i.e. 10 separate covariances). While this assumption may be overly complex and estimate many more variances and covariances than may be necessary, it is a more realistic assumption than the assumptions made by the UNIVARIATE tests.
As we compare the UNIVARIATE and MANOVA results, we are delighted that the MANOVA tests are the appropriate tests because those tests are significant. But which test should we report. Should we pick Roy's test because that has the best significance value? The output states that the Wilk's Lambda is exact, so we would pick that as the value to report. So, based on the Wilks' Lambda, we can say that there was a significant interaction between cond and time (p=.0259).
Let's see how we could do these analyses in SPSS.
data list free file="repeat.txt" /id cond ib1 to ib5 ms1 to ms5 sr1 to sr5. GLM sr1 to sr5 by cond /wsfactors = time 5.
Excerpts of the SPSS output is shown below (you can view the entire output at an1b.htm). Although the output looks different from the SAS output in its formatting, if you look more closely you can see that the results shown in the output are quite identical to those from SAS. Like the SAS output, the SPSS output shows Multivariate (MANOVA) tests for the within subjects effects, and Univariate tests for the within subjects effects.
The tests from Analysis 1 did assess whether there was an interaction between cond and time, but it did not assess this very precisely. There are many possible patterns of data that could produce such an interaction, many of which would not be consistent with the main hypothesis. To perform a more precise test, we will do a trend analysis of the time variable, and look at the interactions of these trend components of time by cond. We would predict that the higher orders of time (e.g., time^3 or time^4) would interact with cond (because we predict the Treatment B group will oscillate more over time than the other groups). This analysis is shown below. The polynomial / summary asks for the trend analysis on time.
PROC GLM DATA=repeat; class cond; model sr1-sr5 = cond ; repeated time 5 polynomial / summary ; Run;
The relevant results are shown below (you can see the entire results in an2.txt . The tests under time_1 refer to the linear effect of time. The results labeled Mean with F=19.78 indicate an overall linear effect of time. The results labeled cond with F=.69 indicate that the interaction of the linear effect of time by cond was not significant. Looking at the higher order effects (time_3 and time_4) we see that time_4 by cond is significant (F=4.96), indicating an interaction of time^4 by cond, as we predicted. We can think of time^4 as the tendency for the scores over time to exhibit 3 bends (like an M or like a W ), so a time^4 by cond interaction means that the conditions differ in their tendency to show 3 bends. This test more precisely reflects a test of our hypothesis than Analysis #1; however, it is still ambiguous which groups are showing the greater tendency to oscillate. In Analysis #3, we will more precisely compare among the groups as stated by our hypotheses. But first, let us show how this analysis can be obtained in SPSS.
Repeated Measures Analysis of Variance Analysis of Variance of Contrast Variables time_N represents the nth degree polynomial contrast for time Contrast Variable: time_1 Source DF Type III SS Mean Square F Value Pr > F Mean 1 80.6095238 80.6095238 19.78 <.0001 cond 2 5.6333333 2.8166667 0.69 0.5071 Error 39 158.9571429 4.0758242 Contrast Variable: time_2 Source DF Type III SS Mean Square F Value Pr > F Mean 1 7.4081633 7.4081633 1.86 0.1800 cond 2 6.7448980 3.3724490 0.85 0.4358 Error 39 154.9897959 3.9740973 Contrast Variable: time_3 Source DF Type III SS Mean Square F Value Pr > F Mean 1 0.00952381 0.00952381 0.01 0.9331 cond 2 2.74761905 1.37380952 1.03 0.3667 Error 39 52.04285714 1.33443223 Contrast Variable: time_4 Source DF Type III SS Mean Square F Value Pr > F Mean 1 0.08707483 0.08707483 0.08 0.7856 cond 2 11.50272109 5.75136054 4.96 0.0121 Error 39 45.26734694 1.16070120
Let's see how we could do these analyses in SPSS. If we examine the output from the first analysis (an1b.htm), we would see this output that we overlooked in Analysis #1. As you can see by comparing these results from those shown above in SAS, these are the same results. By default, SPSS includes a trend analysis on the repeated measures factor (in this case, time) saving us from having to specifically request this analysis.
The tests from Analysis 2 indicated that there were differences in the amount of oscillation over time by cond, but we are specifically interested in seeing if Treatment B causes more oscillation than the Control and Treatment A (and we also predict the amount of oscillation should be the same for Control and Treatment A). We can form comparisons on cond to make these specific comparisons. Ordinarily, you would think to use the contrast statement for this analysis; however, the contrasts would apply only to cond, but we want them to interact with time (and in particular, with the different trend components of time). We will create a new variable called cond1 that is 1 if it is Treatment B, and 0 if not. Comparisons that involve cond1 will directly compare Treatment B to the other 2 groups. Likewise, we will create a new variable called cond2 that will be 1 if Treatment A, 0 if Control, and missing if Treatment B. The variable cond2 directly compares Treatment A with the Control group. This is illustrated below.
data repeat2; set repeat; if cond = 1 or cond = 2 then cond1 = 1; if cond = 3 then cond1 = 2; if cond = 1 then cond2 = 1; if cond = 2 then cond2 = 2; if cond = 3 then cond2 = .; run; PROC GLM DATA=repeat2; class cond1; model sr1-sr5 = cond1 ; repeated time 5 polynomial / summary nouni ; Run; PROC GLM DATA=repeat2; class cond2; model sr1-sr5 = cond2 ; repeated time 5 polynomial / summary nouni ; Run; quit;
The results from the first proc glm are shown below.
The GLM Procedure Repeated Measures Analysis of Variance Analysis of Variance of Contrast Variables time_N represents the nth degree polynomial contrast for time Contrast Variable: time_1 Source DF Type III SS Mean Square F Value Pr > F Mean 1 59.2011905 59.2011905 14.87 0.0004 cond1 1 5.3440476 5.3440476 1.34 0.2535 Error 40 159.2464286 3.9811607 Contrast Variable: time_2 Source DF Type III SS Mean Square F Value Pr > F Mean 1 2.9600340 2.9600340 0.76 0.3878 cond1 1 6.4362245 6.4362245 1.66 0.2053 Error 40 155.2984694 3.8824617 Contrast Variable: time_3 Source DF Type III SS Mean Square F Value Pr > F Mean 1 0.00119048 0.00119048 0.00 0.9766 cond1 1 0.14404762 0.14404762 0.11 0.7471 Error 40 54.64642857 1.36616071 Contrast Variable: time_4 Source DF Type III SS Mean Square F Value Pr > F Mean 1 1.94710884 1.94710884 1.71 0.1984 cond1 1 11.23282313 11.23282313 9.87 0.0032 Error 40 45.53724490 1.13843112
Looking in the section time_4 we see the cond1 effect is significant (F=9.87). Below we show a graph of this effect. The green line (Treatment B) clearly shows a greater tendency to have 3 bends (like an M) than the red line (Treatment A and Control combined). This is exactly the pattern that was hypothesized, and the n time_4 we see the cond1 effect tests for this pattern of results.
The results from the second proc glm are shown below.
The GLM Procedure Repeated Measures Analysis of Variance Analysis of Variance of Contrast Variables time_N represents the nth degree polynomial contrast for time Contrast Variable: time_1 Source DF Type III SS Mean Square F Value Pr > F Mean 1 75.0892857 75.0892857 18.26 0.0002 cond2 1 0.2892857 0.2892857 0.07 0.7929 Error 26 106.9214286 4.1123626 Contrast Variable: time_2 Source DF Type III SS Mean Square F Value Pr > F Mean 1 13.5943878 13.5943878 3.06 0.0920 cond2 1 0.3086735 0.3086735 0.07 0.7941 Error 26 115.4540816 4.4405416 Contrast Variable: time_3 Source DF Type III SS Mean Square F Value Pr > F Mean 1 0.08928571 0.08928571 0.06 0.8160 cond2 1 2.60357143 2.60357143 1.61 0.2155 Error 26 42.00714286 1.61565934 time_N represents the nth degree polynomial contrast for time Contrast Variable: time_4 Source DF Type III SS Mean Square F Value Pr > F Mean 1 2.86989796 2.86989796 3.32 0.0801 cond2 1 0.26989796 0.26989796 0.31 0.5813 Error 26 22.50306122 0.86550235
It was predicted that there would be no differences in the trend over time between Treatment A and the Control group. The results shown above are consistent with these predictions -- the time1 by cond, time2 by cond, time3 by cond, and time4 by cond effects are all not significant. This is consistent with the graph shown below. The red line represents the control group, and the green line represents treatment A. While there may be minor differences in the trends shown by these two lines, the trends exhibited by these two lines are not significantly different.
Let's see how to do these analyses in SPSS. The SPSS program below creates cond1 and cond2 like in the SAS program, and uses glm to perform the analyses.
if cond = 1 or cond = 2 cond1 = 1. if cond = 3 cond1 = 2. if cond = 1 cond2 = 1. if cond = 2 cond2 = 2. if cond = 3 cond2 = $sysmis. GLM sr1 to sr5 by cond1 /wsfactors = time 5. GLM sr1 to sr5 by cond2 /wsfactors = time 5.
Excerpts of the results are shown below (you can see the entire results at an3b.htm. Excerpts from the first glm are shown below, and correspond to the SAS results.
Excerpts from the second glm are shown below, and they too correspond to the SAS results.
The tests from Analysis 3 allow us to assess each of the trend components. Perhaps instead we would like to assess the linear trend and the nonlinear (i.e. quadratic + cubic + order 4) trend. We could assess the non-linear trend by using the information from the prior analyses to compute the SSnon-linear by summing SSquad+SScubic+SSorder4. We could then compute the MSnonlinear by taking SSnonlinear / 3. For example, for the first analysis, the SSnonlinear would be 6.436 + .144 + 11.233 = ??? . Then the MSnonlinear would be ??? / 3. We would then divide this by the Error term. Whoops... there are different error terms for each of these effects. It is not clear what the error term would be. Here is another way you could do this, using SPSS manova (note that we use SPSS manova because we don't believe this kind of analysis is possible using SAS proc glm or SPSS glm).
Let's first illustrate how you can do analysis 3 using SPSS manova, as shown below.
MANOVA sr1 to sr5 by cond(1,3) /WSFACTORS = time(5) /CONTRAST(time) = POLYNOMIAL /CONTRAST(cond)=SPECIAL(1 1 1 1 1 -2 1 -1 0) /WSDESIGN=time(1) time(2) time(3) time(4) /DESIGN=cond(1) cond(2).
You can see the entire results in an3c.htm. Here is the cond by time(order 4)interaction. The cond(1) by time(4) is the first analysis, comparing Control and Treatment A versus Treatment B by time(order 4). You can see these results match the SAS and SPSS results. (You may note that the cond(2) by time(4) has the same SS, but a different F value. This is because this test uses the error term from all of the data, whereas the SAS and SPSS analyses use an error term based only on the data from Treatment A and Control.)
Tests involving 'TIME(4)' Within-Subject Effect. Tests of Significance for T5 using UNIQUE sums of squares Source of Variation SS DF MS F Sig of F WITHIN+RESIDUAL 45.27 39 1.16 TIME(4) .09 1 .09 .08 .786 COND(1) BY TIME(4) 11.23 1 11.23 9.68 .003 COND(2) BY TIME(4) .27 1 .27 .23 .632
Now, let's show how you can get the linear and nonlinear tests. We want to collect together the non-linear effects of time (the second, third, and fourth orders of time). We use the partition subcommand to divide up the effects into the linear effect (i.e., the first order) and the remaining 3 df (i.e., the second, third, and fourth orders of time, the non-linear effects). Now, when we refer to time(1) we are referring to the linear effect of time, and when we refer to time(2) we are referring to the non-linear effect of time. The commands are shown below.
MANOVA sr1 to sr5 by cond(1,3) /WSFACTORS = time(5) /CONTRAST(time) = POLYNOMIAL /PARTITION(time) = (1,3) /WSDESIGN=time(1) time(2) /CONTRAST(cond)=SPECIAL(1 1 1 1 1 -2 1 -1 0) /DESIGN=cond(1) cond(2).
You can see the entire analysis in an4.txt. Below, we show the test of greatest interest, comparing Treatment B versus Treatment A and Control by the non-linear effect of time. As we hypothesized, this result is significant. (Note that we used the MANOVA results, not the univariate results.)
* * * * * * A n a l y s i s o f V a r i a n c e -- design 1 * * * * * * EFFECT .. COND(1) BY TIME(2) Multivariate Tests of Significance (S = 1, M = 1/2, N = 17 1/2) Test Name Value Exact F Hypoth. DF Error DF Sig. of F Pillais .29242 5.09695 3.00 37.00 .005 Hotellings .41327 5.09695 3.00 37.00 .005 Wilks .70758 5.09695 3.00 37.00 .005 Roys .29242 Note.. F statistics are exact.
As we discussed in Analysis #1, SAS proc glm and SPSS glm give you two sets of results, Univariate and Multivariate (MANOVA) results. These results are based on different assumption about the structure of the covariance of the scores across time. The Univariate model uses a structure called Compound Symmetry which estimates the covariance as shown below. There is a single Variance (represented by s2) for all 5 of the trials, and there is a single covariance (represented by s1) for each of the pairs of trials.
s2 s1 s2 s1 s1 s2 s1 s1 s1 s2 s1 s1 s1 s1 s2
The Multivariate (MANOVA) model uses a structure called Unstructured which estimates the covariance as shown below. Each trial has its own variance (e.g., s12 is the variance of trial 1) and each pair of trials has its own covariance (e.g., s21 is the covariance of trial 1 and trial2).
s12
s21 s22
s31 s32 s32
s41 s42 s43 s42
s51 s52 s53 s54 s52
As you might imagine, these are not the only possible covariance structures, but they are the only covariance structures that are available to you if you use SAS proc glm, SPSS glm (or SPSS manova). However, SAS proc mixed allows you to select among a wide variety of covariance structures, allowing you to choose the appropriate covariance structure for your data. First, let's illustrate how to replicate the results from analysis #1 using proc mixed.
The first thing we need to do in order to use proc mixed is to get our data into the proper shape. In our prior analyses, our data was in a wide format where there was one record per subject. SAS proc mixed; however, wants the data in a long format, where there is one score per observation (yielding 5 observations per subject). We can use the program below to reshape the data from wide to long.
data replong ;
set repeat2 ;
array aib(5) ib1-ib5;
array ams(5) ms1-ms5;
array asr(5) sr1-sr5;
do time = 1 to 5;
ib = aib(time);
ms = ams(time);
sr = asr(time);
output;
end;
keep id cond time ib ms sr;
run;
proc print data=replong(obs=10);
run;
We see the output of the proc print below.
Obs id cond time ib ms sr 1 1 1 1 1415.94 64.71 8 2 1 1 2 1457.06 51.04 7 3 1 1 3 1484.38 54.40 5 4 1 1 4 1583.71 86.96 5 5 1 1 5 1603.07 90.23 6 6 2 3 1 1385.95 36.51 14 7 2 3 2 1550.62 50.29 15 8 2 3 3 1528.70 47.32 17 9 2 3 4 1541.53 43.08 17 10 2 3 5 1582.64 46.45 16
We can now use proc mixed to perform one of the Analyses from Analysis #1. This analysis assumes a compound symmetric covariance structure.
PROC MIXED DATA=replong; CLASS id cond time ; MODEL sr = cond time cond*time ; REPEATED time / SUBJECT=id TYPE=CS ; RUN;
The entire results can be seen in an5a.txt, and excerpts are shown below. If you compare these results to Analysis 1, you will see that the results match the Univariate results from SAS proc glm and SPSS glm. This makes sense, since all of these analyses are assuming the covariance is compound symmetric.
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
cond 2 39 5.70 0.0067
time 4 156 8.36 <.0001
cond*time 8 156 1.26 0.2669
Now, let's assume the covariance is Unstructured using proc mixed.
PROC MIXED DATA=replong; CLASS id cond time ; MODEL sr = cond time cond*time ; REPEATED time / SUBJECT=id TYPE=UN R RCORR ; RUN;
The entire results can be seen in an5b.txt, and excerpts are shown below. The results for cond*time are similar to the Multivariate (MANOVA) results from SAS proc glm and SPSS glm. It makes sense that the results should be similar since all of these analyses are assuming the covariance is Unstructured. It makes sense that the results are not exactly the same because the SAS proc glm and SPSS glm results are from Multivariate tests, whereas the SAS proc mixed results are Univariate tests, so the results would not be exactly the same.
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
cond 2 39 5.70 0.0067
time 4 39 5.30 0.0017
cond*time 8 39 2.70 0.0180
Let's further consider these two covariance structures. The compound symmetric covariance seems overly simplistic, and would likely be incorrect. The variance of the scores could change across trials (but the compound symmetric covariance assumes the variance is constant over trials). Also, the covariances among the trials might not all be the same, for example, the covariance may be stronger for trials that are closer together than for trials that are further apart (but the compound symmetric covariance assumes the covariances among trials are all the same).
If the compound symmetric covariance is overly simple, the unstructured covariance seems overly complex. In our example we have 5 trials, so the unstructured covariance has 5 variances and 10 covariances ( (5*4)/2 ), for a total of 15 variances and covariances being estimated. Our sample size is only 42, so we are estimating 15 variances/covariances based on just 42 subjects. We could be concerned that with such a small sample size, the variances/covariances may not be very stable, and the results would be based on these unstable estimates. Using proc mixed, we can seek a compromise, trying to find a covariance structure that strikes a balance between the simplicity of the compound symmetric covariance and the complexity of the unstructured covariance.
Let's start by showing how we can compare 2 covariance structures. If we return to the compound symmetry results (an5a.txt) it contains this information regarding the covariance structure.
Fit Statistics Res Log Likelihood -448.5 Akaike's Information Criterion -450.5 Schwarz's Bayesian Criterion -452.2 -2 Res Log Likelihood 897.0
Likewise, the unstructured results (an5b.txt) contains this information regarding the covariance structure.
Fit Statistics Res Log Likelihood -419.1 Akaike's Information Criterion -434.1 Schwarz's Bayesian Criterion -447.2 -2 Res Log Likelihood 838.3
We can use these Fit Statistics to assess these covariance structures. Larger values of Akaike's Information Criterion and larger values of Schwarz's Bayesian Criterion indicate better fit of the covariance structure, so the unstructured model has a better fit than the compound symmetry model but we would expect this given that it has so many more covariance parameters. Since these are nested models (one is a subset of the other) we can can test whether the additional parameters significantly improve the fit using the -2 Res Log Likelihood. We can take the differences in the Log Likelihood (897-838.3=58.7) and this is distributed as a chi-square with 15-2=13 df. This is significant even at the 0.01 level (exceeding 27.68), indicating that the additional parameters in the unstructured covariance matrix do lead to a better fit than the compound symmetric model.
Because we used the r corr options on the model statement, the covariance matrix and correlation matrix are included in the output, shown below.
Estimated R Matrix for id 1
Row Col1 Col2 Col3 Col4 Col5
1 7.5348 6.0385 6.8223 6.7802 7.1300
2 6.0385 9.9377 11.2564 9.6062 9.1410
3 6.8223 11.2564 16.5733 13.5586 12.6832
4 6.7802 9.6062 13.5586 13.4103 12.3535
5 7.1300 9.1410 12.6832 12.3535 13.4103
Estimated R Correlation Matrix for id 1
Row Col1 Col2 Col3 Col4 Col5
1 1.0000 0.6978 0.6105 0.6745 0.7093
2 0.6978 1.0000 0.8771 0.8321 0.7918
3 0.6105 0.8771 1.0000 0.9095 0.8508
4 0.6745 0.8321 0.9095 1.0000 0.9212
5 0.7093 0.7918 0.8508 0.9212 1.0000
Let's see if we can find a more parsimonious covariance structure than un. The top matrix is the covariance matrix, and we can see that the variances on the diagonal don't appear to be the same (7.5 9.9 16.5 13.4 and 13.4). This is one of the reasons the cs model did not fit well. Also, if we look at the correlations, we see a little bit of a pattern that the correlations are stronger the more adjacent the scores (e.g., trial4 and trial5 is .92, where trial1 and trial5 is .71). This looks like it could be modeled with a autoregressive covariance structure with heterogeneous variances. A traditional autoregressive structure should be familiar to those who have used time series analysis. The ar(1) structure is shown below. The variance is constant across trials. The covariance depends on the number of "steps" between trials, if there is one step, then it is sr and if it is two steps the covariance is sr2 and so forth.
s2
sr s2
sr2 sr s2
sr3 sr2 sr s2
sr4 sr3 sr2 sr s2
The autoregressive with heterogeneous variances model is like the one above, except instead of a common variance down the diagonal, there are separate variances down the diagonal. Looking at the variances from the data is why we chose the heterogeneous variances, because the variances are so different. The SAS proc mixed program for this is shown below.
PROC MIXED DATA=replong; CLASS id cond time ; MODEL sr = cond time cond*time ; REPEATED time / SUBJECT=id TYPE=ARH(1) ; RUN;
The results are shown below (and can be seen in an5c.txt ).
Fit Statistics Res Log Likelihood -432.3 Akaike's Information Criterion -438.3 Schwarz's Bayesian Criterion -443.5 -2 Res Log Likelihood 864.6
We can summarize the fit statistics in the table below.
Cov Parms AIC SBC -2RLL CS 2 -450.5 -452.2 897.0 UN 15 -434.1 -447.2 838.3 ARH 6 -438.3 -443.5 864.6
Comparing the CS and ARH models, we can see that the ARH model has a smaller AIC, SBC, and lower -2RLL. Previously, we tested the CS versus UN and found the UN to fit better. Likewise, we can compare the UN to the ARH structure by taking the difference in the -2RLL (864.6 - 838.3 = 26.3) and testing this against a Chi-Square distribution with 15-6=9 df, which is significant even at the 0.01 level (exceeding 21.66), indicating that the additional parameters in the unstructured covariance make for a significantly better fit than the ARH model. Perhaps there is a simpler covariance structure, but for these data, the unstructured covariance matrix seems to be the best fit.
Let's try analyzing the variable ib focusing on selecting the appropriate covariance structure. Let's start by using an un covariance structure and show the covariance.
PROC MIXED DATA=replong; CLASS id cond time ; MODEL ib = cond time cond*time ; REPEATED time / SUBJECT=id TYPE=UN R RCORR; RUN;
The fit statistics and covariances are shown below.
Fit Statistics
Res Log Likelihood -1137.8
Akaike's Information Criterion -1152.8
Schwarz's Bayesian Criterion -1165.8
-2 Res Log Likelihood 2275.6
Estimated R Matrix for id 1
Row Col1 Col2 Col3 Col4 Col5
1 47771 43575 42022 43310 40147
2 43575 44096 42149 43455 40317
3 42022 42149 43516 43266 40567
4 43310 43455 43266 46672 42137
5 40147 40317 40567 42137 41551
Estimated R Correlation Matrix for id 1
Row Col1 Col2 Col3 Col4 Col5
1 1.0000 0.9494 0.9217 0.9172 0.9011
2 0.9494 1.0000 0.9622 0.9579 0.9419
3 0.9217 0.9622 1.0000 0.9601 0.9540
4 0.9172 0.9579 0.9601 1.0000 0.9569
5 0.9011 0.9419 0.9540 0.9569 1.0000
Looking at the covariance matrix, it seems that a cs structure might fit this data, so let's try a cs structure and compare it to the un.
PROC MIXED DATA=replong; CLASS id cond time ; MODEL ib = cond time cond*time ; REPEATED time / SUBJECT=id TYPE=CS ; RUN;
Here are the fit statistics for the cs structure.
Fit Statistics Res Log Likelihood -1149.9 Akaike's Information Criterion -1151.9 Schwarz's Bayesian Criterion -1153.6 -2 Res Log Likelihood 2299.7
Looking at the covariance matrix, it seems that a cs structure might fit this data, so let's try a CS structure and compare it to the un. We take the difference in the -2 Res Log Likelihoods (2299.7 - 2275.6 = 24.1) and take the difference in the number of parameters (15-2=13). This is significant at the 0.05 level (since it exceeds 22.36), indicating that the improvement in fit in the UN model is significant, although just barely significant at the 0.05 level.
Perhaps an AR(1) model would fit better without needing so many parameters.
PROC MIXED DATA=replong; CLASS id cond time ; MODEL ib = cond time cond*time ; REPEATED time / SUBJECT=id TYPE=AR(1) ; RUN;
Here are the fit statistics for the AR(1) model. The fit this model (-2RLL=2297.4) is not much better than the CS model (-2RLL=2299.7).
Fit Statistics Res Log Likelihood -1148.7 Akaike's Information Criterion -1150.7 Schwarz's Bayesian Criterion -1152.4 -2 Res Log Likelihood 2297.4
With the very high correlations among the scores, perhaps an ARMA (autoregressive moving average) model would fit better than the AR model. We try the ARMA model below.
PROC MIXED DATA=replong; CLASS id cond time ; MODEL ib = cond time cond*time ; REPEATED time / SUBJECT=id TYPE=ARMA(1,1) ; RUN;
The fit statistics are shown below.
Fit Statistics Res Log Likelihood -1142.3 Akaike's Information Criterion -1145.3 Schwarz's Bayesian Criterion -1147.9 -2 Res Log Likelihood 2284.7
The ARMA model fits much better than the CS model. The ARMA model has 3 parameters, compared to the CS model with 2 parameters. The difference in -2 Res Log Likelihood is (2299.7 - 2284.7 = 15) which is significant at the 0.001 level (exceeding 10.82). Comparing the ARMA model to the UN model shows that the fit of the ARMA model is not significantly worse than the UN model (2284.7-2275.6=9.1 with 15-3=12 df, p > .5). For the ib variable, the ARMA covariance structure seems to give the best fit, a better fit than the CS model, and it has fewer parameters than the UN model and does not fit significantly worse than the UN model. So, the ARMA model seems to be yield the most parsimonious fit. In this case, proc mixed offers a considerable advantage over SAS proc glm and SPSS glm, because those tools only allow the CS and UN structures.
UCLA Researchers are invited to our Statistical Consulting Services
We recommend others to our list of Other Resources for Statistical Computing Help
These pages are Copyrighted (c) by UCLA Academic Technology Services