R Tutorials Top Banner



Model Formulae

If you have not yet read the Model Formulae tutorial, you should do so now.

Some Explanations

The naming of experimental designs has led to a great deal of confusion on the part of students as well as tutorial writers. The first design I will discuss in this tutorial is the single factor within subjects design, also called the single factor repeated measures design. This is a design with one response variable, and each "subject" or "experimental unit" is measured multiple times on that response variable. For example, pre-post designs, in which a group of subjects is measured both before and after some treatment, fit into this category. This design is also referred to in some sources as a "treatment × subjects" design (read "treatment by subjects"). You can also think of the design as having treatments nested within subjects. Or you could also think of it as a block design in which subjects correspond to blocks. You see my point.

In my field (psychology), we are used to referring to our experimental units as "subjects," because they are almost always either human beings or animals. I will retain that terminology in this tutorial. This is why, in the first study discussed below, grocery items are referred to as subjects--because the groceries correspond to what would ordinarily be human subjects in a typical psychology experiment. For a change, it is people in other fields who will have to adjust to strange terminology!

Single Factor Designs With Repeated Measures

Several years ago a friend and I decided to compare prices in local grocery stores. We made a list of ten representative grocery items and then went to four local stores and recorded the price on each item in each store. The most convenient way to table the data--and the first way that comes to mind--is as follows (numbers are prices in dollars)...

subject			storeA	storeB	storeC	storeD
lettuce			1.17	1.78	1.29	1.29
potatoes		1.77	1.98	1.99	1.99		
milk			1.49	1.69	1.79	1.59
eggs			0.65	0.99	0.69	1.09
bread			1.58	1.70	1.89	1.89
cereal			3.13	3.15	2.99	3.09
ground.beef		2.09	1.88	2.09	2.49
tomato.soup		0.62	0.65	0.65	0.69
laundry.detergent	5.89	5.99	5.99	6.99
aspirin			4.46	4.84	4.99	5.15
Excuse for the moment the absence of nicities in this table, like stubs and spanners, because shortly you will be happy I did the table this way!

First, I remember when a can of tomato soup was 12¢, but that's another story. We should notice laundry detergent costs a lot more than soup. This isn't interesting to us. In other words, there is a great deal of variability--perhaps 99% or more of the total variability--attributable entirely to individual differences between subjects, in which we are not the least bit interested. If we treated this as a straightforward oneway randomized ANOVA, all of this variability would go into the error term, and any variability due to store would be totally swamped. The advantage of the repeated measures analysis is that it allows us to parcel out variability due to subjects.

Second, we are now faced with getting these data into R. Try this...

> url = "http://ww2.coastal.edu/kingw/statistics/R-tutorials/text/groceries.txt"
> groceries = read.table(url, header=T)
> groceries
Whether or not that will work will depend upon whose servers are up and how much access your server is granted to coastal.edu. Ordinarily, you should be able to read data tables from a website, but who knows how these things sometimes work? The facilities to do so are provided in R.

You can also try it this way. Click on this link, and when the data table comes up in your browser, copy and paste it to a text editor, save the file as "groceries.txt", drop it into your working directory, and then...

> groceries = read.table("groceries.txt", header=T)
> groceries
Or, you can go back to that link, copy the entire table, and then paste it at the "0:" prompt...
> groceries = read.table(stdin(), header=T)
0: subjectstoreAstoreBstoreCstoreD
1: lettuce1.171.781.291.29
2: potatoes1.771.981.991.99
3: milk1.491.691.791.59
4: eggs0.650.990.691.09
5: bread1.581.701.891.89
6: cereal3.133.152.993.09
7: ground.beef2.091.882.092.49
8: tomato.soup0.620.650.650.69
9: laundry.detergent5.895.995.996.99
10: aspirin4.464.844.995.15
> groceries
1            lettuce1.171.781.291.29
2           potatoes1.771.981.991.99
3               milk1.491.691.791.59
4               eggs0.650.990.691.09
5              bread1.581.701.891.89
6             cereal3.133.152.993.09
7        ground.beef2.091.882.092.49
8        tomato.soup0.620.650.650.69
9  laundry.detergent5.895.995.996.99
10           aspirin4.464.844.995.15
As you can see, that didn't work for me because R didn't pick up the tabs as being white space. What can I tell ya? Sometimes that works, and sometimes it doesn't, and I don't know why! More R idiosyncrasies! (Hey! It's free!) If the same thing happened to you, try this. Copy and paste this script into R (WARNING: this will erase all items in your working directory)...
### start copying with this line
groceries = data.frame(
colnames(groceries) = c("subject","storeA","storeB","storeC","storeD")
### end copying with this line
And if that didn't do it, well, start typing!

This is not a proper data frame. Why not? Because our response variable--our ONE response variable--is the price of the items, and this is listed in four columns. Very bad! We need each variable in ONE column of the data frame. This sort of problem occurs so often, you'd think some nice R coder would have written a utility to rearrange such tables. (Remember the pleasure we had rearranging the "anorexia" data frame in an earlier tutorial?) Well, almost...

> stack(groceries)                     # Funny! :^)
I have not shown the output, but this is ALMOST what we want. It would have worked beautifully if this were a between subjects design. Here it is but a start...
> groceries2 = stack(groceries)
> subject = rep(groceries$subject,4)        # create the "subject" variable
> groceries2[3] = subject                   # add it to the new data frame
> rm(subject)                               # clean up your workspace
> colnames(groceries2) = c("price", "store", "subject")    # rename the columns
> groceries2                                               # take a look
   price  store           subject
1   1.17 storeA           lettuce
2   1.77 storeA          potatoes
3   1.49 storeA              milk
4   0.65 storeA              eggs
5   1.58 storeA             bread
6   3.13 storeA            cereal
7   2.09 storeA       ground.beef
8   0.62 storeA       tomato.soup
9   5.89 storeA laundry.detergent
10  4.46 storeA           aspirin
11  1.78 storeB           lettuce
12  1.98 storeB          potatoes
13  1.69 storeB              milk
14  0.99 storeB              eggs
15  1.70 storeB             bread
16  3.15 storeB            cereal
17  1.88 storeB       ground.beef
18  0.65 storeB       tomato.soup
19  5.99 storeB laundry.detergent
20  4.84 storeB           aspirin
21  1.29 storeC           lettuce
22  1.99 storeC          potatoes
23  1.79 storeC              milk
24  0.69 storeC              eggs
25  1.89 storeC             bread
26  2.99 storeC            cereal
27  2.09 storeC       ground.beef
28  0.65 storeC       tomato.soup
29  5.99 storeC laundry.detergent
30  4.99 storeC           aspirin
31  1.29 storeD           lettuce
32  1.99 storeD          potatoes
33  1.59 storeD              milk
34  1.09 storeD              eggs
35  1.89 storeD             bread
36  3.09 storeD            cereal
37  2.49 storeD       ground.beef
38  0.69 storeD       tomato.soup
39  6.99 storeD laundry.detergent
40  5.15 storeD           aspirin
NOW we have a proper data frame for this analysis. I belabor this point because it is a critical one!

At which store should we shop?

> with(groceries2, tapply(price, store, sum))
storeA storeB storeC storeD 
 22.85  24.65  24.36  26.26
For the items in the sample, it looks like storeA had the lowest prices. But as this is only a random sample of items we might have wanted to buy, we have to ask, will this difference hold up in general?

Once the data frame is set up correctly (all values of the response variable in ONE column, and another column holding the subject variable--i.e., which identifies which values of the response go with which subjects), the trick to doing a repeated measures ANOVA is in getting the model formula right. The Error term must reflect that we have "treatments nested within subjects." That is to say, in principle at least, we can see the "store" effect within each and every "subject" (grocery item). The ANOVA is properly done this way...

> aov.out = aov(price ~ store + Error(subject/store), data=groceries2)
> summary(aov.out)

Error: subject
          Df  Sum Sq Mean Sq F value Pr(>F)
Residuals  9 115.193  12.799

Error: subject:store
          Df  Sum Sq Mean Sq F value  Pr(>F)  
store      3 0.58586 0.19529  4.3442 0.01273 *
Residuals 27 1.21374 0.04495                  
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
"Store" and "subject" are our sources of variability. The treatment we are interested in is "store" (that's what we want to see the effect of), and this treatment effect is visible within each subject (i.e., nested within each subject). So the proper Error term is "subject/store", which is read as "store within subject." Notice once all that "subject" variability is parceled out, we do have a significant "store" main effect.

The Tukey HSD test will not run on this model object, so I propose pairwise t-tests with adjusted p-values to see which stores are significantly different (note: the output will be a table of p-values for each possible pairwise comparison)...

> with(groceries2, pairwise.t.test(price, store,
+                    p.adjust.method="holm", paired=T))

        Pairwise comparisons using paired t tests 

data:  price and store 

       storeA storeB storeC
storeB 0.17   -      -     
storeC 0.17   0.69   -     
storeD 0.07   0.49   0.33  

P value adjustment method: holm
The syntax for pairwise.t.test( ) is a bit out of line with other R functions, so this is how it works. The first argument is the response variable. The second argument is the grouping or explanatory variable. There is a comma between these two arguments--the function does not accept a formula. There are many p.adjust.methods available. See ?p.adjust for details. Since these should be paired (i.e., dependent) t-tests, I set the "paired=" option to TRUE. We see no significant differences in our post hoc tests, but storeA vs. storeD was at least close, p = .07.

The Friedman Rank Sum Test

There is a nonparametric version of the oneway ANOVA with repeated measures. It is executed this way...

> friedman.test(price ~ store | subject, data=groceries2)

        Friedman rank sum test

data:  price and store and subject 
Friedman chi-squared = 13.3723, df = 3, p-value = 0.003897
The formula reads "price as a function of store given subject". Here, "subject" is being treated as a blocking variable.

Two Way Mixed Factorial Designs

Clean out your workspace, and then enter the data for this analysis by copying and pasting this script into your R Console...

### begin copying with this line
# create the data vectors for the four conditions
# create the response variable vector
# create the test condition labels
# create the diet condition labels
# create the subject labels
# create the subject vector
# create the dataframe
# clean up
### end copying with this line
You should now have a data frame called "hill" in your workspace...
> ls()
[1] "groceries"  "groceries2" "hill"      
> str(hill)
'data.frame':   36 obs. of  4 variables:
 $ subject: Factor w/ 18 levels "A","B","C","D",..: 1 2 3 4 5 6 7 8 9 1 ...
 $ diet   : Factor w/ 2 levels "chicken","pasta": 1 1 1 1 1 1 1 1 1 1 ...
 $ test   : Factor w/ 2 levels "post","pre": 2 2 2 2 2 2 2 2 2 1 ...
 $ SSS    : num  18 13 18 15 22 32 31 24 15 15 ...
These data were collected by Ben Hill (by now, Ben Hill, Ph.D.) as part of his senior research project conducted in our department during the Fall 1999 semester. His interest was in the role of brain serotonin in creating sensation seeking behavior. Since he couldn't manipulate brain serotonin directly (our IRB would not permit it!), he chose to do so by way of diet, relying on the finding that brain serotonin is elevated after a meal rich in carbohydrates. Therefore, he had 18 subjects come to the lab. All subjects were given the Sensation Seeking Survey. Half the subjects were then given an evening meal consisting primarily of chicken, while the other half were given a meal consisting primarily of pasta. An hour later, the Sensation Seeking Survey was readministered. In the data frame, "subject" identifies the subject so that his or her "pre" and "post" eating scores can be kept paired for the analysis, "diet" identifies the type of meal the subject was given, "test" contains the pre-post meal information, and SSS contains the survey scores.

In this experiment, we have two explanatory variables, "diet", which is a between-subjects variable, and "test", which is a within-subjects variable. Thus, in lingo familiar to social scientists, the design is called a 2 × 2 mixed factorial design with repeated measures on "test". It's worth making sure you understand how the data frame is structured, so let's have a look...

> hill
   subject    diet test SSS
1        A chicken  pre  18
2        B chicken  pre  13
3        C chicken  pre  18
4        D chicken  pre  15
5        E chicken  pre  22
6        F chicken  pre  32
7        G chicken  pre  31
8        H chicken  pre  24
9        I chicken  pre  15
10       A chicken post  15
11       B chicken post  13
12       C chicken post  17
13       D chicken post  15
14       E chicken post  24
15       F chicken post  31
16       G chicken post  31
17       H chicken post  25
18       I chicken post  17
19       J   pasta  pre  17
20       K   pasta  pre  30
21       L   pasta  pre  18
22       M   pasta  pre  13
23       N   pasta  pre  23
24       O   pasta  pre  27
25       P   pasta  pre  27
26       Q   pasta  pre  24
27       R   pasta  pre  23
28       J   pasta post  19
29       K   pasta post  31
30       L   pasta post  18
31       M   pasta post  13
32       N   pasta post  24
33       O   pasta post  27
34       P   pasta post  26
35       Q   pasta post  28
36       R   pasta post  26
Note that all the SSS scores are in ONE column! The other columns completely identify the experimental conditions associated with that score, including which subject it came from.

By way of data summary...

> with(hill, tapply(SSS, list(diet,test), mean))
            post      pre
chicken 20.88889 20.88889
pasta   23.55556 22.44444
Drat! There goes R arranging our factor levels in alphabetical order again, making our means table look backwards. Let's pretty that up...
> hill$test = factor(hill$test, levels=c("pre","post"))
> with(hill, tapply(SSS, list(diet,test), mean))
             pre     post
chicken 20.88889 20.88889
pasta   22.44444 23.55556
Better! A nice graph might be in order as well...
> with(hill, boxplot(SSS ~ diet + test))    # output not shown
> with(hill, boxplot(SSS ~ test + diet))    # compare this one with the last one
> title(main="Ben Hill's SSS Data")
> title(ylab="SSS Scores")
Sensation Seeking
The SSS scores in the pasta group did go up as hypothesized, but there is a lot of within groups variability there, too. Will the effect turn out to be statistically significant?

Once again, the trick is in getting the model formula correct. In this case, we have two explanatory variables, and we want to see all possible main effects and interactions. The "test" variable is "within subjects"...

> aov.out = aov(SSS ~ diet * test + Error(subject/test), data=hill)
> summary(aov.out)

Error: subject
          Df  Sum Sq Mean Sq F value Pr(>F)
diet       1   40.11   40.11  0.5094 0.4857
Residuals 16 1259.78   78.74               

Error: subject:test
          Df  Sum Sq Mean Sq F value Pr(>F)
test       1  2.7778  2.7778  2.1739 0.1598
diet:test  1  2.7778  2.7778  2.1739 0.1598
Residuals 16 20.4444  1.2778
The effect we're looking for would be shown by the diet × test interaction. And sadly, it's not significant.

Other Designs

For reference, here are model formulae for a couple other common designs...

Two factor design with repeated measures on both factors:
DV ~ IV1 * IV2 + Error(subject/(IV1*IV2))

Three factor mixed design with repeated measures on IV2 and IV3:
DV ~ IV1 * IV2 * IV3 + Error(subject/(IV2*IV3))
And so on.

Return to the Table of Contents