Advanced Data Analysis, II

36-402  December 19, 1997  Spring, 1996

The education data

> reading <- read.table("reading.dat", head=T)
> reading
   Group Pre1 Post1
 1 Basal    4     5
 2 Basal    6     9
 3 Basal    9     5
 4 Basal   12     8
 5 Basal   16    10
 6 Basal   15     9
 7 Basal   14    12
 8 Basal   12     5
 9 Basal   12     8
10 Basal    8     7
11 Basal   13    12
12 Basal    9     4
....
> attach(reading, pos=1)
> ls()
[1] ".Last.value" "Group"       "Post1"       "Pre1"
> motif()

Now we want some plots of these data. First off a side-by-side boxplot is handy.

> boxplot(split(Pre1, Group))

and

> par(mfrow=c(2,2))

> qqnorm(split(Pre1,Group)$Basal)
> qqnorm(Pre1[Group=="DRTA"])
> qqnorm(Pre1[Group=="Strat"])

The bosplots seem to say that there might be some differences, and the qq plots (on 22 points each) look quite good.

Now for some simple summaries of the data.

> for (x in split(Pre1, Group)) {
+ print(summary(x)); print(var(x))}
 Min. 1st Qu. Median Mean 3rd Qu. Max.
    4    8.25   11.5 10.5      12   16
[1] 8.833333
 Min. 1st Qu. Median  Mean 3rd Qu. Max.
    6       8      9 10.14      12   19
[1] 11.171
 Min. 1st Qu. Median  Mean 3rd Qu. Max.
    4    6.25    8.5 9.136      12   14
[1] 11.171

The varinaces look about the same (the last two groups are exactly the same). Proceed to a simple model (one-way ANOVA).

> reading.aov <- aov(Pre1 ~ Group)
> summary(reading.aov)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
Group      2   21.9394 10.96970 1.055613 0.3540521
Residuals 63  654.6818 10.39177
> anova(reading.aov)

Analysis of Variance Table
Response: Pre1
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
Group      2   21.9394 10.96970 1.055613 0.3540521
Residuals 63  654.6818 10.39177

What about the post data?

> post.aov <- aov(Post1 ~ Group)
> summary(post.aov); anova(post.aov)
          Df Sum of Sq  Mean Sq  F Value       Pr(F)
Group      2  108.1212 54.06061 5.317437 0.007346757
Residuals 63  640.5000 10.16667
> Analysis of Variance Table
Response: Post1
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value       Pr(F)
Group      2  108.1212 54.06061 5.317437 0.007346757
Residuals 63  640.5000 10.16667

Some plots to help with this analysis.....

What did we just do?

So what have we just done.

Recall what the post part of the education data looks like. (I will not make a habit or re-writing data!).

Basal DRTA Strat
5  7  11
9  5  4
5  13  4
8  5  4
10  14  3
9  14  8
12  10  12
5  13  14
8  12  12
7  11  7
12  8  5
4  7  9
4  10  1
8  8  13
6  8  7
9  10  5
3  12  7
5  10  11
4  11  15
2  7  9
5  8  6
7  12  4

We calculated an ANOVA table as

Source d.f. SS MS F p
Group 2 108 54.06 5.31 0.0073
Residuals 63 640 10.16

“Read” this summary to say—”There appears to be a significant group effect” (Of course we still need to worry about residuals, influence, other diagnostics, but not at this instant).

For now, lets consider how we did the computations and what it all means.

Version 1

Y is the response variable (Post score) and Group is the indicator variable of the groups. We can proceed by doing just a simple linear regression like

Y = α + β_{1}X_{1} + β_{2}X_{2} + β_{2}X_{2} + e

where X_{1}, X_{2}, X_{3} are the indicator variables for the groups.

Of course, because of the linearity constraints, we can only do any two of X_{1}, X_{2}, X_{3}. Now, there are at least two ways of doing this computation.

Marginally by hand

> mymod <- lm(Post1 ~ Group)
> summary(mymod)
Call: lm(formula = Post1 ~ Group)
Residuals:
    Min     1Q Median    3Q   Max
 -6.773 -2.682 0.2273 2.295 7.227
Coefficients:
               Value Std. Error  t value Pr(>|t|)
(Intercept)   8.0758   0.3925    20.5762   0.0000
     Group1   1.5455   0.4807     3.2151   0.0021
     Group2  -0.1515   0.2775    -0.5460   0.5870
Residual standard error: 3.189 on 63 degrees of freedom
Multiple R-Squared: 0.1444
F-statistic: 5.317 on 2 and 63 degrees of freedom, the p-value is 0.007347
Correlation of Coefficients:
       (Intercept) Group1
Group1 0
Group2 0           0
> anova(mymod)
Analysis of Variance Table
Response: Post1
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value       Pr(F)
Group      2  108.1212 54.06061 5.317437 0.007346757
Residuals 63  640.5000 10.16667

Even more by hand

Create the dummy variables by hand,

   Group Post1 x1 x2 x3
 1 Basal     5  1  0  0
 2 Basal     9  1  0  0
 3 Basal     5  1  0  0
 4 Basal     8  1  0  0
 ....
23  DRTA     7  0  1  0
24  DRTA     5  0  1  0
25  DRTA    13  0  1  0
 ....
45 Strat    11  0  0  1
46 Strat     4  0  0  1
47 Strat     4  0  0  1
48 Strat     4  0  0  1
 ....

and do the analysis from there.

> mymod2 <- lm(Post1 ~ x1 + x2)
> summary(mymod2)
Call: lm(formula = Post1 ~ x1 + x2)
Residuals:
    Min     1Q Median    3Q   Max
 -6.773 -2.682 0.2273 2.295 7.227
Coefficients:
               Value Std. Error  t value Pr(>|t|)
(Intercept)   7.7727   0.6798    11.4339   0.0000
         x1  -1.0909   0.9614    -1.1347   0.2608
         x2   2.0000   0.9614     2.0804   0.0416
Residual standard error: 3.189 on 63 degrees of freedom
Multiple R-Squared: 0.1444
F-statistic: 5.317 on 2 and 63 degrees of freedom, the p-value is 0.007347
Correlation of Coefficients:
   (Intercept)      x1
x1 -0.7071
x2 -0.7071      0.5000
> anova(mymod2)
Analysis of Variance Table

Response: Post1
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value      Pr(F)
x1         1   64.1212 64.12121 6.307004 0.01460276
x2         1   44.0000 44.00000 4.327869 0.04156606
Residuals 63  640.5000 10.16667

Another Example of a contrast

The following data are drawn from Box, Hunter, and Hunter, Statistics for Experimenters. The data are from an experiment about blood coagulation times in animals. A group of 24 animals are given 4 different diets. The base question is “Is the coagulation time effected by diet”. The researchers believe, a-priori, that diets B and C are similar and (different to) diets A and D. Here are the data

A       B       C       D
62      63      68      56
60      67      66      62
63      71      71      60
59      64      67      61
        65      68      63
        66      68      63
                        59

Here are the first few lines of the reformated file, blood.dat

Time   Type
62      A
60      A
.....

Here is the start of an analysis.

> blood <- read.table("blood.dat", header=T)
> blood
   Time Type
 1   62    A
 2   60    A
......
> motif()

> plot.factor(blood)

It seems reasonable to assume that groups B and C are similar and that groups A and D are also similar (but different to B and C). A contrast which would check that is

−A + B + C − D

Lets do the simple anova with contrast.

> anova(aov( Time ~ Type) )
Analysis of Variance Table
Response: Time
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value        Pr(F)
Type       3  238.2857 79.42857 14.83708 3.238591e-05
Residuals 19  101.7143  5.35338
>

This analysis shows us that there is certainly a Type effect (very large F value). The next analysis does the 1 degree of freedom contrast to test the specific question.

> mymod <- aov( Time ~ C(Type, c(-1,1,1,-1),1) + Type)
> anova(mymod)
Analysis of Variance Table
Response: Time
Terms added sequentially (first to last)

                           Df Sum of Sq  Mean Sq  F Value     Pr(F)
C(Type, c(-1, 1, 1, -1), 1)  1  225.8182 225.8182 42.18233 0.0000032
Type                         2   12.4675   6.2338  1.16445 0.3333938
Residuals                   19  101.7143   5.3534

This shows a very compelling answer. Almost all of the Sum of Squares for the treatments is due to that one contrast (226 out of 238), and the F-test for the “rest” of the treatment effect is far from significant.

Finally, at the very least we should do some residual plots. Plotting residuals against fitted values is always a good idea, but in this case, there are only 4 fitted values. The following plots (which are automatically generated by one command) pretty much capture the situation.

> plot(mymod)

> plot.factor(data.frame(resid(mymod), Type))

Thats it for simple contrasts. We now need to worry about how multiple comparisons get in the way of reasonable inference.

An example in Splus and SAS

Multiple Comparison Methods

Recall that if we want a (family) of tests at the α level, then we form confidence intervals of the form

D̂ ± Z × ŝ(D̂)

All that differs in the various tests is the multiplier, Z. The following table summarized the various multipliers.

Test Critical Statistic Remarks
Single Comparison t[1 − α/2; ν]  ν df for Error
Tukey all pairs T = \frac{1}{\sqrt{2}}q[1 − α; r; ν] r groups
Scheffe’s method S^{2} = (r − 1)F [1 − α; r − 1; ν]
Bonferroni method B = t[1 − α/2g; ν]  g comparisons

A Randomized Block Design

The following data are found on page 28, of Planning of Experiments by D. R .Cox, Wiley, 1958, and were previously described in Cochran and Cox’s (a different Cox) book.

The treatments are five elves of potash (26, 54, 7, 108, 144 lb K_{2}O per acre applied to a cotton crop. On observation analyzed was a measure of single-fiber strength in arbitrary units, obtained as an average of a number of test on the cotton from each plot. There were three blocks each containing five plots. We have no details of the arrangement of plots within blocks. A summary of the data (which ignores the order in which the data were collected) is

T_{1} T_{2} T_{3} T_{4} T_{5}
Block 1 7.62 8.14 7.76 7.17 7.46
Block 2 8.00 8.15 7.73 7.57 7.68
Block 3 7.93 7.87 7.74 7.80 7.21

We are interested in knowing if the treatments differ. We are prepared to live with block to block variations, but don’t really care about them.

I entered the data into a computer file as follows,

Yield    Tmt    Block
7.62     T1     B1
8.14     T2     B1
7.76     T3     B1
7.17     T4     B1
7.46     T5     B1
8.00     T1     B2
8.15     T2     B2
7.73     T3     B2
7.57     T4     B2
7.68     T5     B2
7.93     T1     B3
7.87     T2     B3

7.74     T3     B3
7.80     T4     B3
7.21     T5     B3

and then read the data into Splus

> coxdat <- read.table("coxdata.dat", head=T)
> coxdat
   Yield Tmt Block
 1  7.62  T1    B1
 2  8.14  T2    B1
......

and now, to a simple analysis (assume I’ve presented all the plots—I did look at them and Treatment 3 is a little less variable than one might expect.).

> mymod <- aov(Yield ~ Block + Tmt)
> anova(mymod)
Analysis of Variance Table
Response: Yield
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
Block      2   0.09712 0.048560 1.111594 0.3749855
Tmt        4   0.73244 0.183110 4.191599 0.0403682
Residuals  8   0.34948 0.043685
>

It is quite clear that there is a reasonable treatment effect, and an ignorable Block effect (more on that in a moment), but which treatment is the best? In the absence of any particular interesting contrasts we might want to test all pairwise contrasts to see which treatments are better than the others. The plot shows us that the highest treatment level doesn’t seem to generate the highest yield.

The multiple comparison procedures

With 5 treatments (and hence 10 possible pairwise comparisons) we can’t do simple t-tests and the simple Bonferroni adjusts would be too severe, so lets try both the Tukey test and the Scheffe test. To get started we need to calculate the treatment means. (We could do this with a fancy use of contrasts, but brute force is just as good.) This is yet another one of those “cool Splus tricks”—you should look up the tapply and apply functions.

> mymeans<-tapply(Yield, Tmt, mean)
> mymeans
   T1       T2       T3       T4   T5
 7.85 8.053333 7.743333 7.513333 7.45

The rest of the calculation is essentially by hand. Lets compare level T2 (the biggest) with level T5 (the smallest). We expect this to be significant.

> mymeans[2] - mymeans[5]
        T2
 0.6033333

A 90% confidence interval for this difference, using Tukey’s studentized range intervals, is

D̂ ± Z × ŝ(D̂)

where

D = .60

,

ŝ(D̂) = \sqrt{MSE × (\frac{1}{3} + \frac{1}{3})} = \sqrt{0.043685 ∗ 2/3} = 0.17

and

Z = \frac{1}{\sqrt{2}}q[0.90; 5; 8] = 2.95

Hence the interval is

0.6 ± 2.95 × 0.17 = 0.6 ± 0.50

.

From the calculation it is clear that no pairwise difference that is less than 0.5 will be significant, which rules out everything else expect T2 against T4 or T5.

If we had done the test just as a pre-known contrast, using the t-test, the confidence 90% confidence region would have been

0.6 ± 1.86 × 0.17 = 0.6 ± 0.31

. Thus, in addition, the comparison of T2 with T3 and T1 with T4 or T5 would all have been significant. (But we would have overstated things).

How about as an F-test? For Scheffe’s method we need,

S^{2} = (r − 1)F [1 − α; r − 1; ν]

. If we are only going to do the pairwise differences, then ŝ(D̂ is the same as before. Now

> qf(.90, 4, 8)
[1] 2.806426

so, S = \sqrt{(r − 1)F [1 − α; r − 1; ν]} is

> sqrt(4*qf(.90, 4, 8))
[1] 3.350478

Now the multiplier is 3.35, so the confidence interval is

0.6 ± 3.35 × 0.17 = 0.6 ± 0.57

So with this method we only get a significant difference between T2 and T5.

An analysis in SAS

Here is a SAS program to do the same thing we did in Splus.

options linesize=78;
data cox;
        input yield tmt $ block $ ;
        cards;
7.62     T1     B1
8.14     T2     B1
7.76     T3     B1
7.17     T4     B1
7.46     T5     B1
8.00     T1     B2
8.15     T2     B2
7.73     T3     B2
7.57     T4     B2
7.68     T5     B2
7.93     T1     B3
7.87     T2     B3
7.74     T3     B3
7.80     T4     B3
7.21     T5     B3
;
proc anova;
        class tmt block;
        model yield = tmt block;
        means tmt/bon tukey scheffe alpha=0.10;

(I will explain what the various lines mean in class—you should annotate this). The log file from that SAS run was

......
NOTE: AUTOEXEC processing beginning; file is /usr/local/lib/sas/autoexec.sas.
NOTE: SAS initialization used:
      real time           0.19 seconds
      cpu time            0.18 seconds

NOTE: AUTOEXEC processing completed.
1          options linesize=78;
2          data cox;
3               input yield tmt $ block $ ;
4               cards;
NOTE: The data set WORK.COX has 15 observations and 3 variables.
NOTE: DATA statement used:
      real time           0.28 seconds
      cpu time            0.06 seconds

20         ;
21         proc anova;
22              class tmt block;
23              model yield = tmt block;
24              means tmt/bon tukey scheffe alpha=0.10;
NOTE: The PROCEDURE ANOVA printed pages 1-5.
NOTE: PROCEDURE ANOVA used:
      real time           0.54 seconds
      cpu time            0.09 seconds

NOTE: The SAS System used:
      real time           1.06 seconds
      cpu time            0.37 seconds
.....

and the susbstantive output follows. Notice that the significant differences we found are confirmed with SAS.

                                The SAS System                               1
                                             09:19 Wednesday, January 31, 1996
                        Analysis of Variance Procedure
                           Class Level Information

                      Class    Levels    Values
                      TMT           5    T1 T2 T3 T4 T5
                      BLOCK         3    B1 B2 B3
                   Number of observations in data set = 15
                                The SAS System                               2
                                             09:19 Wednesday, January 31, 1996
                        Analysis of Variance Procedure
Dependent Variable: YIELD
                                    Sum of           Mean
Source                  DF         Squares         Square   F Value     Pr > F
Model                    6      0.82956000     0.13826000      3.16     0.0677
Error                    8      0.34948000     0.04368500
Corrected Total         14      1.17904000
                  R-Square            C.V.       Root MSE           YIELD Mean
                  0.703589        2.706677       0.209010             7.722000
Source                  DF        Anova SS    Mean Square   F Value     Pr > F
TMT                      4      0.73244000     0.18311000      4.19     0.0404
BLOCK                    2      0.09712000     0.04856000      1.11     0.3750
                                The SAS System                               3
                                             09:19 Wednesday, January 31, 1996
                        Analysis of Variance Procedure
           Tukey's Studentized Range (HSD) Test for variable: YIELD
         NOTE: This test controls the type I experimentwise error rate, but
               generally has a higher type II error rate than REGWQ.

 Alpha= 0.1  df= 8  MSE= 0.043685
Critical Value of Studentized Range= 4.169
Minimum Significant Difference= 0.503
Means with the same letter are not significantly different.
Tukey Grouping              Mean      N  TMT
        A            8.0533      3  T2
        A
B       A            7.8500      3  T1
B       A
B       A            7.7433      3  T3
B
B                    7.5133      3  T4
B
B                    7.4500      3  T5

          The SAS System                               4
                       09:19 Wednesday, January 31, 1996
  Analysis of Variance Procedure
Bonferroni (Dunn) T tests for variable: YIELD
NOTE: This test controls the type I experimentwise error rate, but
generally has a higher type II error rate than REGWQ.
 Alpha= 0.1  df= 8  MSE= 0.043685
    Critical Value of T= 3.36
Minimum Significant Difference= 0.5726
Means with the same letter are not significantly different.
Bon Grouping              Mean      N  TMT
       A            8.0533      3  T2
       A
B       A            7.8500      3  T1
B       A
B       A            7.7433      3  T3
B       A
B       A            7.5133      3  T4
B

B                    7.4500      3  T5

         The SAS System                               5
                      09:19 Wednesday, January 31, 1996
 Analysis of Variance Procedure
Scheffe's test for variable: YIELD
NOTE: This test controls the type I experimentwise error rate but
generally has a higher type II error rate than REGWF for all
pairwise comparisons
Alpha= 0.1  df= 8  MSE= 0.043685
  Critical Value of F= 2.80643
Minimum Significant Difference= 0.5718
Means with the same letter are not significantly different.
Scheffe Grouping              Mean      N  TMT
        A            8.0533      3  T2
        A
B       A            7.8500      3  T1
B       A
B       A            7.7433      3  T3
B       A
B       A            7.5133      3  T4
B
B                    7.4500      3  T5

Randomized Block Experiments

The following data are drawn from Box, Hunter, and Hunter, Statistics for Experimenters. 4 different methods (treatments) are used in the manufacture of penicillin. 5 different blends of the raw material are available for the process. (We want to use different blends to ensure that that whatever method we use is good over a wide variety of situations). The data are

block  treatment
A B C D
blend 1 89^{1} 88^{3} 97^{2} 94^{4}
blend 2 84^{4} 77^{2} 92^{3} 79^{1}
blend 3 81^{2} 87^{1} 87^{4} 85^{3}
blend 4 87^{1} 92^{3} 89^{2} 84^{2}
blend 5 79^{3} 81^{4} 80^{1} 88^{2}

The superscripts represent the order in which the experiments were performed. We are interested in the treatment effects—we don’t really care about the block (blend) effects.

We need to consider both models for the data, and ways of looking at the data. First, lets think about how to enter the data. I re-expressed the data in the following format

Yld Blend Tmt Order
89 blend1  A    1
88 blend1  B    3
97 blend1  C    2
94 blend1  D    4
.....

and then tried to do some simple plots within S,

> corn <- read.table("bhh-corn.dat", header=T)
> corn
   Yld  Blend Tmt Order
 1  89 blend1   A     1
 2  88 blend1   B     3
 3  97 blend1   C     2
 4  94 blend1   D     4
 . . .

> plot.factor(corn)
> interaction.plot(Blend, Tmt, Yld)
> interaction.plot(Tmt, Blend, Yld)

These plots appear below.

What is an appropriate model for these data?

How about

y_{ij} = μ + β_{i} + τ_{j} + ϵ_{ij}

where i identifies the blend and j identifies the tmt.

Our next step is to set up an ANOVA. Here is a hint of what we might do

> attach(corn, pos=1)
> ls()
[1] "Blend" "Order" "Tmt"   "Yld"
> mymod <- aov(Yld ~ Blend + Tmt)
> anova(mymod)
Analysis of Variance Table
Response: Yld
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
Blend      4       264 66.00000 3.504425 0.0407462
Tmt        3        70 23.33333 1.238938 0.3386581
Residuals 12       226 18.83333

and

> anova( aov(Yld ~ Blend + Tmt + Order) )
Analysis of Variance Table
Response: Yld
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
Blend      4  264.0000 66.00000 3.531415 0.0435472
Tmt        3   70.0000 23.33333 1.248480 0.3392314
Order      1   20.4167 20.41667 1.092420 0.3183604
Residuals 11  205.5833 18.68939

Things not to forget

Remember that we are doing linear analyses just like regression, so many of the “tricks” you learnt, like transformations of the response variable and diagnostics about influence, are appropriate here. The careful data analyst will take all these into account with real data.

More examples

The following are some gedanken examples and comments about multifactor experiments, followed by some data.

Example 1 A company investigated the effect of selling price and type of promotion campaign on sales of one if its products. Three selling prices (59, 60, and 64 cents) were studies, as were two types of promotional campaigns (radio, and newspaper). Consider the prices to be factor A and the campaign to be factor B. The experiment was done using each combination of A and B (that is there were 6 runs of the experiment). That isn’t strictly necessary, but does make the analysis a little easier. Each combination was tried in 2 communities, for a total of 12 experimental runs.

This is an experiment because the experimenters exercised control of the assignment of factors (A and B) to treatment units. This is an example of a completely randomized design, with two factors.

Example 2 A steel company studied the effects of carbon content and tempering temperature on the strength of steel. Carbon content was investigated at a high and low level. Tempering temperature was also studied at two levels, high and low. Thus we have a 2 × 2 design with 4 possible treatment combinations.

The 4 treatments were assigned to 12 production batches in a randomized fashion. This is again an experiment, and again a completely randomized design.

Example 3 A researcher studied the effects of family income (under $10,000, 10,000–29,999, 30,000–49,999, over $50,000) and stage in life family life cycle (4 stages from just married to just retired) on appliance purchases. Here we have 4 × 4 = 16 possible treatments.

The researcher selected 20 families with the required income and life-cycle characteristics for each treatment combination, yielding 320 families.

This is an observational study—the data were collected without assigning the factor levels to experimental unit. This data would generally be analyzed as though it were a completely randomized design. (Lots of things can go wrong. One obvious problem is that the chances of being in a family go down as one gets older, …)

If one only has one observation per treatment combination, then one can extract “main effects” for the factors and “interactions”, but there is nothing left to estimate random variability (error). With multiple observations per cell, one can estimate main effects and interactions an still have something left with which to estimate error.

One does not have to do the complete replication of the experiment (all factors at all levels). For example, suppose we wished to investigate an industrial process that had three factors with 6, 5, and 4 levels. That would yield 6 × 5 × 4 = 120 combinations. It is probably not possible to do 120 runs of the process. However, by carefully selecting which combinations to run, one can still recover the main effects, while sacrificing the ability to estimate some interactions.

We will mainly focus on completely replicated experiments, unless someone wants to see something else.

Some more data

Here is some more data. Again, data on blood pressure. There are two types of treatment to consider, Biofeedback (Yes or No) and Drugs (three types, A, B, and C). An experiment is performed and there are 5 replications per treatment combination. The raw data is

Drug A Drug B Drug C
Bio Y 170  186  180
175  194  187
165  201  199
180  215  170
160  219  204
Bio N 173  189  202
194  194  228
197  217  190
190  206  206
176  199  224

The easiest way to lay out this data is in the following form

Resp    Biofeedback    Drug
170        Y            A
175        Y            A
165        Y            A
180        Y            A
160        Y            A
186        Y            B
194        Y            B
201        Y            B
215        Y            B
219        Y            B
180        Y            C
187        Y            C
199        Y            C
170        Y            C
204        Y            C
173        N            A
194        N            A
197        N            A
......

There are many questions to ask of the data, starting with

• What are the main effects of the two treatments, and do they matter?

• Are there any interactions between the treatments?

A good place to start is with an interaction plot of the full data

> blood2 <- read.table("blood2.dat", header=T)
> blood2
   Resp Biofeedback Drug
 1  170           Y    A
 . . . .
> plot.factor(blood2)
> attach(blood2, pos=1)
> ls()
[1] "Biofeedback" "Drug"        "Resp"
> interaction.plot(Drug, Biofeedback, Resp)
> interaction.plot(Biofeedback, Drug, Resp)

The interaction plots are very interesting. It certainly looks like there might be an interaction between the two factors. Lets try a very simple analysis to confirm/deny this.

> ls()
[1] "Biofeedback" "Drug"        "Resp"
> mymod <- aov( Resp ~ Biofeedback * Drug )
> anova(mymod)
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
                 Df Sum of Sq Mean Sq  F Value     Pr(F)

Biofeedback       1      1080 1080.00  6.93419 0.0145635
Drug              2      3420 1710.00 10.97913 0.0004113
Biofeedback:Drug  2       780  390.00  2.50401 0.1028767
Residuals        24      3738  155.75
>

The ANOVA table suggests that both the Bio and Drug treatments are effective (why?) and that there is borderline evidence for the interaction(why?).

What to do now?

• Do some more analyses in class

• Wonder about what component of the interaction is important?

The Plan

• What does interaction mean?

• GLIM notation for main-effects, interactions (Venables and Ripley, p. 153)

• Tests for Main Effects and Interactions

• Contrasts

• Unequal sample sizes in 2-way ANOVA

• Multifactor Studies

• Analysis of Covariance

• Some Standard Designs (Latin Square and Split Plot). See also Venables and Ripley p. 163 et. seq. for a discussion of ANOVA models.

Things to talk about

We have many options about where to go next, but the two likely areas are

• Categorical Data Analysis, including logistic regression, and the so-called generalized linear models. This covers the analysis of data with binary responses, or counts as the response variable. If we do this, it will take between 3 and 5 weeks.

• Robust statistical methods and non-linear regression. This depends on your interests. It should take a couple of weeks.

• Modern regression methods, including additive models and projection-pursuit. De-

pending on how carefully we do this, it can take from 2 to 5 weeks. Which do you prefer?

The Plan

• What does interaction mean?

• GLIM notation for main-effects, interactions

• Tests for Main Effects and Interactions

• Contrasts

• Unequal sample sizes in 2-way ANOVA

• Multifactor Studies

• Analysis of Covariance

• Some Standard Designs (Latin Square and Split Plot).

Interaction

The concept of interaction is both simple and tricky at the same time. Consider a two-way layout of data, indexed by rows i and columns j. Lets call the data Y_{i,j}. A natural additive model for the data is

Y_{i,j} = μ + α_{i} + β_{j} + ϵ_{i,j}

where the errors have some (usually normal) distribution with mean zero and constant variance. Lets assume that μ = 0 and assign arbitrary row and column effects.

14 0 16 10
24 20 26 20
-26 -40 -24 -30
4 -10 6

This data is perfectly additive—the entries in the table are exactly the row effect plus the column effect. Interaction is the term that describes the departure from the additive model. When we only have one entry per cell, then we usually assign interaction to the error (ϵ_{i,j} above, but when there is more than one observation per cell we can get an estimate of interaction (how does the cell mean differ from the pure additive model) and the pure error (what is the variability of the observations around the cell mean).

Even with only one observation per cell, we might sometimes choose to dissect the error into an interaction component and an error component. This is most commonly done when there is a specific type of interaction that we are looking for.

For example, consider the following numbers (I hesitate to call them data as I made them up).

> jj
   V1 V2 V3
1  18  A  a
2  -6  A  b
3  20  A  c
4  29  B  a
5  16  B  b
6  13  B  c
7 -28  C  a
8 -36  C  b
9 -17  C  c
> attach(jj,pos=1)
> mymod <- aov(V1 ~ V2 + V3 )

So we have some numbers in a 3 × 3 layout, and we fit the simple model with row and column effects. The ANOVA table is

Analysis of Variance Table
Response: V1
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
V2         2  3640.667 1820.333 22.51959 0.0066532
V3         2   422.000  211.000  2.61031 0.1881914
Residuals  4   323.333   80.833

Plowing ahead, the researcher believes there is cause to worry about a linear times linear interaction—an effect which is linear in rows and columns. To test this theory she creates a (one degree of freedom) effect to add to the model. The particular effect is

> as.numeric(V2)
[1] 1 1 1 2 2 2 3 3 3
> LinXlin <- as.numeric(V2) * as.numeric(V3)

> LinXlin
[1] 1 2 3 2 4 6 3 6 9

which is easily seen to be linear in both components. The resulting ANOVA and summary is

> mymod <- aov(V1 ~ V2 + V3 + LinXlin)
> anova(mymod)
Analysis of Variance Table
Response: V1
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
V2         2  3640.667 1820.333 18.01815 0.0213049
V3         2   422.000  211.000  2.08853 0.2702472
LinXlin    1    20.250   20.250  0.20044 0.6847201
Residuals  3   303.083  101.028

Notice, all we did was add the new effect to the model and the SS for that effect was extracted from the SS for residuals and assigned to the named effect. In these data there doesn’t seem to be much of a linear times linear effect. The model we fit was

Y_{i,j} = μ + α_{i} + β_{j} + τ × LL + ϵ_{i,j}

where LL is the linear times linear vector. Note that LL can not be written as a linear combination of two vectors.

How else might we break apart the SS for residuals? In the first analysis there were 4 degrees of freedom for error. There are MANY ways we can partition those 4 degrees of freedom, but a common method is into linear and quadratic components.

Below, we create linear and quadratic components of the interaction.

> LL <- as.numeric(V2) * as.numeric(V3)
> LQ <- as.numeric(V2) * (as.numeric(V3))^2
> QQ <- (as.numeric(V2))^2 * (as.numeric(V3))^2
> QL <- (as.numeric(V2))^2 * (as.numeric(V3))
> data.frame(LL,LQ,QL,QQ)
  LL LQ QL QQ
1  1  1  1  1
2  2  4  2  4
3  3  9  3  9
4  2  2  4  4
5  4  8  8 16
6  6 18 12 36
7  3  3  9  9
8  6 12 18 36
9  9 27 27 81

If we now do an ANOVA with all these effects we should find ourselves with no degrees of freedom for error, as we do.

> mymod <- aov(V1 ~ V2 + V3 + LL + LQ + QL + QQ)
> anova(mymod)
Analysis of Variance Table
Response: V1
Terms added sequentially (first to last)
   Df Sum of Sq  Mean Sq
V2  2  3640.667 1820.333
V3  2   422.000  211.000
LL  1    20.250   20.250
LQ  1    44.083   44.083
QL  1   168.750  168.750
QQ  1    90.250   90.250

Just as a check, lets fit the model with no QQ term (that is the model that allows for all interactions except the quadratic by quadratic one). It looks like

> anova(mymod <- aov(V1 ~ V2 + V3 + LL + LQ + QL))
Analysis of Variance Table
Response: V1
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
V2         2  3640.667 1820.333 20.16990 0.1555306
V3         2   422.000  211.000  2.33795 0.4197422
LL         1    20.250   20.250  0.22438 0.7183758
LQ         1    44.083   44.083  0.48846 0.6116710
QL         1   168.750  168.750  1.86981 0.4019821
Residuals  1    90.250   90.250

Get the general idea? Interactions can be your friend, if you know how to manipulate them.

Notation for Effects and Interactions

There are many different notations for main effects and interactions. Many books just use model formulae and subscripts. The S language, and most other computer packages, use some variant of what has come to be known as GLIM notation. (GLIM was an early statistical package for Generalized Linear Interactive Modelling. GLIM is still available, but largely superseded by programs like S). The general idea is that there is a model “calculus” for the formulae. We have already seen some examples.

Consider Y a response vector and R is a factor with n_{R} levels, (Remember that a factor is just a variable whose entries are levels rather than numbers), C is a factor with n_{C} levels and L is a factor with n_{L} levels. Lets work our way through some examples, providing a narrative and symbolic interpretation of the model formulae.

• The formula Y ∼ R + C + L says that the response can be expressed as a linear combination of row levels, column levels, and layer levels. In math notation, this translates into

Y_{ijk} = μ + R_{i} + C_{j} + L_{k}

, plus the usual error. The expected response lives in the linear space of all Y_{ijk} that can be written in that form.

• The formula Y ∼ R + C + L + C : L says that we have main effects for rows, columns, and layers, but also allow for a column times layer interaction. In math that is,

Y_{ijk} = μ + R_{i} + C_{j} + L_{k} + CL_{jk}

• A shorthand form for the above model is Y ∼ R + C ∗ L where the C ∗ L explicitly expands into C + L + C : L

• One could write the model Y ∼ R + C : L where the C : L is meant to refer to the interaction without the corresponding main effects, but this rarely makes sense (and S won’t easily allow it).

• A final example is Y ∼ R ∗ C ∗ L which expands out to Y ∼ R + C + L + R : C + R : L + C : L + R : C : L, and in math notation,

Y_{ijk} = μ + R_{i} + C_{j} + L_{k} + CL_{jk} + RC_{ij} + RL_{ik} + RCL_{ijk}

Tests for Main Effects and Interactions

The best way to look at testing for effects and interactions is to do an example. Here we steal an example from Neter, Wasserman, Kunter.

Two factor study. Bakery sales in display cases

Display Height i  Display Width j Row Display Height
regular wide Total  Average
Bottom  47  46
43  40 176  44
Middle  62  67
68  72 268  67
Top  41  42
39  46 168  42
Column Total  300  312 612
Display width average  50  52  51

Here is the data in a form that I can use

Resp  Width Height
47      1       1
46      2       1
43      1       1
40      2       1
62      1       2
67      2       2
68      1       2
71      2       2
41      1       3
42      2       3
39      1       3
46      2       3

Now lets get down to some analyses.

> castle <- read.table("castle.dat", header=T)
> castle <- data.frame(Resp=castle$Resp, Height=as.factor(castle$Height),
+ Width=as.factor(castle$Width) )
> castle
   Resp Height Width
 1   47      1     1
 2   46      1     2
 3   43      1     1
 4   40      1     2
 5   62      2     1
 6   67      2     2
 7   68      2     1
 8   71      2     2

 9   41      3     1
10   42      3     2
11   39      3     1
12   46      3     2
>

That just made our factors, be factors!

> attach(castle, pos=1)
> par(mfrow=c(2,2))
> plot.factor(castle)
> interaction.plot(Width, Height, Resp)
> interaction.plot(Height, Width, Resp)

And here are the plots

Now onto a quick analysis.

> ls()
[1] "Height" "Resp"   "Width"
> anova(mymod <- aov(Resp ~ Height + Width + Height:Width) )
Analysis of Variance Table

Response: Resp
Terms added sequentially (first to last)
             Df Sum of Sq  Mean Sq  F Value     Pr(F)
Height        2      1544 772.0000 74.70968 0.0000575
Width         1        12  12.0000  1.16129 0.3226055
Height:Width  2        24  12.0000  1.16129 0.3746966
Residuals     6        62  10.3333

It certainly does not look like the interaction is significant, so it makes sense to drop it from the model

> anova(mymod <- aov(Resp ~ Height + Width) )
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
          Df Sum of Sq Mean Sq  F Value     Pr(F)
Height     2      1544  772.00 71.81395 0.0000077
Width      1        12   12.00  1.11628 0.3215794
Residuals  8        86   10.75

We still don’t have a significant effect for Width (how would we have known that?), so lets drop that term as well.

> anova(mymod <- aov(Resp ~ Height) )
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value        Pr(F)
Height     2      1544 772.0000 70.89796 3.099832e-06
Residuals  9        98  10.8889
>

And this looks like a reasonable model that is consistent with the plots.

Bunches of Examples

The following data are the survival times of groups of four animals randomly allocated to three poisons and four treatments. Then experiment was part of an investigation to combat the effects of certain toxic agents. The is a 4-fold replication of a 3 × 4 factorial design. Do an analysis of these data.

Treatment

poison A B C D
I  0.31 0.82 0.43 0.45
0.45 1.10 0.45 0.71
0.46 0.88 0.63 0.66
0.43 0.72 0.76 0.62
II  0.36 0.92 0.44 0.56
0.29 0.61 0.35 1.02
0.40 0.49 0.31 0.71
0.23 1.24 0.40 0.38
III 0.22 0.30 0.23 0.30
0.21 0.37 0.25 0.36
0.18 0.38 0.24 0.31
0.23 0.29 0.22 0.33

As always the first step is to enter the data. This time, I’ll just enter the data and create the factor levels on the fly.

Here is the raw data, in a file called poison.dat

0.31 0.82 0.43 0.45
0.45 1.10 0.45 0.71
0.46 0.88 0.63 0.66
0.43 0.72 0.76 0.62
0.36 0.92 0.44 0.56
0.29 0.61 0.35 1.02
0.40 0.49 0.31 0.71
0.23 1.24 0.40 0.38
0.22 0.30 0.23 0.30
0.21 0.37 0.25 0.36
0.18 0.38 0.24 0.31
0.23 0.29 0.22 0.33

and here is my first pass at reading in the data and creating appropriate factors (note

the use of as.factor.
> poison <- scan("poison.dat")

> poison
 [1] 0.31 0.82 0.43 0.45 0.45 1.10 0.45 0.71 0.46 0.88 0.63 0.66 0.43 0.72 0.76
[16] 0.62 0.36 0.92 0.44 0.56 0.29 0.61 0.35 1.02 0.40 0.49 0.31 0.71 0.23 1.24
[31] 0.40 0.38 0.22 0.30 0.23 0.30 0.21 0.37 0.25 0.36 0.18 0.38 0.24 0.31 0.23
[46] 0.29 0.22 0.33
> length(poison)
[1] 48

Create a vector of the Tmt levels,

> rep(1:4, len=48)
 [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2
[39] 3 4 1 2 3 4 1 2 3 4
> Tmt <- as.factor( rep(1:4, len=48 ) )

and one of the Poison levels

> rep(1:3, rep(16,3))
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3
[39] 3 3 3 3 3 3 3 3 3 3
> Poison <- as.factor(rep(1:3, rep(16,3)))

and finally a data frame to hold all the data neatly together. (You might not want to do this, but my data area is always a mess).

> poison.frame <- data.frame(poison, Tmt, Poison)
> poison.frame
   poison Tmt Poison
 1   0.31   1      1
 2   0.82   2      1
 3   0.43   3      1
 4   0.45   4      1
> rm(poison, Tmt, Poison)

So what sort of model might make sense here? How about

Y_{ij} = μ + τ_{i} + β_{j} + ω_{ij} + ϵ_{ijk}

where k implicitly indexes the replications?

That model allows us to estimate main effects and interactions, and still leaves room for random error? First thing to do is look at some pictures.

> plot.factor(poison.frame)
> attach(poison.frame, pos=1)
> interaction.plot(Poison, Tmt, poison)
> interaction.plot(Tmt, Poison, poison)

The plot of the RAW DATA is

It looks like there are treatment differences and some potential interactions. Why not let the ANOVA do its job and see where we go.

> anova(mymod <- aov(poison ~ Poison + Tmt))
Analysis of Variance Table
Response: poison
Terms added sequentially (first to last)
          Df Sum of Sq   Mean Sq  F Value        Pr(F)
Poison     2  1.033012 0.5165062 20.64329 5.703730e-07
Tmt        3  0.921206 0.3070688 12.27267 6.696971e-06
Residuals 42  1.050862 0.0250205

That analysis looks OK, except we forgot to add the interaction term. Here is a version with interactions.

> anova(mymod <- aov(poison ~ Poison + Tmt + Poison*Tmt))
Analysis of Variance Table

Response: poison
Terms added sequentially (first to last)
           Df Sum of Sq   Mean Sq  F Value     Pr(F)
Poison      2  1.033012 0.5165062 23.22174 0.0000003
Tmt         3  0.921206 0.3070688 13.80558 0.0000038
Poison:Tmt  6  0.250138 0.0416896  1.87433 0.1122506
Residuals  36  0.800725 0.0222424

It appears there are strong main effects, but little or no interaction.

At the very least, we should now try a residual plot of these data.

Those are not stellar looking plots. The Box plots do not look uniform, and there seems to be a relationship between residual and fitted value. What to do?

One could go through a formal process of selecting a transformation, perhaps from the Box-Cox family. I’ve decided to just use an inverse, and sqrt.

> anova(mymod <- aov( (poison)^(-1) ~ Poison + Tmt + Poison*Tmt))
Analysis of Variance Table
Response: (poison)^(-1)

Terms added sequentially (first to last)
           Df Sum of Sq  Mean Sq  F Value     Pr(F)
Poison      2  34.87712 17.43856 72.63475 0.0000000
Tmt         3  20.41429  6.80476 28.34307 0.0000000
Poison:Tmt  6   1.57077  0.26180  1.09042 0.3867329
Residuals  36   8.64308  0.24009
> plot(1:10)
> plot.factor(data.frame(resid(mymod), Poison, Tmt))
> plot(fitted(mymod),resid(mymod))

> anova(mymod <- aov( (poison)^(0.5) ~ Poison + Tmt + Poison*Tmt))
Analysis of Variance Table
Response: (poison)^(0.5)
Terms added sequentially (first to last)
           Df Sum of Sq   Mean Sq  F Value     Pr(F)
Poison      2 0.5609264 0.2804632 34.38932 0.0000000
Tmt         3 0.4306422 0.1435474 17.60123 0.0000003

Poison:Tmt  6 0.0793987 0.0132331  1.62259 0.1692592
Residuals  36 0.2935992 0.0081555
> plot.factor(data.frame(resid(mymod), Poison, Tmt))
> plot(fitted(mymod),resid(mymod))

Note that the inverse transformation seems to be the ticket. Note particularly that any hint of interaction is completely missing on the inverse scale.

A Latin Square Design

Suppose that 4 cars and four drivers are employed in a study of possible differences between four gasoline additives. Even though the cars are identical models, slight systematic differences are likely to occur in their performance, and even though each driver may do her best to drive the car using the experimental protocol, slight systematic differences can occur from driver to driver. It would be desirable to design (and analyze) an experiment that eliminates both the driver-to-drive and automobile-to-automobile effects. The Latin Square arrangement below allows this to be done.

cars

1 2 3 4
I A B D C
drivers II D C A B
III B D C A
IV C A B D

with cars (I,II,III, and IV), drivers (1, 2, 3, and 4) and additive (A,B,C,and D).

Data will be assigned to each combination, so that data will be indexed by three things (car, driver, and additive), each with 4 levels, but there will only be 16 data points. We do not have a complete replication of the experiment.

The model might be written

Y_{ijk} = μ + d_{i} + c_{j} + a_{k} + ϵ_{ijk}

There is not much hope of estimating interactions.

Here is some data from this experiment.

cars

1 2 3 4
I A 21 B 26 D 20 C 25
drivers II D 23 C 26 A 20 B 27
III B 15 D 13 C 16 A 16
IV C 17 A 15 B 20 D 20

Lets do a quick analysis of these data, first entering them, and then an ANOVA and plots. Here is gas.dat

A 21
B 26
D 20
C 25
D 23
C 26
A 20
B 27
B 15
D 13
C 16
A 16
C 17
A 15
B 20
D 20

and here is the data being entered.

> delme <-read.table("gas.dat")
> delme

   V1 V2
 1  A 21
 2  B 26
 3  D 20
 4  C 25
 5  D 23
 6  C 26
 7  A 20
 8  B 27
 9  B 15
10  D 13
11  C 16
12  A 16
13  C 17
14  A 15
15  B 20
16  D 20
> gas <- delme$V1
> resp <- delme$V2
> car <- rep(1:4, len=16)
> car
 [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
> driver <- (1:4,rep(4,4))
Syntax error: "," used illegally at this point:
driver <- (1:4,
> driver <- rep(1:4,rep(4,4))
> driver
 [1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
> gas.frame <- data.frame(resp, gas, car=as.factor(car), driver=as.factor(driver) )

So here is my constructed data set

> gas.frame
   resp gas car driver
 1   21   A   1      1
 2   26   B   2      1
 3   20   D   3      1
 4   25   C   4      1
 . . . . . .

Now onto an analysis.

> anova(mymod <- aov(resp ~ gas + car + driver) )
Analysis of Variance Table
Response: resp

Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq F Value     Pr(F)
gas        3        40 13.33333       5 0.0451975
car        3        24  8.00000       3 0.1169598
driver     3       216 72.00000      27 0.0006987
Residuals  6        16  2.66667

Whoops, that was constructed data, look at those residuals!!!

Yet more design examples

In this class we will consider several examples of designs, and essentially finish our discussion of ANOVA. Some of the calculations I’ll do (particularly figuring out alias patterns in S) take much more computer power than they should.

The material is copied from Chapter 13 of Box, Hunter, and Hunter, Statistics for Experimenters.

Example 13.1

Here we are asked to figure out what the alias pattern is. My version of this follows. First we enter the data (this is by now a common thing for use to do). I added a Response because S seemed to want that in the the subsequent calculations.

> rr
  V1 V2 V3 V4 V5 Resp
1  0  0  0  1  0    1
2  1  0  0  1  1    2
3  0  1  0  0  1    3
4  1  1  0  0  0    4
5  0  0  1  0  1    5
6  1  0  1  0  0    6
7  0  1  1  1  0    7
8  1  1  1  1  1    8

Then attach the dataset and ask for the alias pattern

> attach(rr)
> alias(Resp ~ (V1+V2+V3+V4+V5)^5 )
Model
Resp ~ (V1 + V2 + V3 + V4 + V5)^5
Complete
               (Intercept) V1 V2 V3 V4 V5 V1:V2 V1:V3
         V1:V4                          1
         V1:V5                       1
         V2:V3                       1
         V2:V4                    1
         V2:V5                                   1
         V3:V4                 1
         V3:V5                             1
         V4:V5              1
      V1:V2:V3                          1
      V1:V2:V4                                   1
      V1:V2:V5                    1
      V1:V3:V4                             1

      V1:V3:V5                 1
      V1:V4:V5  1
      V2:V3:V4  1
      V2:V3:V5              1
      V2:V4:V5                             1
      V3:V4:V5                                   1
   V1:V2:V3:V4              1
   V1:V2:V3:V5  1
   V1:V2:V4:V5                 1
   V1:V3:V4:V5                    1
   V2:V3:V4:V5                          1
V1:V2:V3:V4:V5                       1

Lets follow through on this. We have a 1/4 fraction of the full design so we expect (if the design was well done) that each estimable effect will be aliased with three other effects. The output tells us that the estimable effects (the first row) are the main effects and a couple of two-way interactions. All the other two-way, three-way and higher interactions are confounded/aliased with the estimable effects. So, V1 is the same as V4:V5 V2:V3:V5 and V1:V2:V3:V4. This is a special kind of design, consructed with considerable care.

Example 13.2

It should be clear from the text that one can “analyse” these data without much work, but some supporting calculations are as follows.

> e2
  V1 V2 V3 V4 Resp
1  0  0  0  0   20
2  1  0  0  1   14
3  0  1  0  1   17
4  1  1  0  0   10
5  0  0  1  1   19
6  1  0  1  0   13
7  0  1  1  0   14
8  1  1  1  1   10
> anova(mymod <- aov(Resp ~ V1 + V2 + V3 + V4, data=e2))
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
V1         1    66.125 66.12500 144.2727 0.0012415
V2         1    28.125 28.12500  61.3636 0.0043321
V3         1     3.125  3.12500   6.8182 0.0796050

V4         1     1.125  1.12500   2.4545 0.2151699
Residuals  3     1.375  0.45833
> anova(mymod <- aov(Resp ~ (V1 + V2 + V3 + V4)^2, data=e2))
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
      Df Sum of Sq Mean Sq
V1     1    66.125  66.125
V2     1    28.125  28.125
V3     1     3.125   3.125
V4     1     1.125   1.125
V1:V2  1     0.125   0.125
V1:V3  1     1.125   1.125
V1:V4  1     0.125   0.125

To treat it as a replicated design on just the first two variables, the appropriate ANOVA is

> anova(mymod <- aov(Resp ~ V1 + V2, data=e2))
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
          Df Sum of Sq Mean Sq  F Value       Pr(F)
V1         1    66.125  66.125 58.77778 0.000601443
V2         1    28.125  28.125 25.00000 0.004104716
Residuals  5     5.625   1.125

Factorial Designs—Examples

Factorial designs with each factor at two levels are extremely important in practice. Many experiments, especially industrial ones, have factors at 2 levels. Things like (off,on), (high, low), (treatment, placebo).

To perform a general factorial design, select a fixed number of levels for each factor in the design. If the levels are l_{1}, l_{2}, …, l_{k}, then such an experiment is called an l_{1} × l_{2} × ⋯ × l_{k} factorial design. When l_{k} = 2 for all k then it is called a 2^{k} factorial design.

Factorial designs, and in particular 2^{k} factorial design are important because:

• They require relatively few runs per factor studied.

• It is possible to run fractional designs that are even more economical—at some expense in accuracy/estimability.

• They are easy to interpret, and some simple effects can be seen with common sense and a little arithmetic. (Hence, it is easier to convince your boss of what is happening)

There are many standard notations for a (say) 2^{3} design. If the three factors are T, C, and K, each at 2 levels, then here are just some of the ways such a design is represented.

Run T C K  T C K
1 − − − 1 0 0 0
2 + − − t 1 0 0
3 − + − c 0 1 0
4 + + − tc 1 1 0
5 − − + k 0 0 1
6 + − + tk 1 0 1
7 − + + ck 0 1 1
8 + + + tck 1 1 1

Here is data from a particular experiment, coded in a couple of different ways. The original data were part of a pilot plant investigation of a process and have been simplified for this example. There are two quantitative predictors (temperature and concentration) and a qualitative predictor (catalyst). The response is the chemical yield.

test  temp concentration catalyst yield
condition (C)  (%)  (A or B) (grams)
number  T  C  K  y

Original units of variables

1  160  20  A  60
2  180  20  A  72
3  160  40  A  54
4  180  40  A  68
5  160  20  B  52
6  180  20  B  83
7  160  40  B  45
8  180  40  B  80

Coded units of variables

1  −  −  −  60
2  +  −  −  72
3  −  +  −  54
4  +  +  −  68
5  −  −  +  52
6  +  −  +  83
7  −  +  +  45
8  +  +  +  80

with

Temperature Concentration Catalyst
−  +  −  +  − +
160 180 20  40  A B

Lets try a quick analysis of these data, in the first form. Here is the data

Run Temp  Conc Caty Yield
1   160   20   A   60
2   180   20   A   72
3   160   40   A   54
4   180   40   A   68
5   160   20   B   52
6   180   20   B   83
7   160   40   B   45
8   180   40   B   80

and reading the data into S

> chem <- read.table("chemical.dat", header=T)
> chem
  Run Temp Conc Caty Yield
1   1  160   20    A    60
2   2  180   20    A    72
3   3  160   40    A    54
4   4  180   40    A    68
5   5  160   20    B    52
6   6  180   20    B    83
7   7  160   40    B    45
8   8  180   40    B    80
> mychem <- data.frame(Run=chem$Run, Temp=as.factor(chem$Temp),
+ Conc=as.factor(chem$Conc), Caty=chem$Caty, Yield=chem$Yield)
> chem <- mychem
> rm(mychem)
> chem
  Run Temp Conc Caty Yield
1   1  160   20    A    60
2   2  180   20    A    72
3   3  160   40    A    54
4   4  180   40    A    68
5   5  160   20    B    52
6   6  180   20    B    83
7   7  160   40    B    45
8   8  180   40    B    80

Note the “fancy” footwork to make sure we have factors.

> plot.factor(chem, y=Yield)
> interaction.plot(Temp, Conc, Yield)
> interaction.plot(Temp, Caty, Yield)
> interaction.plot(Conc, Caty, Yield)

51

And now the simple ANOVA.

> mymod <- aov(Yield ~ Temp + Conc + Caty + Temp*Caty, data=chem)
> anova(mymod)
Analysis of Variance Table
Response: Yield
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq F Value     Pr(F)
Temp       1    1058.0 1058.000   634.8 0.0001371
Conc       1      50.0   50.000    30.0 0.0119669
Caty       1       4.5    4.500     2.7 0.1988922
Temp:Caty  1     200.0  200.000   120.0 0.0016286
Residuals  3       5.0    1.667
> plot(mymod)

If one really wanted to look at all the two-way effects one might do

> mymod <- aov(Yield ~ Temp + Conc + Caty + Temp*Caty +
+ Temp*Conc + Caty*Conc, data=chem)
> anova(mymod)

Response: Yield
Terms added sequentially (first to last)
          Df Sum of Sq Mean Sq F Value     Pr(F)
Temp       1    1058.0  1058.0    2116 0.0138374
Conc       1      50.0    50.0     100 0.0634510
Caty       1       4.5     4.5       9 0.2048328
Temp:Caty  1     200.0   200.0     400 0.0318045
Temp:Conc  1       4.5     4.5       9 0.2048328
Caty:Conc  1       0.0     0.0       0 1.0000000
Residuals  1       0.5     0.5

One might also wish to do a QQ plot of the model effects

> mymod <- aov(Yield ~ Temp + Conc + Caty + Temp*Caty, data=chem)
> qqnorm(mymod)

Fractional Factorial Designs—Examples

Factorial designs are great, but the number of runs required for a full 2^{k} design increases exponentially in k. When k is large (say bigger than 4 or 5) it is often a good idea to run a fraction of the design—and recover a fraction of the possible information. Hopefully one designs the experiment to recover a useful fraction of the information.

Consider a two-level design with seven variables. A complete factorial design would require 2^{7} = 128 runs, which might be a prohibative number. These 128 degress of freedom would be partioned into main effects and interactions, as follows

average  main  2-fact  3-fact  4-fact  5-fact  6-fact  7-fact
  1       7      21     35      35      21      7       1

The high order interactions are probably not very interesting. If we were to drop them from the model (assume they are zero!), it would be very much like assuming that the higher order terms of a Taylor Series expansion are zero. (Actually, a Fourier series expansion is the better model, but you get the idea). The aim of a fractional factorial design is to run the experiment so that the high order interactions (which we might assume to be zero) are confounded with things that we would like to know.

The best way to proceed is via an example. First lets generate the design and add the data.

> reactor <- fac.design(c(2,2,2,2,2))
> attach(reactor, pos=1)
> Resp <- c(61,53,63,61,53,56,54,61,69,61,94,93,66,60,95,98,
+ 56,63,70,65,59,55,67,65,44,45,78,77,49,42,81,82)
> length(Resp)
[1] 32
> detach(1, save="newreact")
> reactor <- newreact
> rm(newreact)

Now lets just print out the data to look at.

> reactor
    A  B  C  D  E Resp
 1 A1 B1 C1 D1 E1   61
 2 A2 B1 C1 D1 E1   53
 3 A1 B2 C1 D1 E1   63
 4 A2 B2 C1 D1 E1   61
 5 A1 B1 C2 D1 E1   53
 6 A2 B1 C2 D1 E1   56
 7 A1 B2 C2 D1 E1   54
 8 A2 B2 C2 D1 E1   61
 9 A1 B1 C1 D2 E1   69

10 A2 B1 C1 D2 E1   61
11 A1 B2 C1 D2 E1   94
12 A2 B2 C1 D2 E1   93
13 A1 B1 C2 D2 E1   66
14 A2 B1 C2 D2 E1   60
15 A1 B2 C2 D2 E1   95
16 A2 B2 C2 D2 E1   98
17 A1 B1 C1 D1 E2   56
18 A2 B1 C1 D1 E2   63
19 A1 B2 C1 D1 E2   70
20 A2 B2 C1 D1 E2   65
21 A1 B1 C2 D1 E2   59
22 A2 B1 C2 D1 E2   55
23 A1 B2 C2 D1 E2   67
24 A2 B2 C2 D1 E2   65
25 A1 B1 C1 D2 E2   44
26 A2 B1 C1 D2 E2   45
27 A1 B2 C1 D2 E2   78
28 A2 B2 C1 D2 E2   77
29 A1 B1 C2 D2 E2   49
30 A2 B1 C2 D2 E2   42
31 A1 B2 C2 D2 E2   81
32 A2 B2 C2 D2 E2   82

Now suppose that we had not run the entire experiment, but only some of the runs, in particular, runs

> runs
 [1]  2  3  5  8  9 12 14 15 17 20 22 23 26 27 29 32

So our data might have been,

> half.reactor <- reactor[runs,]
> half.reactor
    A  B  C  D  E Resp
 2 A2 B1 C1 D1 E1   53
 3 A1 B2 C1 D1 E1   63
 5 A1 B1 C2 D1 E1   53
 8 A2 B2 C2 D1 E1   61
 9 A1 B1 C1 D2 E1   69
12 A2 B2 C1 D2 E1   93
14 A2 B1 C2 D2 E1   60
15 A1 B2 C2 D2 E1   95
17 A1 B1 C1 D1 E2   56
20 A2 B2 C1 D1 E2   65
22 A2 B1 C2 D1 E2   55

23 A1 B2 C2 D1 E2   67
26 A2 B1 C1 D2 E2   45
27 A1 B2 C1 D2 E2   78
29 A1 B1 C2 D2 E2   49
32 A2 B2 C2 D2 E2   82

The variables are feed rate (10, 15), catalyst (1%, 2%), agitation rate (100,120), temp (140, 180), concentration (3%, 6%) and the response is the percentage of the ingredients that reacted.

How would we analyse these data? First off we might like to know what effects are “aliased” with what other effects. Aliasing means that a contrast for one effect is exactly the same as a contract for another. Here is the alias pattern for these data.

> alias(half.reactor, ~ (A+B+C+D+E)^4)
Model
 ~ (A + B + C + D + E)^4
Complete
        (Intercept)  A  B  C  D  E A:B A:C A:D A:E B:C B:D B:E C:D C:E D:E
  A:B:C                                                                 1
  A:B:D                                                             1
  A:B:E                                                         1
  A:C:D                                                     1
  A:C:E                                                 1
  A:D:E                                             1
  B:C:D                                         1
  B:C:E                                     1
  B:D:E                                 1
  C:D:E                             1
A:B:C:D                          1
A:B:C:E                       1
A:B:D:E                    1
A:C:D:E                 1
B:C:D:E              1

It is important to note that the main effects and two-factor interactions are aliased with three- and four- factor interactions. This is a very desirable type of design (but how did we generate it?).

What about a quick (for now) analysis of these data?

> attach(half.reactor, pos=1)
> mymod <- aov(Resp ~ (A+B+C+D+E) + B*D + D*E )
> anova(mymod)
Analysis of Variance Table
Response: Resp

Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
A          1     16.00   16.000   2.3594 0.1630778
B          1   1681.00 1681.000 247.8894 0.0000003
C          1      0.00    0.000   0.0000 1.0000000
D          1    600.25  600.250  88.5161 0.0000134
E          1    156.25  156.250  23.0415 0.0013553
B:D        1    462.25  462.250  68.1659 0.0000348
D:E        1    361.00  361.000  53.2350 0.0000842
Residuals  8     54.25    6.781
>

What would this have looked like with the full data???

> anova(mymod)
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
A          1    15.125   15.125   1.3432 0.2578723
B          1  3042.000 3042.000 270.1499 0.0000000
C          1     3.125    3.125   0.2775 0.6031664
D          1   924.500  924.500  82.1018 0.0000000
E          1   312.500  312.500  27.7521 0.0000211
B:D        1  1404.500 1404.500 124.7290 0.0000000
D:E        1   968.000  968.000  85.9648 0.0000000
Residuals 24   270.250   11.260

The answer is only slightly different!

Yet more design examples

Hera are some more examples of designs, and essentially an end to our discussion of ANOVA. Some of the calculations I’ll do (particularly figuring out alias patterns in S) take much more computer power than they should.

The material is copied from Chapter 13 of Box, Hunter, and Hunter, Statistics for Experimenters.

Example 13.1

Here we are asked to figure out what the alias pattern is. My version of this follows. First we enter the data (this is by now a common thing for us to do). I added a Response because S seemed to want that in the the subsequent calculations.

> rr
  V1 V2 V3 V4 V5 Resp
1  0  0  0  1  0    1
2  1  0  0  1  1    2
3  0  1  0  0  1    3
4  1  1  0  0  0    4
5  0  0  1  0  1    5
6  1  0  1  0  0    6
7  0  1  1  1  0    7
8  1  1  1  1  1    8

Then attach the dataset and ask for the alias pattern

> attach(rr)
> alias(Resp ~ (V1+V2+V3+V4+V5)^5 )
Model
Resp ~ (V1 + V2 + V3 + V4 + V5)^5
Complete
               (Intercept) V1 V2 V3 V4 V5 V1:V2 V1:V3
         V1:V4                          1
         V1:V5                       1
         V2:V3                       1
         V2:V4                    1
         V2:V5                                   1
         V3:V4                 1
         V3:V5                             1
         V4:V5              1
      V1:V2:V3                          1
      V1:V2:V4                                   1
      V1:V2:V5                    1
      V1:V3:V4                             1
      V1:V3:V5                 1
      V1:V4:V5  1
      V2:V3:V4  1
      V2:V3:V5              1
      V2:V4:V5                             1
      V3:V4:V5                                   1
   V1:V2:V3:V4              1
   V1:V2:V3:V5  1
   V1:V2:V4:V5                 1
   V1:V3:V4:V5                    1
   V2:V3:V4:V5                          1
V1:V2:V3:V4:V5                       1

Lets follow through on this. We have a 1/4 fraction of the full design so we expect (if the design was well done) that each estimable effect will be aliased with three other effects. The output tells us that the estimable effects (the first row) are the main effects and a couple of two-way interactions. All the other two-way, three-way and higher interactions are confounded/aliased with the estimable effects. So, V1 is the same as V4:V5 V2:V3:V5 and V1:V2:V3:V4. This is a special kind of design, consructed with considerable care.

Example 13.2

It should be clear from the text that one can “analyse” these data without much work, but some supporting calculations are as follows.

> e2
  V1 V2 V3 V4 Resp
1  0  0  0  0   20
2  1  0  0  1   14
3  0  1  0  1   17
4  1  1  0  0   10
5  0  0  1  1   19
6  1  0  1  0   13
7  0  1  1  0   14
8  1  1  1  1   10
> anova(mymod <- aov(Resp ~ V1 + V2 + V3 + V4, data=e2))
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F)
V1         1    66.125 66.12500 144.2727 0.0012415
V2         1    28.125 28.12500  61.3636 0.0043321
V3         1     3.125  3.12500   6.8182 0.0796050
V4         1     1.125  1.12500   2.4545 0.2151699
Residuals  3     1.375  0.45833
> anova(mymod <- aov(Resp ~ (V1 + V2 + V3 + V4)^2, data=e2))
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
      Df Sum of Sq Mean Sq
V1     1    66.125  66.125
V2     1    28.125  28.125
V3     1     3.125   3.125
V4     1     1.125   1.125

V1:V2  1     0.125   0.125
V1:V3  1     1.125   1.125
V1:V4  1     0.125   0.125

To treat it as a replicated design on just the first two variables, the appropriate ANOVA is

> anova(mymod <- aov(Resp ~ V1 + V2, data=e2))
Analysis of Variance Table
Response: Resp
Terms added sequentially (first to last)
          Df Sum of Sq Mean Sq  F Value       Pr(F)
V1         1    66.125  66.125 58.77778 0.000601443
V2         1    28.125  28.125 25.00000 0.004104716
Residuals  5     5.625   1.125

Discrete—Categorical—Qualitative—Frequency—Data

References

• S. E. Fienberg, The Analysis of Cross-Classified Categorical Data, MIT Press, 1980. (a.k.a., the Xmas-Tree book)

• R. L. Plackett, The Analysis of Categorical Data, Second Edition, Griffin, 1981.

• D. R. Cox, Analysis of Binary Data, several editions, Chapman and Hall.

• S. Haberman, Analysis of Qualitative Data, Volumes 1 and 2, Academic Press.

• Y. M. M. Bishop, S. E. Fienberg, and P. Holland, Discrete Multivariate Analysis, MIT Press, 1975.

• S. Haberman, The Analysis of Frequency Data, University of Chicago Press, 1974. This is a very difficult, theoretical book.

• W. N. Venables and B. D. Ripley, Modern Applied Statistics with S-Plus, Chapter 7 on Generalized Linear Models.

• A. Agresti, Lost title, recent....

The basic data is a table (of list) of counts, i.e., we will consider cases where the response is categorical and the explanatory variables are also categorical (also known as factors). The following table shows how this fits into the general scheme of statistics, and in particular where logistic regression (which we will talk about) and regression/ANOVA (which we have talked about) fit in.

Explanatory Variables

Response Categorical Continuous Both
Categorical Contingency Tables Logistic Regression Loglinear Model
Continuous ANOVA/MANOVA Linear Regression ANACOVA
Both  ?/Hard  ?/Hard  ?/Hard

History

This subject has a long history, perhaps even longer than the history of regression. There have been lots of famous disputes about appropriate methods that have led to lots of interesting statistics. Look at any of the books.

Motivation

Consider the following table. These data, as are many of the examples, are stolen from the Xmas-Tree book.

The response variable (Heart Disease) is binary and the other variables have been broken up into 4 categories each.

Coronary Serum  Systolic
Heart Cholesterol Blood Pressure (mm Hg)
Disease (mg/100 cc) < 127 127–146 147–166 167+
< 200  2  3  3  4
Present 200–219  3  2  0  3
220–259  8  11  6  6
>= 260  7  12  11  11
< 200  117  121  47  22
Absent  200–219  85  98  43  20
220–259  119  209  68  43
>= 260  67  99  46  33

Exercise Understand the data in this table—which are the most and least risky groups.

One dimensional tables—Poisson and Multinomial Distributions

Recall the properties of the Poisson and Binomial distributions.

The Poisson distribution with mean λ has density

Poi(λ) : Pr(R = r) = (e^{−λ}λ^{r})/r!

and that E(R) = λ and Var(R) = λ.

The binomial distribution with parameter p, Bi(1, p) has properties, Pr(R = 1) = p, Pr(R = 0) = 1 − p, E(R) = p, and Var(R) = p(1 − p).

Now assume we have a response variable R which can take on 5 possible values. The example I have in mind is a classification of the fingerprints of the left first fingers of males. Here is some sample data (drawn from the book by Plackett).

Arch Small Loop Large Loop Composite Whorl 41  139  53  28  68

There are (at least) two ways I could have generated the above data.

1. I decided I wanted 329 observations, and then found them and classified them (i.e., the sample size was fixed).

2. I collection observations until I got tired (ran out of time/money/) and ended up with 329 people

Either way, if I selected people at random we should be able to assume that the individual observations are independent. The Independence is important but very difficult to check. The following scenario is an example where independence may not be appropriate. Suppose we actually selected 100 families (an got 329 men in the sample). If there is a correlation between (say) race an finger print type, then our independence assumption may not be valid.

What questions are of interest? We may wish to know what is the proportion of our population with (say) small loops? Do you have any good estimates? What makes your estimate good?

Sampling Schemes and Likelihoods

Assume sampling method 1 (that is a fixed sample size), and classification probabilities of p_{1},p_{2},p_{3},p_{4}, and p_{5}, naturally with the constraint that ∑ p_{i} = 1.

Then the probability that I observe the data that I did (this is called the likelihood) is

P(R_{1} = 41, …, R_{5} = 68) = \binom{329}{41, 139, 53, 28, 68} · p_{1}^{41} · p_{2}^{239} · p_{3}^{53} · p_{4}^{28} · p_{5}^{68}

where \binom{329}{41,139,53,28,68} is the multinomial coefficient which is defined to be equal to

\frac{329!}{41!139!53!28!69!}

If we now ask what is the maximum likelihood estimate of, say p_{1}, then it is an easy exercise to show that

p̂_{1} = \frac{41}{139}

Similarly for the other p̂_{i}.

What happens under sampling method 2, i.e., where the total number of observations is a random variable? The multinomial distribution is not longer appropriate, but each of the counts can be viewed as one observation on a Poisson random variable. Let m_{1}–M_{5} denote the mean values of Poisson distributions. If we assume that the data are generated independently and the distribution of each category is Poisson, then

P(R_{1} = 41, …, R_{5} = 68) = \frac{e^{−m_{1}}m_{1}^{41}}{41!} \frac{e^{−m_{2}}m_{2}^{139}}{139!} \frac{e^{−m_{3}}m_{3}^{53}}{53!} \frac{e^{−m_{4}}m_{4}^{28}}{28!} \frac{e^{−m_{5}}m_{5}^{68}}{68!}
= exp −(m_{1} + m_{2} + m_{3} + m_{4} + m_{5})\frac{∏ m_{i}^{r_{i}}}{∏ r_{i}!}

It is (should be) obvious that m̂_{i} = r_{i} and the MLE of the proportion in each category remains the same.

It is a fact that unless we wish to make inferences about population sizes, then it doesn’t matter which sampling method we use. Mathematically this means that,

P_{Poi}(R_{1}, …, R_{5}|R_{1} + R_{2} + R_{3} + R_{4} + R_{5} = 329) = P_{Mult}(R_{1}, …, R_{5})

We haven’t finished with this simple example yet. We need an estimate of Var(p̂_{i}) so that we have a scale of comparison. Recall that p̂_{i} = \frac{r_{i}}{n} and note that in our two sampling cases:

1.

Var(R_{i}) = np(1 − p)

so

V̂ar(p̂_{i}) = \frac{p̂(1 − p̂)}{n}

2. Var(R_{1}) = m_{1}, but p_{1} = \frac{R_{1}}{R_{1}+···+R_{5}} . It turns out that

Var(p_{i}) = \frac{(\frac{m_{i}}{R_{+}} )(1 − \frac{m_{i}}{R_{+}} )}{R_{+}}

which is equivalent to the multinomial case—because of the conditional distribution argument.

Another Example

We examine 3343 families with 4 or more children and find out how many boys in the first 4 children. The data are

0 1 2 3 4 183 789 1250 875 246

So far we have not considered models. What is an appropriate model for these data? If we assume (a little erroneously) that Pr(Boy) = 0.5, and that births within a family are independence then we might generate the following probabilities,

0 1  2 3 4
\frac{1}{2}^{4} 4\frac{1}{2}^{4} 6 \frac{1}{2}^{4} 4\frac{1}{2}^{4} \frac{1}{2}^{4}
\frac{1}{16}  \frac{1}{4}  \frac{3}{8}  \frac{1}{4} \frac{1}{16}

which would give us an expected table of counts of

0  1  2  3  4 208.9 835.75 1253.6 837.75 208.9

and then we need to ask how well does the data fit the model.

One test of this is called Pearson’s X^{2} which is

X^{2} = ∑ \frac{(O − E)^{2}}{E}
= \frac{(183 − 208.9)^{2}}{208.9} + ⋯ + \frac{(246 − 208.9)^{2}}{208.9}
= 14.27

Asymptotically (i.e., as the sample size gets big) this is distributed as a χ_{4}^{2} random variable. Why 4 df? Note that χ_{4}^{2}(0.05) = 9.49, and χ_{4}^{2}(0.01) = 13.28.

2 × 2 Tables

Example 1 Linus Pauling (Nobel laureate and Vitamin C fame) reports the following double blind data on French skiers and their propensity to get colds.

Cold No Cold
Placebo  31  109 140
Vit C (1 gm/day)  17  122 139

How was (could have) these data been collected? Two groups were selected, one was given treatment and the other the placebo. The aim was to compare the probability of getting a cold in each group. In other words we have two binomials, one for Placebo and one for Vitamin C and we wish to compare responses. Here getting a cold is the response variable and the treatment is the explanatory variable.

Example 2 The following data is the cross-classification of Perch Height by Perch Diameter of 164 lizards.

Height Diameter (in)
(feet) <= 4.0 > 4.0
> 4.75  32  11
<= 4.75  86  35

How were these data collected? Is the Poisson or Multinomial distribution appropriate? What sort of question are appropriate for these data?

• Is the proportion at height 1 > height 2?

• Is the proportion at diameter 1 > diameter 2?

• Is height independent of diameter?

Whether we can answer these questions depends on how the data were sampled?

Two Binomial—Skiers and Vitamin C

We have two groups of skiers—those that were given Vitamin C and those that were not. We wish to ask the question: Is the probability of getting a cold the same in both groups?

If the data is,

Cold No Cold
Placebo  31  109 140
Vit C (1 gm/day)  17  122 139

then a reasonable model for the data might be

Cold No Cold
Placebo  p_{1} (1 − p_{1})
Vit C (1 gm/day)  p_{2} (1 − p_{2})

and the question becomes: Is p_{1} = p_{2} = p? A good estimate of the common p is \frac{(31+17)}{279}. But how do we test if this common value of

p adequately describes the data?

Consider the generic table of data,

\begin{matrix}x_{1,1} & x_{1,2} & n_{1} \\ x_{2,1} & x_{2,2} & n_{2} \end{matrix}

Method 1 Note that p̂_{1} = \frac{x_{1,1}}{n_{1}} and p̂_{2} = \frac{x_{2,1}}{n_{2}} and we wish to test if

p_{1} − p_{2} = 0

Note that E(p̂_{1} − p̂_{2}) = p_{1} − p_{2} and that

Var(p̂_{1} − p̂_{2}) = Var(p̂_{1}) + Var(p̂_{2})
= \frac{p_{1}(1 − p_{1})}{n_{1}} + \frac{p_{2}(1 − p_{2})}{n_{2}}

If we assume that p_{1} = p_{2} then

Var(p̂_{1} − p̂_{2}) = (\frac{1}{n_{1}} + \frac{1}{n_{2}} )(1 − p)p

If we further assume that p̂_{1} − p̂_{2} is approximately normal with mean 0 and variance Var(p̂_{1} − p̂_{2}), then we can test the hypothesis with the statistic

\frac{p̂_{1} − p̂_{2}}{\sqrt{(\frac{1}{n_{1}} + \frac{1}{n_{2}} )(1 − p̂)p̂}} ∼ N(0, 1)

For the skiers example, this is = 2.19, which is quite a large number.

We had to assume approximate normality. Numerous numerical (and theoretical) studies have shown this is a quite valid assumption even when the counts in each “cell” of the table become quite small.

Method 2 An alternative test for p_{1} = p_{2} can be obtained by using Pearson’s X^{2},

X^{2} = ∑ \frac{(O − E)^{2}}{E}

where O are the observed values and E are the expected values under the model. For our type of data the expected table is

\begin{matrix}n_{1}p̂ & n_{1}(1 − p̂) \\ n_{2}p̂ & n_{2}(1 − p̂) \end{matrix}

which is equal to

24.1 115.9
23.9 115.1

and thus X^{2} = 4.81 = (2.19)^{2}. For this example, the X^{2} statistic is the square of a N(0, 1)—so its distribution is χ_{1}^{2}. In general, X^{2} is a useful statistics that is distributed as a χ_{∗}^{2}. There is some uninteresting algebra which show the equivalence of the two tests.

Why is it on 1 degree of freedom? Note that the original data had only 2 degrees of freedom because the sample sizes were predetermined. We used the data to estimate one parameter, p, which left 1 degree of freedom for the test.

This test for equality of Binomials is more generally referred to as a test of Homogeneity of Proportions. Note the importance of identifying the Response and Explanatory variables.

Poisson or Multinomial—Lizards on Perches

We have one group of lizards which are cross-classified according to the perch height and diameter. Sampling may have been either multinomial (we knew how many lizards there were) or Poisson (we didn’t know). Notice we have 2 Response variables, and not explanatory variables, so appropriate questions about the data concern the interrelationships between the responses.

Height Diameter (in)
(feet) <= 4.0 > 4.0
> 4.75  32  11
<= 4.75  86  35

The data and probabilities of observation are, Consider the generic table of data,

\begin{matrix}x_{1,1} & x_{1,2} & n_{1} \\ x_{2,1} & x_{2,2} & n_{2} \end{matrix}  probabilities  \begin{matrix}p_{1,1} & p_{1,2} & p_{1,+} \\ p_{2,1} & p_{2,2} & p_{2,+} \\ p_{+,1} & p_{+,2} & p_{+,+} \end{matrix}

We may like to ask questions like “Is Diameter independent of Height?”. For independence we know that marginal probabilities need to be multiplied. If we take the marginal Height and Diameter probabilities as given then the independence hypothesis say that

p̂ = \begin{matrix}p_{+,1}p_{1,+} & p_{+,2}p_{1,+} \\ p_{+,1}p_{2,+} & p_{+,2}p_{2,+} \end{matrix}

should be the observed probabilities under independence.

In other words, the table of fitted values under the independence model is

\begin{matrix}\frac{x_{+,1}x_{1,+}}{x_{+,+}} & \frac{x_{+,2}x_{1,+}}{x_{+,+}} \\ \frac{x_{+,1}x_{2,+}}{x_{+,+}} & \frac{x_{+,2}x_{2,+}}{x_{+,+}} \end{matrix} = 30.9 12.1 87.1 33.9

The value of X^{2} for these data is 0.18 which should be compared to a χ_{1}^{2}—not very large.

Loglinear Models

Notice that for both of the 2 × 2 table examples the fitted values were of the form

m̂_{i,j} = \frac{x_{i,+} · x_{+,j}}{x_{+,+}}

Now, take logarithms of this expression and notice that

log(m̂_{i,j}) = − log(x_{+,+}) + log(x_{i,+}) + log(x_{+,j})

In other words, the fitted values are linear functions of the margins on the logarithmic scale.

Now lets consider the model for the fitted values which says,

log(m_{i,j}) = u + u_{1(i)} + u_{2(j)} + u_{12(i,j)}

with the constraints that ∑_{i} u_{1(i)} = 0, and ∑_{j} u_{2(j)} = 0, and ∑_{i} u_{12(i,j)} = ∑_{j} u_{12(i,j)} = 0. These models are the analogue of ANOVA models and are “natural” models for Pois-

son/Multinomial data. In particular, the notions of independence and conditional independence can be easily considered using these models.

Another way of viewing the loglinear model is to say that it represents the logarithm of the expected values (m_{i,j}) in terms of main effects and interactions—which is just what the ANOVA models did.

A third, and most general, way of describing the models is to say that the loglinear models posit that

log(m) ∈ M

where M is a linear space which describes the hypothesis of interest. That is about the same as saying that an ANOVA model is just a linear model, and just as in the ANOVA situation it ignores the usefulness of various forms for the linear space M.

The first formulation, in terms of the u-term parameterization, is one useful way of describing the linear space M.

To make this a little more concrete consider a 3 × 3 table,

\begin{matrix}x_{1,1} & x_{1,2} & x_{1,3} \\ x_{2,1} & x_{2,2} & x_{2,3} \\ x_{3,1} & x_{3,2} & x_{3,3} \end{matrix}

where the two categorizing variables are both responses. (as in the lizard example). Lets assume that the expected values of the response (under some model to be determined) are

m = \begin{matrix}m_{1,1} & m_{1,2} & m_{1,3} \\ m_{2,1} & m_{2,2} & m_{2,3} \\ m_{3,1} & m_{3,2} & m_{3,3} \end{matrix}

The model which asserts independence of the two categories is

ln(m) ∈ span \{ \begin{matrix}1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \\ & u \end{matrix} , \begin{aligned}1 & 1 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ & u_{1(1)} \end{aligned} , \begin{aligned}0 & 0 & 0 \\ 1 & 1 & 1 \\ 0 & 0 & 0 \\ & u_{1(2)} \end{aligned} , \begin{aligned}0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 1 & 1 \\ & u_{1(3)} \end{aligned} , \begin{aligned}1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ & u_{2(1)} \end{aligned} , \begin{aligned}0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ & u_{2(2)} \end{aligned} , \begin{matrix}0 & 0 & 1 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \\ & u_{2(3)} \end{matrix} \}

That is, any set of fitted values which can be written as a loglinear combination of the above vectors satisfies independence. One can also see from this discussion where the constraints of our particular parameterization come from. (I will probably write out a few more of these gory spanning vectors, but not too many. They are just too hard for me to type!).

In general, we will oscillate between talking about abstract loglinear models, specific parameterizations, and the compromise model formulae of a package like S-Plus. Mostly (hopefully all the time, but I’ll hedge a little), the model embodies some underlying notion and one should be careful to talk about the underlying model and not over-interpret the specific parameter estimates.

Now that we have written a model in this general (loglinear) form, an obvious question arises.

• Why not model m ∈ M, a linear model?

or on some other scale?

The main reason is that if the loglinear model is written as a model for probabilities (which come directly from the fitted values), then we can be sure that the estimated probabilities will lie between 0 and 1. Without a lot of extra constraints and complications, linear models can not guarantee this. A more (or maybe less) compelling reason is that under the loglinear model the likelihood function factorizes in a natural way.

Bigger Two-way Tables

Consider the following cross-classification of 1019 school children into a teachers rating of schoolwork and an “independent” rating of homework conditions. (Independent in the nontechnical sense.) We would like to ask whether the two responses are independent of each other.

Homework Conditions

Rating A B C D E
A  141 67 114 79 39 440
B  131 66 143 72 35 447
C  36 14 38 28 16 132
308 147 295 179 90 1019

Under independence the fitted values are generated by multiplying the row margins but the column margins and dividing by the grand total. This doesn’t have to be obvious, but it is true.

Looking at the data we see that as homework conditions get worse (A to E) more people get a bad teachers rating. Is this significant, or am I imagining it? The table of fitted values, under independence is.

133 63.5 121.4 77.3 38.9
135.1 64.5 129.4 78.5 39.5
39.9 19.0 38.2 23.2 11.6

with X^{2} = 8.21 on (3 − 1) × (5 − 1) = 8 degrees of freedom. The value of X^{2} tells us that there is not strong reason to reject the assumption of independence. What about looking at the residuals (generally a good thing to do), but which residuals? (there are at least two commonly used residuals …).

The observations x_{i,j} − m̂_{i,j} have different variances (which is approximately m̂_{i,j}) so it makes sense to look at the Pearson standardized residuals

r_{i,j} = \frac{x_{i,j} − m̂_{i,j}}{\sqrt{m̂_{i,j}}}

For these data these numbers are

0.69 0.44 -0.67 0.19 0.02
-0.35 0.19 1.20 -0.73 -0.72
-0.62 -1.15 0.0 1.0 1.29

None of these residuals is very large (compare to N(0, 1) nor does there appear to be any strong pattern to them. It seems that it is pretty safe to assume that the homework conditions and teachers ratings are independent.

Informal Exercise

1. Verify that the table of fitted values can be written as log(m_{i,j}) = μ+μ_{1(i)}+μ_{2(j)}+μ_{12(i,j)} with the constraints that ∑_{i} μ_{1(i)} = 0, and ∑_{j} μ_{2(j)} = 0.

2. Use the residuals to calculate X^{2}.

Example The following example is drawn from Problem 2.2 in the Xmas-Tree book. The data below, from the National Crime Survey, represent repeat-victimization information (i.e., successive pairs of victimization for households in the survey). If, over th course of the period in question, a household is victimized m times, then these m victimizations will lead to m−1 entries in the table, one for each successive pair of victimizations. Here A is assault; B is burglary; HL is household larceny; MV is motor vehicle theft; PL is personal larceny; PP/PS is pocket picking and purse snatching; Ra is rape; and Ro is robbery.

Second Victimization

First  Ra  A Ro PP/PS  PL  B  HL MV Totals
Ra  26 50 11  6  82 39  48 11  273
A  65 2997 238  85 2553 1083 1349 216 8586
Ro  12 279 197  36 459 197 221 47 1448
PP/PS 3 102 40  61 243 115 101 38  703
PL  75 2628 413  229 12137 2658 3689 687 22526
B  52 1117 191  102 2649 3210 1973 301 9595
HL  42 1251 206  117 3757 1962 4646 391 12372
MV  3 221 51  24 678 301 367 269 1914
Total 278 8645 1347  660 22558 9565 12394 1960

There are many things to ask of these data, including:

• Test these data for homogeneity of proportions

• Examine a table of the square roots of the contribution from each cell to X^{2}, i.e. (Observed − Expected)/\sqrt{Expected} Which cell stands out?

• Provide an interpretation to your findings above.

• Why might you not have bothered to do the cull calculation of X^{2} after you found the contribution of the (1, 1) cell?

Here is a cursory set of answers to these questions. A quick examination of the data shows that the row and colun margins are very similar and the diagonal cells of teh table (corresponding to repeat victimizations of the same type) clearly dominate the table.

Thus prior to formal modelling we have good reason to expect that the model of homogeneity of row proportions will not be appropriate.

Under the model of homogeneity of row proportions, with data x_{i,j}, and estimated expected values

m̂_{i,j} = (x_{i,+} · x_{+,j})/x_{+,+}

then the fitted values are as follows.

Reiss <- fac.design(c(8,8),factor.names=c("Second","First"))
Reiss$Resp <- scan("reiss.dat")
> Reiss[1:10,]
Second First Resp
1  Scnd1 Frst1   26
2  Scnd2 Frst1   50
3  Scnd3 Frst1   11
4  Scnd4 Frst1    6
5  Scnd5 Frst1   82
6  Scnd6 Frst1   39
7  Scnd7 Frst1   48
8  Scnd8 Frst1   11
9  Scnd1 Frst2   65
10  Scnd2 Frst2 2997

Now, fit the model and look at the Analysis of Deviance statistics,

> mymod<-glm(Resp ~ First + Second, family=poisson, data=Reiss)
> anova(mymod)
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev
  NULL                    63   127032.9
 First  7 58883.15        56    68149.7
Second  7 59387.44        49     8762.3

What do the fitted values look like?

> t(matrix(round(fitted(mymod),digits=3),8,8))
        [,1]     [,2]    [,3]    [,4]     [,5]     [,6]     [,7]    [,8]
[1,]   1.322   41.111   6.406   3.139  107.275   45.487   58.940   9.321
[2,]  41.579 1292.978 201.462  98.712 3373.857 1430.576 1853.692 293.145
[3,]   7.012  218.056  33.976  16.647  568.990  241.262  312.619  49.438
[4,]   3.404  105.866  16.495   8.082  276.243  117.132  151.776  24.002
[5,] 109.036 3390.716 528.316 258.863 8847.631 3751.555 4861.137 768.745
[6,]  46.465 1444.924 225.137 110.312 3770.342 1598.693 2071.532 327.594
[7,]  59.913 1863.117 290.297 142.239 4861.560 2061.389 2671.078 422.407
[8,]   9.269  288.232  44.910  22.005  752.104  318.906  413.227  65.348
>

Notice the “trick” to get the data to print out in a nice tabular format.

The goodness-of-fit statistics for the model of homogeneity of proportions is X^{2} = 11, 189 and G^{2} = 8, 762, both of with should be asymptotically distributed as a χ_{49}^{2} random variable. The likelihood ratio statistic comes directly from the output,

> anova(mymod)
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev
  NULL                    63   127032.9
 First  7 58883.15        56    68149.7

Second  7 59387.44        49     8762.3
>

The Pearson statistic requires a custom function. I wrote a little function, called Pearson, to print the statistics and degrees of freedom.

> Pearson(mymod)
Degrees of Freedom: 64 Total; 49 Residual
Residual Pearson X^2 11189.88

The Pearson function itself if,

> Pearson
function(mymod)
{
# This is not completely bullet proof.
        nobs <- length(mymod$residuals)
        rdf <- mymod$df.resid
        if(is.null(rdf))
                rdf <- nobs - rank
        cat("Degrees of Freedom:", nobs, "Total;", rdf, "Residual\n")
        cat("Residual Pearson X^2", format(sum((mymod$y - fitted(mymod))^2/
                fitted(mymod))), "\n")
        invisible(mymod)
}

It also makes a lot of sense to examine the residuals from this model. There are many different residuals one could consider, but a natural candidate is the so-called pearson residual which is the component of Pearson’s X^{2},

(Observed − Expected)/\sqrt{Expected}

For these data, the numbers are

> t(matrix(round( residuals(mymod,type="pearson"), dig=3), 8,8))
        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]
[1,]  21.457   1.386   1.815   1.615  -2.440  -0.962  -1.425   0.550
[2,]   3.632  47.388   2.574  -1.380 -14.132  -9.189 -11.722  -4.506
[3,]   1.883   4.127  27.967   4.743  -4.611  -2.850  -5.182  -0.347
[4,]  -0.219  -0.376   5.787  18.613  -2.000  -0.197  -4.121   2.857
[5,]  -3.259 -13.098  -5.017  -1.856  34.970 -17.854 -16.812  -2.948
[6,]   0.812  -8.627  -2.275  -0.791 -18.262  40.299  -2.165  -1.469
[7,]  -2.314 -14.181  -4.947  -2.116 -15.842  -2.189  38.213  -1.528
[8,]  -2.059  -3.960   0.909   0.425  -2.702  -1.003  -2.274  25.192

The contribution to the overall X^{2} from just the first component of the residuals is (21.457)^{2} = 460.40, which is much bigger than any upper bound of a χ_{49}^{2}.

Digression on Measures of Association

Some people have asked about things like correlation for tables and discrete data. There is a literature on this (it isn’t very interesting) and even entire books. For these short remarks I’ll concentrate on the 2 × 2 case. It is easier!

First, some general comments. Consider an I × J table of counts. We would like to have a number (like the correlation

coefficient) to summarize the dependence between the two classifying quantities. Desirable properties of the measure would be:

• if p_{i,j} = p_{i,+} · p_{+,j} then the association should be 0. This corresponds to independence, which we like to think of as uncorrelated.

• if p′_{i,j} = t_{i} · s_{j} · p_{i,j} then the association for the table p′ should be the same as the association for the table p. In other words, we change the marginal distributions then the association should not change.

• If the rows of columns are permuted the measure of association should not change.

• The measure should assume values between −1 and 1 or between 0 and 1. This seems to be a historical requirement in an attempt to draw analogies with the correlation coefficient.

Association in 2 × 2 tables

Consider the 2 × 2 table,

\begin{matrix}x_{1,1} & x_{1,2} \\ x_{2,1} & x_{2,2} \end{matrix}

Here are some measures of association.

1. The Cross-product ratio. Let

α = \frac{x_{1,1} · x_{2,2}}{x_{1,2} · x_{2,1}}

First, verify that α = 0 iff the table represents independence. Note that

log α = log x_{1,1} − log x_{1,2} − log x_{2,1} + log x_{2,2}

= (1, −1, −1, 1) · ( \begin{matrix}log & x_{1,1} \\ log & x_{1,2} \\ log & x_{2,1} \\ log & x_{2,2} \end{matrix} )

This is just like the interaction contrast term for an ANOVA model with 2 factors at 2 levels. Seems to make sense.

Also, consider

\begin{matrix}p_{1,1} & p_{1,2} \\ p_{2,1} & p_{2,2} \end{matrix}

as a table of probabilities. Then

α = \frac{p_{1,1} · p_{1,2}}{p_{2,1} · p_{2,2}}

i.e., it is the ration of the odds of being in either category. This measure is the one that best fits into the loglinear model framework.

2. Yule’s Q (named for Quetelet). This is defined as

Q = \frac{p_{1,1} · p_{2,2} − p_{1,2} · p_{2,1}}{p_{1,1} · p_{2,2} + p_{1,2} · p_{2,1}}
= \frac{α − 1}{α + 1}

Another measure proposed by Yule is

Y = \frac{\sqrt{p_{1,1} · p_{2,2}} − \sqrt{p_{1,2} · p_{2,1}}}{\sqrt{p_{1,1} · p_{2,2}} + \sqrt{p_{1,2} · p_{2,1}}}

In some circumstances Yule’s Q has an easy interpretation. Consider the table to be

Response at time 1

Response at p_{1,1}  p_{1,2}
time 2  p_{2,1}  p_{2,2}

Now look at Q as

Q = \frac{p_{1,1}p_{2,2}}{p_{1,1}p_{2,2} + p_{1,2}p_{2,1}} − \frac{p_{1,2}p_{2,1}}{p_{1,1}p_{2,2} + p_{1,2}p_{2,1}}

which can be read as the conditional probability that of two observations one will have the same response minus the conditional probability that of two observations one will have a different response.

3. Correlation coefficient like measures. If we assume that the categories have values of 0 and 1, then the Pearson correlation is

ρ = \frac{p_{1,1} · p_{2,2} − p_{1,2} · p_{2,1}}{\sqrt{p_{1,+}p_{2,+}p_{+,2}p_{+,2}}}

Note that under the model of independence X^{2} = ρ^{2} . There is a whole class of measures based on ρ^{2}, e.g., P = \sqrt{\frac{ρ^{2}}{ρ^{2}+1}} and R_{1} = \frac{ρ^{2}}{x_{+,+}} , etc.

End Digression on Measures of Association.

Three-dimensional Tables

Consider data which might come in a table that looks like,

\begin{matrix}x_{1,1,1} & x_{1,1,2} & x_{1,1,3} \\ x_{1,2,1} & x_{1,2,2} & x_{1,2,3} \end{matrix}  \begin{matrix}x_{2,1,1} & x_{2,1,2} & x_{2,1,3} \\ x_{2,2,1} & x_{2,2,2} & x_{2,2,3} \end{matrix}

where i indexes the faces of a cube, j indexes the rows and k indexes the columns.

The first question to ask is what sort of sampling schemes could have generated data like these? and consequently, what sort of (conditional) independence statements can one make?

There are many ways that these data could have been sampled/generated. I’ll list some of these here.

1. Multinomial sampling. Assume that the total number of observations in the table x_{+,+,+} is fixed and we just cross-classify these into I × J × K categories. (This implies 3 response variables)

2. Poisson Each cell in the table is an independent Poisson observation, so the grand total is not fixed. (This too implies 3 response variables)

3. Product-Multinomial. Lets assume that the category i is an explanatory variable, categories j and k are responses. That is we know x_{i,+,+} for i = 1, 2. Thus the two faces of the table are multinomials. We are interested in comparisons between the two populations, such as how does the interaction change from one population to the other.

4. Product-Multinomial. Both i and j are explanatory variables, so we know x_{i,j,+}. That is we sample I × J multinomials, each with K cells.

Each of these cases corresponds to a different probability distribution. Each of the following looks different, but they are all very similar, modulo some conditioning arguments. The probability distributions are:

1.

P(x_{i,j,k}) = \binom{x_{+,+,+}}{x_{1,1,1}, ⋯, x_{I,J,K}} ∏_{i,j,k}(p_{i,j,k}^{x_{i,j,k}})

with the constraint that ∑_{i,j,k} p_{i,j,k} = 1. That is, we have one big multinomial with probabilities, p_{i,j,k}.

2.

P(x_{i,j,k}) = ∏_{i,j,k} \frac{e^{−m_{i,j,k}}m_{i,j,k}^{x_{i,j,k}}}{(x_{i,j,k}!)}

That is we have I × J × K independent Poissons each with mean m_{i,j,k}.

3.

P(x_{i,j,k}) = ∏_{i} [\binom{x_{i,+,+}}{x_{i,1,1}, ⋯, x_{i,J,K}} ∏_{j,k}(p_{i,j,k}^{x_{i,j,k}})]

where ∑_{j,k} p_{i,j,k} = 1 for each i. That is, we have lots of independent multinomials.

P(x_{i,j,k}) = ∏_{i,j} [\binom{x_{i,j,+}}{x_{i,j,1}, ⋯, x_{i,j,K}} ∏_{k}(p_{i,j,k}^{x_{i,j,k}})]

where ∑_{k} p_{i,j,k} = 1 for each combination of i and j. That is, we have lots and lots of independent multinomials.

It is not an accident that each of these likelihoods can be written in the form

k(x_{i,j,k}) × ∏_{i,j,k} e^{−m_{i,j,k}}m_{i,j,k}^{x_{i,j,k}}

with some appropriate constraints on the m_{i,j,k}.

Thus each of the situations is essentially the same except the likelihoods require that the Sampling Constraints be recognized. That is for the last case where we have I × J independent multinomials we must make sure that we get

∑_{i} p̂_{i,j,k} = 1, ∀i, j

Models

Let us consider the loglinear model

log(m_{i,j}) = u + u_{1(i)} + u_{2(j)} + u_{3(k)} + u_{12(i,j)} + u_{13(i,k)} + u_{23(j,k)} + u_{123(i,j,k)}

the model that include the (so-called) (i) main effects, (ii) two-factor interactions, (iii) threefactor interaction. The linear space spanned by this model is (quite complicated but) written out below.

u spanned by \{ \begin{matrix}1 & 1 & 1 \\ 1 & 1 & 1 \end{matrix} , \begin{matrix}1 & 1 & 1 \\ 1 & 1 & 1 \end{matrix} \} depends on no index

u_{1(i)} spanned by \{ \begin{matrix}1 & 1 & 1 \\ 1 & 1 & 1 \end{matrix} , \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \} , \{ \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} , \begin{matrix}1 & 1 & 1 \\ 1 & 1 & 1 \end{matrix} \} , depends only on i index

u_{2(j)} spanned by \{ \begin{matrix}1 & 1 & 1 \\ 0 & 0 & 0 \end{matrix} , \begin{matrix}1 & 1 & 1 \\ 0 & 0 & 0 \end{matrix} \} , \{ \begin{matrix}0 & 0 & 0 \\ 1 & 1 & 1 \end{matrix} , \begin{matrix}0 & 0 & 0 \\ 1 & 1 & 1 \end{matrix} \} , depends only on j index

u_{12(i,j)} spanned by \begin{matrix}\{ \begin{matrix}1 & 1 & 1 \\ 0 & 0 & 0 \end{matrix} , \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \} , \{ \begin{matrix}0 & 0 & 0 \\ 1 & 1 & 1 \end{matrix} , \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \} , \\ \{ \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} , \begin{matrix}1 & 1 & 1 \\ 0 & 0 & 0 \end{matrix} \} , \{ \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} , \begin{matrix}0 & 0 & 0 \\ 1 & 1 & 1 \end{matrix} \} \end{matrix} , depends only on i, j index

and so on for the other indices.

Obviously there need to be some constraints, of the form ∑_{i} u_{1(i)} = 0, and ∑_{j} u_{2(j)} = 0, and ∑_{i} u_{12(i,j)} = ∑_{j} u_{12(i,j)} = 0. ∑_{i,j} u_{12(i,j)} = 0.

These models are (from experience) useful ways of modelling data. They are also (from theory) useful ways of describing independence and conditional independence relationships.

FACT Let M be the model we are considering (defined by parameters, or whatever), then the sufficient statistics (i.e., minimal summary of the data) are P_{M}(x) and the likelihood equations are:

1. P_{M}(x) = P_{M}(m̂)

2. log m̂ ∈ M

EXAMPLE Consider the independence model

log(m_{i,j,k}) = u + u_{1(i)} + u_{2(j)} + u_{3(k)}

then the sufficient statistics are x_{i,+,+}, x_{+,j,+}, and x_{+,+,k}, and the likelihood equations are

x_{i,+,+} = m̂_{i,+,+}
x_{+,j,+} = m̂_{+,j,+}
x_{+,+,k} = m̂_{+,+,k}

and

log m̂_{i,j,k} = u + u_{1(i)} + u_{2(j)} + u_{3(k)}

The solution to these equations is easily (once you know the answer) seen to be

m̂_{i,j,k} = \frac{x_{i,+,+}x_{+,j,+}x_{+,+,k}}{x_{+,+,+}x_{+,+,+}}

The indices on the u-term parameterization give the form of the sufficient statistics. (That is, with some experience one can just read these things off).

This model was OK if we had three response variables, and could be interpreted as the independence model for all three factors. What if (say) variable 1 was an explanatory variable—that is the margin corresponding to that variable was fixed by the sampling scheme? If that were the case, then the first condition of the likelihood equations, x_{i,+,+} = m̂_{i,+,+} is not a model constraint but is rather a sampling constraint. In other words, we can only consider models which include the sampling constraint or in other words, include the u_{1} term. We can not talk about a model that is like

log(m_{i,j,k}) = u + u_{2(j)} + u_{3(k)}

because the fitted values need not (and probably will not) satisfy the sampling constraint.

What about the product-multinomial case where both i and j are explanatory. This implies that all valid models must include at least the u_{1(i)} + u_{2(j)} + u_{12(i,j)} terms. An example of an acceptable model is

log(m_{i,j,k}) = u_{1(i)} + u_{2(j)} + u_{3(k)} + u_{12(i,j)} + u_{13(i,k)}

with corresponding likelihood equations of

x_{i,j,+} = m̂_{i,j,+}
x_{i,+,k} = m̂_{i,+,k}

and

log(m_{i,j,k}) = u_{1(i)} + u_{2(j)} + u_{3(k)} + u_{12(i,j)} + u_{13(i,k)}

which implies MLEs of

m̂_{i,j,k} = \frac{x_{i,j,+}x_{i,+,k}}{x_{i,+,+}}

Models need not have closed-form MLEs—we were just lucky here.

Categorical Data Examples

Here is a Splus view of Problem 5.1, to complement the written write-up of the problem.

> y <- scan("t58.dat")
> r5 <- rep(rep(1:2, times=c(1,1)), length=16)
> r1 <- rep(rep(1:2, times=c(2,2)), length=16)
> anat <- rep(rep(1:2, times=c(4,4)), length=16)
> sub <- rep(rep(1:2, times=c(8,8)), length=16)
> data.frame(y,sub,anat,r1,r5)
     y sub anat r1 r5
 1 187   1    1  1  1
 2  15   1    1  1  2
 3 256   1    1  2  1
 4  42   1    1  2  2
 5  42   1    2  1  1
 6  40   1    2  1  2
 7  34   1    2  2  1
 8  62   1    2  2  2
 9 177   2    1  1  1
10  14   2    1  1  2
11 194   2    1  2  1
12  27   2    1  2  2
13  30   2    2  1  1
14  63   2    2  1  2
15  52   2    2  2  1
16 121   2    2  2  2
>
# Write a function to generate the design matrix elements more easily
> gm <- function(string, times, length) {
+ as.factor(rep(rep(string, times=times), length=length))}
> gm(1:2,8,16)
 [1] 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2

Now we have the data, lets try fitting some models.

> mymod <- glm(y ~ sub*anat*r1*r5 -sub:anat:r1:r5,fam=poisson)
> anova(mymod)
Analysis of Deviance Table
Poisson model
Response: y
Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev
       NULL                    15   936.6933

        sub  1   0.0000        14   936.6933
       anat  1 164.8921        13   771.8012
         r1  1  35.8515        12   735.9497
         r5  1 263.6337        11   472.3160
   sub:anat  1  26.0619        10   446.2540
     sub:r1  1   0.0000         9   446.2540
    anat:r1  1   1.6969         8   444.5571
     sub:r5  1  15.8852         7   428.6719
    anat:r5  1 401.8329         6    26.8391
      r1:r5  1  11.6760         5    15.1630
sub:anat:r1  1   8.3821         4     6.7810
sub:anat:r5  1   3.9663         3     2.8147
  sub:r1:r5  1   1.6838         2     1.1309
 anat:r1:r5  1   0.7276         1     0.4032

Its also a good idea to make a Pearson X2 function,

> Pearson <- function(mymod){
+ sum( (mymod$y - fitted(mymod))^2/fitted(mymod))}
> Pearson(mymod)
[1] 0.4032997

Now, look at the coefficients,

> summary(mymod, corr=F)
Call: glm(formula = y ~ sub * anat * r1 * r5 - sub:anat:r1:r5, family = poisson)
Deviance Residuals:
           1         2          3         4         5         6          7
 -0.07424731 0.2688111 0.06365658 -0.155701 0.1582356 -0.159483 -0.1727309
         8        9         10          11        12        13        14
 0.1298921 0.076601 -0.2655282 -0.07290024 0.1982652 -0.183648 0.1288457
        15         16
 0.1419844 -0.0922108
Coefficients:
                   Value Std. Error     t value
(Intercept)  11.02839521  1.6399111  6.72499569
        sub   1.24740863  0.8477184  1.47148937
       anat  -2.30346258  1.1085618 -2.07788370
         r1   0.36697573  0.8246952  0.44498346
         r5  -6.23963203  1.4216230 -4.38909047
   sub:anat  -1.62266440  0.5922377 -2.73988712
     sub:r1  -0.68462590  0.3790851 -1.80599506
    anat:r1  -1.05301778  0.6149272 -1.71242682
     sub:r5  -0.04290055  0.7340333 -0.05844497
    anat:r5   2.31075418  0.7006763  3.29789109

      r1:r5   1.45785420  0.6690421  2.17901701
sub:anat:r1   0.85622629  0.2913198  2.93912792
sub:anat:r5   0.55397168  0.3011419  1.83957041
  sub:r1:r5  -0.37356244  0.3095982 -1.20660398
 anat:r1:r5  -0.26518704  0.3121409 -0.84957491
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 936.6933 on 15 degrees of freedom
Residual Deviance: 0.403235 on 1 degrees of freedom
Number of Fisher Scoring Iterations: 3

It appears that the smallest t-value belongs to the last effect.

> mymod <- update(mymod, . ~ . - anat:r1:r5)
> summary(mymod, corr=F)
Call: glm(formula = y ~ sub + anat + r1 + r5 + sub:anat + sub:r1 + anat:r1 + sub:r5 +
        anat:r5 + r1:r5 + sub:anat:r1 + sub:anat:r5 + sub:r1:r5, family =
        poisson)
Deviance Residuals:
          1           2           3          4           5          6
 0.02494643 -0.08736951 -0.02129867 0.05274965 -0.05246496 0.05405967
          7           8         9         10         11        12         13
 0.05866558 -0.04321903 0.1822275 -0.6112582 -0.1725532 0.4790516 -0.4293209
        14        15         16
 0.3080217 0.3399899 -0.2179512
Coefficients:
                   Value Std. Error     t value
(Intercept)  10.24919539  1.3497469  7.59342048
        sub   1.21355168  0.8448334  1.43643905
       anat  -1.75927849  0.8994832 -1.95587701
         r1   0.84016530  0.6085070  1.38069952
         r5  -5.58723602  1.1659988 -4.79180263
   sub:anat  -1.61087973  0.5884003 -2.73772768
     sub:r1  -0.65801478  0.3787459 -1.73735177
    anat:r1  -1.40074511  0.4589962 -3.05175752
     sub:r5  -0.02661837  0.7190576 -0.03701842
    anat:r5   1.87681625  0.4765376  3.93844352
      r1:r5   1.07722895  0.4880211  2.20734086
sub:anat:r1   0.85404893  0.2894753  2.95033465
sub:anat:r5   0.55650858  0.3007260  1.85055021
  sub:r1:r5  -0.39545876  0.3056948 -1.29363919

(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 936.6933 on 15 degrees of freedom
Residual Deviance: 1.130875 on 2 degrees of freedom
Number of Fisher Scoring Iterations: 3

and so on removing the last item again (why), and so on until

> mymod <- update(mymod, . ~ . - sub:anat:r5)
> summary(mymod, corr=F)
Call: glm(formula = y ~ sub + anat + r1 + r5 + sub:anat + sub:r1 + anat:r1 + sub:r5 +
        anat:r5 + r1:r5 + sub:anat:r1, family = poisson)
Deviance Residuals:
           1           2          3        4         5          6        7
 0.008463085 -0.02979897 -0.4232902 1.117573 0.7708976 -0.7305349 0.335136
          8         9         10        11         12        13       14
 -0.2409824 0.2258727 -0.7479464 0.2694421 -0.6869914 -1.332779 1.040749
        15         16
 0.1658572 -0.1075457
Coefficients:
                  Value Std. Error    t value
(Intercept)  11.1385141  0.9310536  11.963344
        sub   0.5948815  0.5576064   1.066848
       anat  -3.2231285  0.6748951  -4.775747
         r1   1.2711671  0.5358550   2.372222
         r5  -5.9069071  0.4052989 -14.574200
   sub:anat  -0.6026888  0.4005110  -1.504800
     sub:r1  -0.9444908  0.3338057  -2.829463
    anat:r1  -1.1818001  0.3946812  -2.994316
     sub:r5   0.1881420  0.1474611   1.275875
    anat:r5   2.7006491  0.1480232  18.244770
      r1:r5   0.5033607  0.1508749   3.336279
sub:anat:r1   0.6927903  0.2394423   2.893350
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 936.6933 on 15 degrees of freedom
Residual Deviance: 6.780989 on 4 degrees of freedom
Number of Fisher Scoring Iterations: 3

and now the sub:r5 looks to be removed,

> mymod <- update(mymod, . ~ . - sub:r5)

> summary(mymod, corr=F)
Call: glm(formula = y ~ sub + anat + r1 + r5 + sub:anat + sub:r1 + anat:r1 + anat:r5 +
        r1:r5 + sub:anat:r1, family = poisson)
Deviance Residuals:
          1          2          3         4        5        6         7
 0.09700898 -0.3321493 -0.2465499 0.6322148 1.141603 -1.04393 0.7696995
        8         9         10         11         12        13       14
 -0.53309 0.1081195 -0.3710333 0.08684518 -0.2289323 -1.613918 1.295144
         15        16
 -0.2107835 0.1401241
Coefficients:
                  Value Std. Error    t value
(Intercept)  11.0083736  0.9251979  11.898398
        sub   0.7110673  0.5499720   1.292915
       anat  -3.3670856  0.6697008  -5.027746
         r1   1.2774576  0.5352994   2.386436
         r5  -5.6655003  0.3551315 -15.953248
   sub:anat  -0.5241251  0.3957200  -1.324485
     sub:r1  -0.9489361  0.3337764  -2.843029
    anat:r1  -1.2073411  0.3955115  -3.052607
    anat:r5   2.7247353  0.1471578  18.515741
      r1:r5   0.5093611  0.1508618   3.376343
sub:anat:r1   0.7059996  0.2391792   2.951761
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 936.6933 on 15 degrees of freedom
Residual Deviance: 8.405387 on 5 degrees of freedom
Number of Fisher Scoring Iterations: 3

Thats about as far as we can go. We try one step further and end up in trouble. One should also draw the picture as the well.

> summary(mymod2, corr=F)
Call: glm(formula = y ~ sub + anat + r1 + r5 + sub:anat + sub:r1 + anat:r1 + anat:r5 +
        r1:r5, family = poisson)
Deviance Residuals:
       Min         1Q     Median        3Q      Max
 -2.382332 -0.5972633 0.04882098 0.5460737 2.112859
Coefficients:

                   Value Std. Error     t value
(Intercept)  13.24730054  0.5327145  24.8675445
        sub  -0.75823360  0.2367596  -3.2025467
       anat  -5.09584953  0.3437450 -14.8245056
         r1  -0.13912795  0.2366472  -0.5879129
         r5  -5.66550075  0.3551402 -15.9528548
   sub:anat   0.59604781  0.1175572   5.0702773
     sub:r1  -0.02001409  0.1112065  -0.1799723
    anat:r1  -0.11536188  0.1436549  -0.8030488
    anat:r5   2.72473474  0.1471426  18.5176455
      r1:r5   0.50936175  0.1508436   3.3767543
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 936.6933 on 15 degrees of freedom
Residual Deviance: 17.12928 on 6 degrees of freedom
Number of Fisher Scoring Iterations: 3

A Categorical Data Example

Lets do the Red Dye number 2 example from Fienberg in Splus. You might try to do the same thing in SAS ( and ask questions about PROC CATMOD when you get stuck). Here is the problem and data.

The following table contains one version of the results of an experiment in which the food additive Red Dye No. 2 was fed to four different groups of rats at varying dosages. Some of the rats died during the experiment (i.e., between 0 and 131 weeks). The remainder were sacrificed at the end of 131 weeks. All rats were examined for tumors. Try to understand the data, and in particular what are the explanatory and response variables, and find some good models to describe the data.

Age at Death 0–131 weeks Terminal sacrifice
Dosage Group Low High Low  High
Cancer Present  4  7  0  7
Absent  26  16 14  14
Totals  30  23 14  21

This is a 2 × 2 table, with 3 categorizing variables. It seems reasonable to regard the dosage level as an explanatory variable with age at death and presence or absence of chance as both response variables. One could argue that the age at death variable might be an explanatory variable, but the number of observations in this group is not fixed (in advance) so for most purposes it makes sense to regard it as a response variable.

The distinction of response or explanatory variable does not change how we read in or manipulate the data, but does change the class of models we might consider.

In the following, I choose to enter the data the fancy way, using the fac.design command. Had I been lazy, I would have just generated the design by hand.

> reddye <- fac.design(c(2,2,2),c("Dose","Age","Cancer"))
> reddye$Resp <- c(4,7,0,7,26,16,14,14)
> reddye
   Dose  Age Cancer Resp
1 Dose1 Age1  Cncr1    4
2 Dose2 Age1  Cncr1    7
3 Dose1 Age2  Cncr1    0
4 Dose2 Age2  Cncr1    7
5 Dose1 Age1  Cncr2   26
6 Dose2 Age1  Cncr2   16
7 Dose1 Age2  Cncr2   14
8 Dose2 Age2  Cncr2   14

Now lets try to explore some models for these data. The command we need to use it glm and the appropriate family is Poisson.

The simplest model we can fit includes just the dose term. If we do not include dose then we risk fitting models which have fitted numbers of animals at the various dose levels different from the oberserved numbers—and that is not a good idea.

Thus to fit the simplest model (in this case as dose is explanatory),

> anova(mymod <- glm( Resp ~ Dose, family=poisson, data=reddye))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
     Df      Deviance Resid. Df Resid. Dev
NULL                          7   49.47734
Dose  1 -4.214286e-09         6   49.47734

and the corresponding Pearson X2 is

> Pearson(mymod)
Degrees of Freedom: 8 Total; 6 Residual
Residual Pearson X^2 42.72691

There is little to be said here. The dosage group is not a predictor of the raw response. We shouldn’t be surprised by this.

The next step is to fit various models to the data. Here goes with an unstructured look at things. We’ll see more structured approaches shortly.

> anova(mymod <- glm( Resp ~ Dose + Age, family=poisson, data=reddye))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
     Df Deviance Resid. Df Resid. Dev
NULL                     7   49.47734
Dose  1 0.000000         6   49.47734
 Age  1 3.707932         5   45.76941
> anova(mymod <- glm( Resp ~ Dose*Age, family=poisson, data=reddye))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
         Df Deviance Resid. Df Resid. Dev
    NULL                     7   49.47734

    Dose  1 0.000000         6   49.47734
     Age  1 3.707932         5   45.76941
Dose:Age  1 2.336721         4   43.43269
>

That is pretty boring. We have learnt that Age might be appropriate, but we have no strategy. Instead, let us try to actually struggle our way through this, like we would an ANOVA.

> anova(mymod <- glm( Resp ~ Dose*Age + Age*Cancer + Dose*Cancer, family=poi>
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev
       NULL                     7   49.47734
       Dose  1  0.00000         6   49.47734
        Age  1  3.70793         5   45.76941
     Cancer  1 32.82534         4   12.94407
   Dose:Age  1  2.33672         3   10.60735
 Age:Cancer  1  0.00739         2   10.59995
Dose:Cancer  1  7.60034         1    2.99961

Even though things are NOT orthogonal, it looks like one might wish to remove the Age times Cancer interaction (leaving the main effects)

> anova(mymod <- update(mymod, . ~ . - Age:Cancer))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev
       NULL                     7   49.47734
       Dose  1  0.00000         6   49.47734
        Age  1  3.70793         5   45.76941
     Cancer  1 32.82534         4   12.94407
   Dose:Age  1  2.33672         3   10.60735
Dose:Cancer  1  7.31735         2    3.29000

Now things are not clear, one may—or may not—remove the Dose times Age interaction. The upper 95% point of χ_{1}^{2} is

> qchisq(.95,1)
[1] 3.841459

so maybe we should just remove it.

> anova(mymod <- update(mymod, . ~ . - Age:Dose))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev
       NULL                     7   49.47734
       Dose  1  0.00000         6   49.47734
        Age  1  3.70793         5   45.76941
     Cancer  1 32.82534         4   12.94407
Dose:Cancer  1  7.31735         3    5.62672
> Pearson(mymod)
Degrees of Freedom: 8 Total; 3 Residual
Residual Pearson X^2 4.228423

Now it is clear we can’t proceed very much. Either of these last two models seems adequate and both are comparatively easy to interpret.

Oh yeah? How does one interpret these models, and what does it all mean. The final model (the one I would probably end up with) has an overall goodness of fit statistic of (likelihood ratio) 5.62 or (Pearson) 4.22, on 3 degrees of freedom. The upper 95% point of the appropriate chi-square is

> qchisq(.95,3)
[1] 7.814728

So, on an overall basis, the model fits very well. The model itself says that Age is indepedent of everything, but that Cancer does depend on Dose. We could look at the actual parameter estimates, but would not learn very much.

However, we are not yet done. We certainly want to look at some plots and diagnostics, even though we only have 8 data points. We think we have a good model, do the diagnostics confuse the issue.

Start with a plot of the model/residuals.

> par(mfrow=c(2,3))
> plot(mymod)

This doesn’t show anything at all untoward.

Another example, Graduate Admissions.

Consider the following data which gives admissions data for the six largest graduate programs at U.C. Berkeley for 1973.

Men  Women
Number Number not Number Number not
Major admitted admitted admitted admitted
A  512  313  89  19
B  353  207  17  8
C  120  205  202  391
D  138  279  131  244
E  53  138  94  299
F  22  351  24  317

Is there an issue of sex bias in admissions, or not?

One can argue about what is a response or explanatory variable, but to me it makes most sense to consider both Major and Sex as explanatory variables and consider just Admit as a response.

Start by reading in the data, then think about how you want to answer the question.

> berk <- fac.design( c(2,2,6), fac=list(Admit=c("Yes","No"),
+ Sex=c("Men","Women"), Major=c("A","B","C","D","E","F")))
> berk$Resp <- scan("/tmp/berk")
> berk
   Admit   Sex Major Resp
 1   Yes   Men     A  512
 2    No   Men     A  313
 3   Yes Women     A   89
 4    No Women     A   19
 5   Yes   Men     B  353
 6    No   Men     B  207
 7   Yes Women     B   17
 8    No Women     B    8
 9   Yes   Men     C  120
10    No   Men     C  205
11   Yes Women     C  202
12    No Women     C  391
13   Yes   Men     D  138
14    No   Men     D  279
15   Yes Women     D  131
16    No Women     D  244
17   Yes   Men     E   53
18    No   Men     E  138
19   Yes Women     E   94
20    No Women     E  299
21   Yes   Men     F   22
22    No   Men     F  351
23   Yes Women     F   24
24    No Women     F  317

The first thing one might want to do is to compare the admission rates for men and women. One might start by comparing the percent of women who are admitted and the percent of men admitted.

> sum(Resp[Sex=="Men"])
[1] 2691
> sum(Resp[Sex=="Men" & Admit=="Yes"])
[1] 1198
>
> sum(Resp[Sex=="Women"])

[1] 1835
> sum(Resp[Sex=="Women" & Admit=="Yes"])
[1] 557
>
> 1198/2691
[1] 0.4451877
> 557/1835
[1] 0.3035422

That is quite revealing. Only 30% of the women are admitted while almost 45% of the men are. Is it time to go and call the lawyers and argue discrimination? How could we have been mislead by this?

The next logical step is to look at the admission rates for the different programs, which are

> applied <- Resp[Admit=="Yes"] + Resp[Admit=="No"]
> applied
   1   3   5  7   9  11  13  15  17  19  21  23
 825 108 560 25 325 593 417 375 191 393 373 341
>
> round(Resp[Admit=="Yes"]/applied, dig=3)
     1     3    5    7     9    11    13    15    17    19    21   23
 0.621 0.824 0.63 0.68 0.369 0.341 0.331 0.349 0.277 0.239 0.059 0.07

Here the first two numbers are the admission rates for Men and Women in Program A, the next two for Program B, etc. There is much less difference here. This is easier to read below,

> dodo <- round(Resp[Admit=="Yes"]/applied, dig=3)
> dodo[(0:5)*2 + 1]
     1    5     9    13    17    21
 0.621 0.63 0.369 0.331 0.277 0.059
> dodo[(0:5)*2 + 2]
     3    7    11    15    19   23
 0.824 0.68 0.341 0.349 0.239 0.07

Now there appears to be much less difference between Men and Women, what is going on? Time to call off the lawyers, or have we made a mistake. Perhaps the best way to proceed from here (after we decide we think we know what is going on), is to try to fit models.

We have two explanatory variables, Sex and Major, so the minimal model must include at least the interaction between these two.

Thus our starting model is

> anova(mymod <- glm(Resp ~ Sex*Major, family=poisson))
Analysis of Deviance Table

Poisson model
Response: Resp
Terms added sequentially (first to last)
          Df Deviance Resid. Df Resid. Dev
     NULL                    23   2650.095
      Sex  1  162.874        22   2487.221
    Major  5  159.522        17   2327.699
Sex:Major  5 1220.615        12   1107.084

and we need to ask if there is any effect due to the Admit variable. But lets be careful here. What we are really asking about is the interaction between admissions rates and sex. We know, or at least should suspect that admissions rates will vary by major, so maybe the null (that is uninteresting) model should be

> anova(mymod <- glm(Resp ~ Sex*Major + Admit*Major, family=poisson))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev
       NULL                    23   2650.095
        Sex  1  162.874        22   2487.221
      Major  5  159.522        17   2327.699
      Admit  1  230.028        16   2097.671
  Sex:Major  5 1220.615        11    877.056
Admit:Major  5  855.321         6     21.736

This is a reasonable model. On an overall level, it does not completely explain that data, pchisq(21,6) = 0.9981654, but it is a great improvement over the null model. Now to answer the real question.

> anova(mymod <- glm(Resp ~ Sex*Major + Admit*Major + Admit*Sex
+ ,family=poisson))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)

            Df Deviance Resid. Df Resid. Dev
       NULL                    23   2650.095
        Sex  1  162.874        22   2487.221
      Major  5  159.522        17   2327.699
      Admit  1  230.028        16   2097.671
  Sex:Major  5 1220.615        11    877.056
Admit:Major  5  855.321         6     21.736
  Admit:Sex  1    1.531         5     20.204

Here we see that the improvement in goodness of fit by adding the interaction between admit and sex is only 1.53, with a p-value of 0.78. This can hardly be considered strong (or any evidence).

So what is going on. If it wasn’t already obvious, the men are applying to programs that have higher admission rates, so their overall admission rate is higher.

What about diagnostics for this model? Our resulting model only has 5 degrees of freedom left, so we don’t have a lot of leeway.

The one really big residual corresponds to the small number of women not admitted to program A. The model predicts that more women should not be admitted. With this discussion you also see some of the problems with residuals in this case. Although we have 24 residuals we know that there are lots of constraints imposed by the model, so the residuals are far from independent from each other.

Goodness of Fit Statistics

Thus far we have only considered one goodness of fit statistic. There are many others, just as there are many measures of association. We will only consider the important other one.

We have looked at

X^{2} = ∑ \frac{(Observed − Expected)^{2}}{Expected}

The other common goodness-of-fit statistic is the log likelihood ration statistic,

G^{2} = 2 × ∑ Observed × log (\frac{Observed}{Expected})

This statistic is really a version of

G^{2} = 2 × log \{\frac{likelihood(data|Model_{1})}{likelihood(data|Model_{2})}\}

where Model_{2} ⊂ Model_{1}. If Model_{1} is R^{n} (i.e., no model, or the full data), then the MLE is just the data itself

(x), and if for Model_{2} the MLE is m̂, then the likelihood ratio statistic above reduces to

G^{2} = 2 × ∑ x_{i,j} × log (\frac{x_{i,j}}{m̂_{i,j}} )

at least for Poisson likelihoods.

Otherwise, if MLE_{1} is m̂^{1} and Otherwise, if MLE_{2} is m̂^{2}, then the general likelihood ratio statistic becomes

G^{2} = 2 × ∑ m̂_{i,j}^{1} × log (\begin{matrix} & m̂_{i,j}^{1} \\ & m̂_{i,j}^{2} \end{matrix} )

This statistic allows you to compare two models, rather than just a model with the fixed data. In other words, this allows us to look at the improvement of one model over another as the reduction in G^{2}.

Either way the distribution of G^{2} tends towards χ_{k}^{2} where k = dimen(Model_{1})−dimen(Model_{2}). The X^{2} statistic can be similarly used to compare models, e.g.,

X_{∗}^{2} = ∑ \frac{(m̂_{i,j}^{1} − m̂_{i,j}^{2})^{2}}{m̂_{i,j}^{2}}

However G^{2} has an advantage. Consider three models, M_{3} ⊂ M_{2} ⊂ M_{1}. Then

G^{2} for comparing M_{3} to M_{2} is 2 × log \{\frac{lik(x|M_{2})}{lik(x|M_{3})}\}  (1)
G^{2} for comparing M_{2} to M_{1} is 2 × log \{\frac{lik(x|M_{1})}{lik(x|M_{2})}\}  (2)
G^{2} for comparing M_{3} to M_{1} is 2 × log \{\frac{lik(x|M_{1})}{lik(x|M_{3})}\}  (3)

Now notice that (1) + (2) = (3), i.e., the G^{2} tests add up, just like SS add up in regression. This is not true for X^{2}.

Both G^{2} and X^{2} are asymptotically equivalent—both have χ_{k}^{2} distributions where k is the number of parameters in M_{1} minus the number of parameters in M_{2}. Both measures are useful, but it appears that X^{2} is a little more accurate when the counts in all cells are quite small.

Another Example, from Raw data

I dug up the following very raw data from a data archive. I would like to study these data, finding interesting questions, and manipulating the data to answer the questions. The subjects were 478 students in grades 4–6 from three school districts in Ingham and Clinton Counties, Michigan. Chase and Dummer (Chase, M. A., and Dummer, G. M. (1992), The Role of Sports as a Social Determinant for Children, Research Quarterly for Exercise and Sport, 63, 418-424 ) stratified their sample, selecting students from urban, suburban, and rural school districts with approximately 1/3 of their sample coming from each district. Students indicated whether good grades, athletic ability, or popularity was most important to them. They also ranked four factors: grades, sports, looks, and money, in order of their importance for popularity. The questionnaire also asked for gender, grade level, and other demographic information.

The variable names are

1. Gender: Boy or girl

2. Grade: 4, 5 or 6

3. Age: Age in years

4. Race: White, Other

5. Urban/Rural: Rural, Suburban, or Urban school district

6. School: Brentwood Elementary, Brentwood Middle, Ridge, Sand, Eureka, Brown, Main, Portage, Westdale Middle

7. Goals: Student’s choice in the personal goals question where options were 1 = Make Good Grades, 2 = Be Popular, 3 = Be Good in Sports

8. Grades: Rank of ”make good grades” (1=most important for popularity, 4=least important)

9. Sports: Rank of ”being good at sports” (1=most important for popularity, 4=least important)

10. Looks: Rank of ”being handsome or pretty” (1=most important for popularity, 4=least important)

11. Money: Rank of ”having lots of money” (1=most important for popularity, 4=least important)

The first few lines of the data

Gender  Grade   Age     Race    Urban/Rural School  Goals   Grades  Sports  Looks   Money
boy     5       11      White   Rural       Elm     Sports  1       2       4       3
boy     5       10      White   Rural       Elm     Popular 2       1       4       3
girl    5       11      White   Rural       Elm     Popular 4       3       1       2

girl    5       11      White   Rural       Elm     Popular 2       3       4       1
girl    5       10      White   Rural       Elm     Popular 4       2       1       3
girl    5       11      White   Rural       Elm     Popular 4       2       1       3
girl    5       10      White   Rural       Elm     Popular 3       4       1       2
girl    5       10      White   Rural       Elm     Grades  3       4       2       1
girl    5       10      White   Rural       Elm     Sports  3       2       1       4
girl    5       10      White   Rural       Elm     Sports  4       3       2       1
girl    5       11      White   Rural       Elm     Sports  2       3       1       4
girl    4       10      White   Rural       Elm     Grades  2       3       4       1

So, now what do we do. This doesn’t look like a table at all. There seems to be one line of data for each subject. So, at the least we have to tabulate, but first we have to decide what we want tablulated.

There are many variables to consider, but for now I want to concentrate on just 5 of them, Gender, Grade, Race, Urban.Rural, and Goals. One of these (Urban.Rural) has to be considered an explanatory variable (recall the earlier discussion about how the data were collected). One can debate whether each of the others is explanatory or response, but there is at least no danger considering the rest as response variables.

Start by reading in the raw data......
> kids <- read.table("ex.dat", head=T)
> kids[1:10,]
   Gender Grade Age  Race Urban.Rural School   Goals Grades Sports Looks Money
 1    boy     5  11 White       Rural    Elm  Sports      1      2     4     3
 2    boy     5  10 White       Rural    Elm Popular      2      1     4     3
 3   girl     5  11 White       Rural    Elm Popular      4      3     1     2
 4   girl     5  11 White       Rural    Elm Popular      2      3     4     1
 5   girl     5  10 White       Rural    Elm Popular      4      2     1     3
 6   girl     5  11 White       Rural    Elm Popular      4      2     1     3
 7   girl     5  10 White       Rural    Elm Popular      3      4     1     2
 8   girl     5  10 White       Rural    Elm  Grades      3      4     2     1
 9   girl     5  10 White       Rural    Elm  Sports      3      2     1     4
10   girl     5  10 White       Rural    Elm  Sports      4      3     2     1
>

We still only have the raw data, and we have no good ideas on how to tabulate. There has to be a better way of doing this, but right now I can’t seem to figure it out. I’m open to suggestions. What I did was use the Splus tabulate command, and then use some tricks from there.

> attach(kids,pos=1)
> table(Gender,Grade,Race,Urban.Rural,Goals)
, , Other, Rural, Grades
     4 5 6
 boy 0 0 0

girl 0 1 2
, , White, Rural, Grades
      4 5  6
 boy  7 6  8
girl 12 9 12
, , Other, Suburban, Grades
     4 5 6
 boy 0 2 2
girl 0 1 0

This has clearly produced something that looks like a matrix (or table), but that is not the form that we need for the glm command. (It is what the loglinear command needs, but we won’t ever use that.

As it isn’t the format we need, I created a little bit of code to turn this into the proper form. If I were smarter, I would clean it up and make it more general, but I’m not.

doit <- function(){
        fofo <-table(Gender,Grade,Race,Urban.Rural,Goals)
        dimd <- attributes(fofo)$dim
        name <- attributes(fofo)$dimnames
        fnames<-list(Gender=name[[1]],Grade=name[[2]],
             Race=name[[3]],Urban=name[[4]],Goals=name[[5]])
        ans <- fac.design(dimd,fnames)
        ans$Resp <- as.vector(fofo)
        ans
}

It doesn’t much matter if you follow this or not, the chief idea is that it takes a table and converts it into the correct type of data frame.

> kids.f <- doit()
> kids.f[1:10,]
   Gender Grade  Race Urban  Goals Resp
 1    boy     4 Other Rural Grades    0
 2   girl     4 Other Rural Grades    0
 3    boy     5 Other Rural Grades    0
 4   girl     5 Other Rural Grades    1
 5    boy     6 Other Rural Grades    0
 6   girl     6 Other Rural Grades    2
 7    boy     4 White Rural Grades    7
 8   girl     4 White Rural Grades   12
 9    boy     5 White Rural Grades    6
10   girl     5 White Rural Grades    9

> assign("kids.f", kids.f,w=2)
> detach(1,save=F)

Whew!, now at least I have the data in a form that I can use, so maybe I can get on with the modelling. Now, try to look at the data as a loglinear model.

> attach(kids.f,pos=1)
> ls()
[1] ".Last.value" "Gender"      "Goals"       "Grade"       "Race"
[6] "Resp"        "Urban"

What follows is not a particularly directed analysis, but we will end up in a reasonable place.

> mymod <- glm(Resp ~ Urban*Gender*Goals*Grade*Race -
         Urban:Gender:Goals:Grade:Race, family=poisson)
Warning messages:
  linear convergence not obtained in  10  iterations. in: glm.fitter(x = X, y =
        Y, w = w, start = start, offset = offset, family = family, ....

Note the error messages about lack of convergence. Why might that happen, and why didn’t it worry me (yet).

> anova(mymod)
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
                         Df Deviance Resid. Df Resid. Dev
                    NULL                   107   679.9861
                   Urban  2   3.2349       105   676.7511
                  Gender  1   1.2055       104   675.5456
                   Goals  2  79.2794       102   596.2662
                   Grade  2  16.2392       100   580.0271
                    Race  1 407.2321        99   172.7949
            Urban:Gender  2  12.1266        97   160.6683
             Urban:Goals  4  18.5706        93   142.0976
            Gender:Goals  2  23.3472        91   118.7504
             Urban:Grade  4  30.9144        87    87.8360
            Gender:Grade  2   0.8923        85    86.9437
             Goals:Grade  4   1.3379        81    85.6058

              Urban:Race  2  10.2736        79    75.3322
             Gender:Race  1   0.3164        78    75.0158
              Goals:Race  2   0.6479        76    74.3679
              Grade:Race  2   1.0471        74    73.3208
      Urban:Gender:Goals  4   3.7229        70    69.5979
      Urban:Gender:Grade  4   4.3112        66    65.2867
       Urban:Goals:Grade  8   6.4401        58    58.8466
      Gender:Goals:Grade  4   5.2033        54    53.6433
       Urban:Gender:Race  2   1.2544        52    52.3889
        Urban:Goals:Race  4   4.9217        48    47.4672
       Gender:Goals:Race  2   2.5781        46    44.8892
        Urban:Grade:Race  4   4.5668        42    40.3224
       Gender:Grade:Race  2   1.2655        40    39.0568
        Goals:Grade:Race  4   4.9344        36    34.1224
Urban:Gender:Goals:Grade  8  10.4653        28    23.6572
 Urban:Gender:Goals:Race  4  11.2730        24    12.3841
 Urban:Gender:Grade:Race  4   2.2919        20    10.0922
  Urban:Goals:Grade:Race  8   3.5111        12     6.5811
 Gender:Goals:Grade:Race  4   6.5796         8     0.0015
Warning messages:
     (like before)

That was actually quite informative. (The model fits well, it might fit better if converged). The 4-way interactions all seem pretty useless, so lets remove them all.

> mymod <- update(mymod, . ~ . - Gender:Goals:Grade:Race -
              Urban:Goals:Grade:Race - Urban:Gender:Grade:Race -
              Urban:Gender:Goals:Race - Urban:Gender:Goals:Grade)
> anova(mymod)
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
                   Df Deviance Resid. Df Resid. Dev
              NULL                   107   679.9861
             Urban  2   3.2349       105   676.7511
            Gender  1   1.2055       104   675.5456
             Goals  2  79.2794       102   596.2662
             Grade  2  16.2392       100   580.0271
              Race  1 407.2321        99   172.7949
      Urban:Gender  2  12.1266        97   160.6683
       Urban:Goals  4  18.5706        93   142.0976

      Gender:Goals  2  23.3472        91   118.7504
       Urban:Grade  4  30.9144        87    87.8360
      Gender:Grade  2   0.8923        85    86.9437
       Goals:Grade  4   1.3379        81    85.6058
        Urban:Race  2  10.2736        79    75.3322
       Gender:Race  1   0.3164        78    75.0158
        Goals:Race  2   0.6479        76    74.3679
        Grade:Race  2   1.0471        74    73.3208
Urban:Gender:Goals  4   3.7229        70    69.5979
Urban:Gender:Grade  4   4.3112        66    65.2867
 Urban:Goals:Grade  8   6.4401        58    58.8466
Gender:Goals:Grade  4   5.2033        54    53.6433
 Urban:Gender:Race  2   1.2544        52    52.3889
  Urban:Goals:Race  4   4.9217        48    47.4672
 Gender:Goals:Race  2   2.5781        46    44.8892
  Urban:Grade:Race  4   4.5668        42    40.3224
 Gender:Grade:Race  2   1.2655        40    39.0568
  Goals:Grade:Race  4   4.9344        36    34.1224

Again, the high order interactions are quite pointless, so be cavalier and remove them all.

> mymod <- update(mymod, . ~
  (Urban+Gender+Goals+Grade+Race)*
  (Urban+Gender+Goals+Grade+Race))
> anova(mymod)
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
             Df Deviance Resid. Df Resid. Dev
        NULL                   107   679.9861
       Urban  2   3.2349       105   676.7511
      Gender  1   1.2055       104   675.5456
       Goals  2  79.2794       102   596.2662
       Grade  2  16.2392       100   580.0271
        Race  1 407.2321        99   172.7949
Urban:Gender  2  12.1266        97   160.6683
 Urban:Goals  4  18.5706        93   142.0976
 Urban:Grade  4  30.9144        89   111.1832
  Urban:Race  2  10.2736        87   100.9097
Gender:Goals  2  23.3472        85    77.5624

Gender:Grade  2   0.8923        83    76.6701
 Gender:Race  1   0.3164        82    76.3537
 Goals:Grade  4   1.3379        78    75.0158
  Goals:Race  2   0.6479        76    74.3679
  Grade:Race  2   1.0471        74    73.3208

There are still some two-way interactions that are very small, so remove them.

> anova(mymod <- update(mymod, . ~ . - Grade:Race - Goals:Race - Goals:Grade))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
             Df Deviance Resid. Df Resid. Dev
        NULL                   107   679.9861
       Urban  2   3.2349       105   676.7511
      Gender  1   1.2055       104   675.5456
       Goals  2  79.2794       102   596.2662
       Grade  2  16.2392       100   580.0271
        Race  1 407.2321        99   172.7949
Urban:Gender  2  12.1266        97   160.6683
 Urban:Goals  4  18.5706        93   142.0976
 Urban:Grade  4  30.9144        89   111.1832
  Urban:Race  2  10.2736        87   100.9097
Gender:Goals  2  23.3472        85    77.5624
Gender:Grade  2   0.8923        83    76.6701
 Gender:Race  1   0.3164        82    76.3537
>
> anova(mymod <- update(mymod, . ~ . - Gender:Race - Gender:Grade))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
             Df Deviance Resid. Df Resid. Dev
        NULL                   107   679.9861
       Urban  2   3.2349       105   676.7511
      Gender  1   1.2055       104   675.5456
       Goals  2  79.2794       102   596.2662
       Grade  2  16.2392       100   580.0271

        Race  1 407.2321        99   172.7949
Urban:Gender  2  12.1266        97   160.6683
 Urban:Goals  4  18.5706        93   142.0976
 Urban:Grade  4  30.9144        89   111.1832
  Urban:Race  2  10.2736        87   100.9097
Gender:Goals  2  23.3472        85    77.5624
> qchisq(0.95,85)
[1] 107.5217
> qchisq(0.95,2)
[1] 5.991465
> qchisq(0.95,4)
[1] 9.487729

OK, now just about everything looks significant, and most surprising the overall model fits quite well! For reasons that will be much clearer later, I prefer to add a term to this model—a term which by itself is not very significant.

> anova(mymod <- update(mymod, . ~ . + Urban*Gender*Goals))
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
                   Df Deviance Resid. Df Resid. Dev
              NULL                   107   679.9861
             Urban  2   3.2349       105   676.7511
            Gender  1   1.2055       104   675.5456
             Goals  2  79.2794       102   596.2662
             Grade  2  16.2392       100   580.0271
              Race  1 407.2321        99   172.7949
      Urban:Gender  2  12.1266        97   160.6683
       Urban:Goals  4  18.5706        93   142.0976
       Urban:Grade  4  30.9144        89   111.1832
        Urban:Race  2  10.2736        87   100.9097
      Gender:Goals  2  23.3472        85    77.5624
Urban:Gender:Goals  4   3.7406        81    73.8219
>

For this model the plots look OK, I can interpret the model and no other terms seem to be important.

But what is the interpretation of the model? I can only answer that by drawing a picture.

It isn’t obvious what this picture means, so let me “read” it to you. Remember that Urban is really an explanatory variable, but for now lets pretend it is a response. The picture says that Gender and Goals are conditionally independent of Race and Grade give the Urban Status. Even more, Gender and Goals are related (that is not independent), but Race is independent of Grade, given Urban Status.

Eventually you will learn how to read these sort of graphs, and you will also understand why I wanted the three-way interaction that I did.

More models and some notation

In order to interpret our results we intend to restrict our approach to a hierarchical family of models. The family we will consider restricts us to the class where if we want to fit an interaction term we must also fit the corresponding lower order interactions and main effects.

Now consider a 3-way table.

\begin{matrix}x_{1,1,1} & x_{1,1,2} & x_{1,1,3} \\ x_{1,2,1} & x_{1,2,2} & x_{1,2,3} \end{matrix}  \begin{matrix}x_{2,1,1} & x_{2,1,2} & x_{2,1,3} \\ x_{2,2,1} & x_{2,2,2} & x_{2,2,3} \end{matrix}

We will consider the following models, abbreviations, and corresponding pictures.

Model  Sufficient Stats Picture  Meaning

u + u_{1} + u_{2} + u_{3}  [1][2][3]  1 indep 2 indep 3

u + u_{1} + u_{2} + u_{3} + u_{12} [12][3] (1 and 2) indep 3

u + u_{1} + u_{2} + u_{3}+  [12][23]  1 indep 3 given 2

u_{12} + u_{13} u + u_{1} + u_{2} + u_{3}+ [12][23][23] none I don’t know u_{12} + u_{13} + u_{23}

The [·] are the indicate the margins of the table that make up the sufficient statistics— thus [1] is the margin formed by summing over the 2 and 3 indices of the table and [12] is the margin formed by summing over the 3 index.

Pictures (like the above) are really only appropriate when all the variables are responses, but they are very useful in all cases.

Model Selection

In order to select a model it is necessary that we be able to compare models. The techniques available allow us to compare, say, M_{5} with M_{2} but that does not necessarily allow for a comparison between M_{5} with M_{6}. This is exactly the same problem we face in linear regression—where we looked at things like RSS, residuals, and other properties to decide on a good model.

Models M_{5} with M_{2} can be compared using either Pearsons X^{2} or the log-likelihood G^{2} statistic. So can models M_{5} with M_{0}. For the G^{2} statistics (and not the X^{2} statistic) it is true that

G^{2}(M_{2}|M_{0}) + G^{2}(M_{5}|M_{2}) = G^{2}(M_{5}|M_{0})

i.e., an additive partition of χ^{2} with the two components being independent χ_{n1}^{2} and χ_{n2}^{2} components.

We will eventually discuss formal selection techniques (like forward selection and backward elimination), and less formal, but more useful, techniques based on the independence diagrams. Often the data will deal us to consider models which involve more information than just the classifying variables. In particular we will consider models where we have (or can construct) a vector(s) of covariates, which leads to ordered categories and other models.

Higher Dimensional Tables

Consider the following 3 × 2 × 2 × 2 contingency table where the classification variables are

1. Water softness (3 levels)

2. Previous use of Brand M laundry detergent

3. Temperature at which the laundry is done

4. Preference of brand M over brand X

and the counts are the number of people in each category. We are primarily interested in (4), but the way the data were collected makes each of (1)–(4) response variables. Thus we aim to model how the response variables interact.

Previous user of M  Previous nonuser of M
Water Brand High Low High Low
softness preference temperature temperature temperature temperature
Soft  X  19  57  29  63
M  29  49  27  53
Medium  X  23  47  33  66
M  47  55  23  50
Hard  X  24  37  42  68
M  43  52  30  42

IF we consider only models that can be written like factorial models (e.g., log(m_{i,j,k,l} = u + u_{1(i)} + ⋯) and include all main effects, there are still 113 possible models. Even if all these models had closed form estimates (which they don’t) look at all the models would be difficult. As always, our prefer to have a reasonable (fitting) model with a small number of parameters over a similar model with more parameters.

Fienberg, in the Xmas tree book gives the following scenario for constructing a model for these data. (This ignores the graphic model stuff we have just been talking about. I will try to do a graphical model in a moment.)

1. The simplest model we consider has no interactions; in the abbreviated notation it is [1][2][3][4]

2. It would be natural for brand preference to be related to use, so our next model includes this interaction, [24][1][3]

3. Since some detergents are designed for cold water and some for hot water, it would not be unreasonable for preference to be related to temperature as well, [24][34][1]

4. Next, we add the interaction between softness and temperature, since the softer the water, the higher the temperature the detergent manufacturers usually suggest, [12][24][34]

5. Finally we include some models that incorporate three-factor interactions, [13][234], and [123][234]

Lets start by getting in the data, and looking at it.

> ries <- fac.design(c(2,2,2,3),list(Temp=c("H","L"),
+ Pre=c("M","X"),Brand=c("X","M"), Soft=c("Soft","Med","Hard")))
> ries$Resp <- scan("ries-smith.dat")
> ries
   Temp Pre Brand Soft Resp
 1    H   M     X Soft   19
 2    L   M     X Soft   57
 3    H   X     X Soft   29
 4    L   X     X Soft   63
 5    H   M     M Soft   29
 6    L   M     M Soft   49
 7    H   X     M Soft   27
 8    L   X     M Soft   53
 9    H   M     X  Med   23
10    L   M     X  Med   47
11    H   X     X  Med   33
12    L   X     X  Med   66
13    H   M     M  Med   47
14    L   M     M  Med   55
15    H   X     M  Med   23
16    L   X     M  Med   50
17    H   M     X Hard   24
18    L   M     X Hard   37
19    H   X     X Hard   42
20    L   X     X Hard   68
21    H   M     M Hard   43
22    L   M     M Hard   52
23    H   X     M Hard   30
24    L   X     M Hard   42

Start by fitting the simplest model in the list above,

> glm(Resp ~ Temp+Pre+Brand+Soft, family=poisson, data=ries)
Call:
glm(formula = Resp ~ Temp + Pre + Brand + Soft, family = poisson, data = ries)
Coefficients:
 (Intercept)      Temp       Pre       Brand      Soft1       Soft2
    3.699209 0.2745538 0.0436785 -0.00793668 0.02687216 0.003092153
Degrees of Freedom: 24 Total; 18 Residual
Residual Deviance: 42.92866

Here is a list of the other models suggested above

[24][1][3]
glm(formula = Resp ~ Temp + Brand * Pre + Soft, family = poisson,
        data= ries)
Degrees of Freedom: 24 Total; 17 Residual
Residual Deviance: 22.34719
[24][34][1]
glm(formula = Resp ~ Soft + Pre * Brand + Temp * Brand, family = poisson,
         data = ries)
Degrees of Freedom: 24 Total; 16 Residual
Residual Deviance: 17.98559
[13][24][34]
glm(formula = Resp ~ Soft * Temp + Pre * Brand + Temp * Brand, family
        = poisson, data = ries)
Degrees of Freedom: 24 Total; 14 Residual
Residual Deviance: 11.88649
[13][234]
glm(formula = Resp ~ Soft * Temp + Pre * Brand * Temp, family = poisson,
      data = ries)
Degrees of Freedom: 24 Total; 12 Residual
Residual Deviance: 8.406872
[123][234]
glm(formula = Resp ~ Pre * Soft * Temp + Pre * Brand * Temp,
        family = poisson, data = ries)
Degrees of Freedom: 24 Total; 8 Residual
Residual Deviance: 5.656044

What was the good model?

In many respects this type of “guided” model selection is very satisfying. But, as with regression a more formal stepwise procedure is often helpful. Splus provides all the usual tools, so lets just let it go.

> step(mymod,scope=list(upper=~.^4,lower=~.^1))
Start:  AIC= 54.9287
 Resp ~ Pre + Soft + Temp + Brand
Single term additions
Model:
Resp ~ Pre + Soft + Temp + Brand

scale:  1
           Df Sum of Sq      RSS       Cp
    <none>              43.87621 55.87621
  Pre:Soft  2   1.07484 42.80136 58.80136
  Pre:Temp  1   1.25286 42.62335 56.62335
 Pre:Brand  1  20.50510 23.37111 37.37111
 Soft:Temp  2   6.07932 37.79688 53.79688
Soft:Brand  2   0.39516 43.48104 59.48104
Temp:Brand  1   4.35622 39.51998 53.51998
Step:  AIC= 36.3472
 Resp ~ Pre + Soft + Temp + Brand + Pre:Brand
Single term deletions
Model:
Resp ~ Pre + Soft + Temp + Brand + Pre:Brand
scale:  1
          Df Sum of Sq      RSS       Cp
   <none>              23.12993 37.12993
Pre:Brand  1  20.37267 43.50260 55.50260
Single term additions
Model:
Resp ~ Pre + Soft + Temp + Brand + Pre:Brand
scale:  1
           Df Sum of Sq      RSS       Cp
    <none>              23.12993 37.12993
  Pre:Soft  2  1.075196 22.05473 40.05473
  Pre:Temp  1  1.253319 21.87661 37.87661
 Soft:Temp  2  6.081331 17.04860 35.04860
Soft:Brand  2  0.395230 22.73470 40.73470
Temp:Brand  1  4.357797 18.77213 34.77213
Step:  AIC= 33.9856
 Resp ~ Pre + Soft + Temp + Brand + Pre:Brand + Temp:Brand
Single term deletions

Model:
Resp ~ Pre + Soft + Temp + Brand + Pre:Brand + Temp:Brand
scale:  1
           Df Sum of Sq      RSS       Cp
    <none>              18.33277 34.33277
 Pre:Brand  1  20.37196 38.70473 52.70473
Temp:Brand  1   4.35039 22.68316 36.68316
Single term additions
Model:
Resp ~ Pre + Soft + Temp + Brand + Pre:Brand + Temp:Brand
scale:  1
           Df Sum of Sq      RSS       Cp
    <none>              18.33277 34.33277
  Pre:Soft  2  1.075230 17.25754 37.25754
  Pre:Temp  1  0.691974 17.64079 35.64079
 Soft:Temp  2  6.081616 12.25115 32.25115
Soft:Brand  2  0.395246 17.93752 37.93752
Step:  AIC= 31.8865
 Resp ~ Pre + Soft + Temp + Brand + Pre:Brand + Temp:Brand + Soft:Temp
Single term deletions
Model:
Resp ~ Pre + Soft + Temp + Brand + Pre:Brand + Temp:Brand + Soft:Temp
scale:  1
           Df Sum of Sq      RSS       Cp
    <none>              11.91664 31.91664
 Pre:Brand  1  20.37144 32.28808 50.28808
Temp:Brand  1   4.35029 16.26693 34.26693
 Soft:Temp  2   6.05905 17.97569 33.97569
Single term additions
Model:
Resp ~ Pre + Soft + Temp + Brand + Pre:Brand + Temp:Brand + Soft:Temp
scale:  1

           Df Sum of Sq      RSS       Cp
    <none>              11.91664 31.91664
  Pre:Soft  2  1.088615 10.82802 34.82802
  Pre:Temp  1  0.691996 11.22464 33.22464
Soft:Brand  2  0.343602 11.57304 35.57304
Call:
glm(formula = Resp ~ Pre + Soft + Temp + Brand + Pre:Brand + Temp:Brand + Soft:
        Temp, family = poisson, data = ries)
Coefficients:
 (Intercept)        Pre      Soft1      Soft2      Temp      Brand  Pre:Brand
    3.682699 0.04343019 0.04342718 0.01564269 0.2782644 0.01658998 -0.1437655
  Temp:Brand   Soft1Temp   Soft2Temp
 -0.06836051 -0.05251834 -0.04906982
Degrees of Freedom: 24 Total; 14 Residual
Residual Deviance: 11.88649

We end up with a model that is at least familiar looking.

We could also proceed by looking at parameter values, again as we did in regression, but it is a little tricky.

>  mymod <- glm(Resp ~ Soft*Pre*Temp*Brand, family=poisson, data=ries)
> summary(mymod)
Call: glm(formula = Resp ~ Soft * Pre * Temp * Brand, family = poisson, data = ries)
Coefficients:
                         Value Std. Error       t value
      (Intercept)  3.674892823 0.03362662  109.28523105
            Soft1  0.034602276 0.04166801    0.83042779
            Soft2  0.012394833 0.02349486    0.52755510
              Pre  0.042297443 0.03362662    1.25785595
             Temp  0.285673485 0.03362662    8.49545670
            Brand  0.015236657 0.03362662    0.45311295
         Soft1Pre -0.039962253 0.04166801   -0.95906306
         Soft2Pre  0.016015575 0.02349486    0.68166297
        Soft1Temp -0.045744735 0.04166801   -1.09783816
        Soft2Temp -0.052760617 0.02349486   -2.24562393
         Pre:Temp  0.025849868 0.03362662    0.76873229
       Soft1Brand  0.012866861 0.04166801    0.30879469
       Soft2Brand -0.001045718 0.02349486   -0.04450837
        Pre:Brand -0.157008089 0.03362662   -4.66916072
       Temp:Brand -0.064072675 0.03362662   -1.90541531

     Soft1PreTemp  0.048167488 0.04166801    1.15598236
     Soft2PreTemp -0.000712269 0.02349486   -0.03031595
    Soft1PreBrand -0.062159775 0.04166801   -1.49178638
    Soft2PreBrand -0.030357357 0.02349486   -1.29208509
   Soft1TempBrand  0.012586592 0.04166801    0.30206845
   Soft2TempBrand  0.007774692 0.02349486    0.33091034
   Pre:Temp:Brand  0.050458667 0.03362662    1.50055726
Soft1PreTempBrand  0.010509090 0.04166801    0.25221001
Soft2PreTempBrand -0.019138419 0.02349486   -0.81457903
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 118.6269 on 23 degrees of freedom
Residual Deviance: 0 on 0 degrees of freedom

Quite clearly the 4 factor interaction, and all the three factor interactions can be removed, because of small t-values, which leaves us with all t-factor interactions.

>   mymod <- glm(Resp ~ (Soft+Pre+Temp+Brand)^2, family=poisson, data=ries)
> summary(mymod)
Call: glm(formula = Resp ~ (Soft + Pre + Temp + Brand)^2, family = poisson, data =
        ries)
Deviance Residuals:
       Min         1Q     Median        3Q      Max
 -1.107065 -0.4930669 0.02855437 0.5093409 1.386992
Coefficients:
                   Value Std. Error      t value
(Intercept)  3.682222962 0.03328281  110.6343930
      Soft1  0.043536098 0.04083575    1.0661272
      Soft2  0.015051992 0.02294280    0.6560661
        Pre  0.036165800 0.03309738    1.0927087
       Temp  0.277223691 0.03296423    8.4098346
      Brand  0.015077502 0.03315577    0.4547474
   Soft1Pre -0.023416451 0.03912050   -0.5985724
   Soft2Pre  0.017893207 0.02258248    0.7923490
  Soft1Temp -0.050438215 0.04091668   -1.2327056
  Soft2Temp -0.049901584 0.02296955   -2.1725101
 Soft1Brand  0.017642716 0.03916083    0.4505195
 Soft2Brand -0.002361722 0.02257536   -0.1046150
   Pre:Temp  0.028578975 0.03322482    0.8601695
  Pre:Brand -0.141752874 0.03193614   -4.4386345

 Temp:Brand -0.064157285 0.03321404   -1.9316312
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 118.6269 on 23 degrees of freedom
Residual Deviance: 9.846211 on 9 degrees of freedom

Now things get a little more interesting. It looks like we can get rid of Pre:Temp and Soft:Brand, but what else?

> summary(mymod <- update(mymod, . ~ . - Pre:Temp - Soft:Brand))
Call: glm(formula = Resp ~ Soft + Pre + Temp + Brand + Soft:Pre + Soft:Temp + Pre:
        Brand + Temp:Brand, family = poisson, data = ries)
Deviance Residuals:
       Min        1Q      Median        3Q      Max
 -1.207259 -0.395826 -0.02912609 0.4428775 1.363834
Coefficients:
                  Value Std. Error      t value
(Intercept)  3.68213446 0.03329713  110.5841383
      Soft1  0.04409441 0.04083153    1.0799107
      Soft2  0.01492096 0.02294642    0.6502525
        Pre  0.04385639 0.03187799    1.3757577
       Temp  0.27826378 0.03293986    8.4476310
      Brand  0.01662144 0.03311757    0.5018918
   Soft1Pre -0.02717637 0.03867930   -0.7026076
   Soft2Pre  0.01692800 0.02229341    0.7593277
  Soft1Temp -0.05225327 0.04080825   -1.2804584
  Soft2Temp -0.04923468 0.02290183   -2.1498144
  Pre:Brand -0.14381446 0.03185316   -4.5149199
 Temp:Brand -0.06847041 0.03277791   -2.0889194
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 118.6269 on 23 degrees of freedom
Residual Deviance: 10.79797 on 12 degrees of freedom

The only likely suspect is Soft:Pre, so remove it.

> summary(mymod <- update(mymod, . ~ . - Pre:Soft))
Call: glm(formula = Resp ~ Soft + Pre + Temp + Brand + Soft:Temp + Pre:Brand

    + Temp:Brand, family = poisson, data = ries)
Deviance Residuals:
       Min         1Q      Median        3Q      Max
 -1.433821 -0.3559708 -0.07358303 0.3974668 1.513481
Coefficients:
                  Value Std. Error      t value
(Intercept)  3.68269868 0.03327911  110.6609783
      Soft1  0.04342718 0.04080584    1.0642393
      Soft2  0.01564269 0.02290053    0.6830710
        Pre  0.04343019 0.03185253    1.3634772
       Temp  0.27826441 0.03293940    8.4477687
      Brand  0.01658998 0.03311667    0.5009554
  Soft1Temp -0.05251834 0.04080584   -1.2870299
  Soft2Temp -0.04906982 0.02290053   -2.1427374
  Pre:Brand -0.14376554 0.03185253   -4.5134735
 Temp:Brand -0.06836051 0.03277528   -2.0857337
(Dispersion Parameter for Poisson family taken to be 1 )
    Null Deviance: 118.6269 on 23 degrees of freedom
Residual Deviance: 11.88649 on 14 degrees of freedom

Not much obvious left to do. This is the same model as before.

What is this model, and what does it mean? What are the diagnostics?

The model is quite easy to read. Soft and Pre and independent of each other, given Temp and Brand. Soft is independent of Brand given Temp, and Pre is independent of Temp given Brand.

The residuals, and a plot, look pretty good.

Next topic, Fixed and Random Zeroes, then Graphical Models.

Fixed and Random Zeroes—A quick look

A zero in a contingency table can mean at least two different things.

1. A Sampling Zero, comes about when the cell in the table has positive probability but we didn’t observe anything there

2. A Fixed Zero, comes about when the cell in the cross-classification is not possible. then over the whole season we should get counts in every cell, whereas

It is easy to give examples of both types. If we did a two-way classification of professional baseball players by

Age  ≤ 30 > 30

Hit homer No homers

on the first day it is likely that the first row (hit homer) will be empty. (What sort of sampling would this be?)

Heart Cancer Old Age In Childbirth
Male  0

Female

For an example of a fixed, or structural zero, consider the following cross-classification of causes of death. There is one cell which obviously (at least to me) must be empty.

The two types of zeros give rise to different problems. If there are too many sampling zeroes, it is possible that the maximum likelihood estimates for a model may not exist. The problem is that there is not enough information to estimate the model. There is a theorem to make this precise.

Structural zeroes (i.e., fixed in advance) do not in themselves cause any estimation difficulties, but they may make it difficult to formulate a model. Reconsider the Sex × Death table and view both variables as responses. Lets make our initial model one that asks if cause of death is independent of sec. This implies that the table of fitted values must be

Heart Cancer Old Age In Childbirth
Male p_{1+}p_{+1} p_{1+}p_{+2} p_{1+}p_{+3}  p_{1+}p_{+4}
Female p_{2+}p_{+1} p_{2+}p_{+2} p_{2+}p_{+3}  p_{2+}p_{+4}

If the upper-right cell is non-zero, then we KNOW that there is something wrong. Clearly we need to develop new models for this situation. One solution to this (particular)

problem is to use the model of quasi-independence.

Let S be the set of cells which are net structural zeros. (E.g., above, S = {(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (2, 4)}. The quasi-independence model would posit that

log(m_{i,j}) = \{ \begin{matrix}u + u_{1(i)} + u_{2(j)} & for(i, j) ∈ S \\ 0 & otherwise \end{matrix}

This ensures that the structural zeroes “stay” zero and the rest of the table displays an independence-like structure.

There are rarely closed-form estiamtes for these models, but they are easy to fit. Just add weights to the glm that weight the missing cell as zero and the other cells as one.

As an example, consider the following data. A colony of 6 monkeys was studied over a period of time and a record was kept of how often each monkey displayed its genitals toward each other monkey. The constraint is that monkeys cannot (at least in this experiment) display to themselves. The data appear below.

Displayers  Watchers Name R S T U V W Total R  – 1 5 8 9 0  23 S  29 – 14 46 4 0  93 T  0 0 – 0 0 0  0 U  2 3 1 – 38 2  46 V  0 0 0 0 – 1  1 W  9 25 4 6 13 –  57 Total 40 29 24 60 64 3 220

Notice that T is bachful and never displays—hence we really need to treat the entire 3rd row as well as the diagonal as structural zeros.

Now lets try to fit the quasi-independence model to these data.

> monkey <- fac.design(c(6,6),list(Watch=c("R","S","T","U","V","W"),Display=c("R","S","T","U","V","W")))
> monkey$Resp <- scan("monkey.dat")
> monkey
   Watch Display Resp
 1     R       R    0
 2     S       R    1
 3     T       R    5
 4     U       R    8
 5     V       R    9
 6     W       R    0
 7     R       S   29
 8     S       S    0
 9     T       S   14
10     U       S   46
11     V       S    4
12     W       S    0
13     R       T    0
14     S       T    0
.........

Now set up the weights,

> wei <- rep(1,length(Resp))
> monkey[Display=="T",]
   Watch Display Resp
13     R       T    0
14     S       T    0
15     T       T    0
16     U       T    0
17     V       T    0
18     W       T    0
> monkey[Display==Watch,]
   Watch Display Resp
 1     R       R    0
 8     S       S    0
15     T       T    0
22     U       U    0
29     V       V    0
36     W       W    0
> wei[Display=="T"] <- 0
> wei[Display==Watch] <- 0
> wei
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
 0 1 1 1 1 1 1 0 1  1  1  1  0  0  0  0  0  0  1  1  1  0  1  1  1  1  1  1  0
 30 31 32 33 34 35 36
  1  1  1  1  1  1  0

Now, try to fit the model

> mymod <- glm(Resp ~ Watch + Display, family=poisson, weight=wei)
> anova(mymod)
Analysis of Deviance Table
Poisson model
Response: Resp
Terms added sequentially (first to last)
        Df Deviance Resid. Df Resid. Dev
   NULL                    35   352.9142
  Watch 16  95.6390        19   257.2753
Display  4 122.1061        15   135.1691

This clearly fits very poorly, the residuals and fitted values provide some insight.

> data.frame(Display, Watch, Resp, fitted(mymod), residuals(mymod))

   Display Watch Resp fitted.mymod. residuals.mymod.
 1       R     R    0    4.61117003       0.00000000
 2       R     S    1    5.25950848      -2.28011881
 3       R     T    5    2.48072581       1.40368049
 4       R     U    8    8.21594826      -0.07567288
 5       R     V    9    6.64804865       0.86505430
 6       R     W    0    0.39577027      -0.88968564
 7       S     R   29   19.18599258       2.08150826
 8       S     S    0   21.88357619       0.00000000
 9       S     T   14   10.32171588       1.08537395
10       S     U   46   34.18462584       1.91855873
11       S     V    4   27.66096481      -5.64376708
12       S     W    0    1.64670687      -1.81477650
13       T     R    0  394.66861611       0.00000000
14       T     S    0  450.15970329       0.00000000
15       T     T    0  212.32455423       0.00000000
16       T     U    0  703.20046854       0.00000000
17       T     V    0  569.00442629       0.00000000
18       T     W    0   33.87385456       0.00000000
19       U     R    2   10.93639583      -3.32821201
20       U     S    3   12.47407192      -3.22457813
21       U     T    1    5.88358252      -2.49456075
22       U     U    0   19.48591390       0.00000000
23       U     V   38   15.76729789       4.73158037
24       U     W    2    0.93865554       0.95032975
25       V     R    0    0.22001588      -0.66334890
26       V     S    0    0.25095050      -0.70844971
27       V     T    0    0.11836455      -0.48654816
28       V     U    0    0.39201311      -0.88545256
29       V     V    0    0.31720286       0.00000000
30       V     W    1    0.01888366       2.44472583
31       W     R    9    9.65764532      -0.21409227
32       W     S   25   11.01552689       3.60687599
33       W     T    4    5.19563795      -0.54687791
34       W     U    6   17.20750128      -3.12601524
35       W     V   13   13.92368867      -0.25035727
36       W     W    0    0.82890217       0.00000000

It appears that the Display is directed toward specific members (not randomly).

Graphical Models—An aid for Higher Dimensional Tables

We have already seen a number of attempts to introduce graphical models. For high dimensional tables (say 6 or higher) it is very difficult to develop models (let alone good models) without some form of extra guidance. If one restricts oneself the the class of graphical models then it is possible to draw useful pictures which can be a major help in selecting and interpreting models. The definition of what makes a graphical model is determined by a triangulation rule.

If A, B, and C are factors then the interactions A*B,A*C, and B*C imply that the A*B*C interaction must also be in the model.

Consider a model with 6 factors (A–F). Draw 6 points. For each two-way interaction in the model, add a line joining the points. For each factor not included in the model, change the filled circle to an unfilled one.

As an example, the model A*B*C + C*D*E has the following picture.

From the picture we can read off that D and E are independent of B and A given C. Can we draw a picture of A*B*C + C*D*E + A*E? A useful way to use these graphs is the choose models by adding and deleting edges to

the graphs.

Some examples, is A+B+C+D

and  is A*B + C*D

and  is A*B*D + C

and I don’t know the picture for A*B +C + C*A + D*B + A*B because of loss of triangularization.

Next topic, Survival Analysis.

Survival Analysis

References: Rupert G. Miller Jr., “Survival Analysis”, Wiley 1981.

John Kalbfleisch and R.L. Prentice, “The Statistical Analysis of Failure Time Data”, Wiley, 1980

Miller, Efron, et. al, “Biostatistics Casebook”, Wiley, 1980.

Venables and Ripley, Chapter 11.

Aim

For both of us (you and I, to learn something about the problems and methods of survival analysis. We will eventually take some data sets and attempt to do at least some rudimentary, and perhaps somewhat complicated, analyses on them.

Problem

We observer a positive valued random variable (the survivor time or failure time) and some co-variates. Typically the value of the random variable is the time to failure of some physical component or the time to death of a biological unit (or the time to relapse of some disease, etc.). There are all sorts of examples of waiting times which can be similarly analyzed—e.g., the length of time required to learn a new skill.

Early work in this area was some of the first (modern) statistics ever recorded. Graunt and Edmund Halley were responsible for the first useful life tables (estimate of the force of mortality) in about 166X. Medical and insurance applications carried the work steadily forward.

The modern era was ushered in by the need to test mechanical components for reliability. Rapid industrial expansion in the early 1940’s was a major impetus.

The most novel feature of failure-time data is that notion of censoring—not all sample units (e.g., patients in a medical trial) will be observed to failure (e.g., death), yet even the “partially” observed data provide useful information on the lifetime distribution.

Example Carcinogenesis. Two groups of rats were given pre-treatments then exposed to a carcinogen. We would like to compare the pre-treatments. Presumably, we could transform the data and perform a two-sample t-test, except there is some censoring in the data. Some of the rats died of causes unrelated to the cancer, or survived the study.

Days to Vaginal Cancer Mortality in Rats

Group 1 143 164 188 188 190 192 206 209 213 216
220 227 230 234 246 265 304 216* 244*
Group 2 142 156 163 198 205 232 232 233 233 233
233 239 240 261 280 280 296 296 323 204*

344*

The items with (*) are those that are censored. This is like a form of missing data except the “missing-ness” of the data contains information about the questions of interest.

Example Accelerated Life Test. Many physical components are designed to be very reliable and have long lifetimes. (and then there are the other components). It is still of interest to know something about the expected lifetime of a component or a comparison between differing designs. When the lifetime is very long, it is often difficult to actually test the component. (Think about a computer chip. If you tested the lifetime in the standard way, the chip would be obsolete before you saw any failures).

In these circumstances and “accelerated life test” is often performed. This consists of an experiment where the components are subjected to more stress than would be expected under normal operating conditions with the aim of extrapolating or comparing.

For a certain type of motor designed to operate ate 130 degrees C, the following severely censored data were obtained.

150C All 10 motorettes without failure at 8064 hours

170C 1764, 2772, 3444, 3542, 3780, 4860, 5196, 3 motorettes without failure at 5448 hours

190C 408, 408, 1344, 1344, 1440, 5 motorettes without failure at 1680 hours

220C 408, 408, 504, 504, 504, 5 motorettes without failure at 528 hours

How would you analyze these data? Suppose you wanted to compare two types of motor?

African Children A Stanford psychiatrist wanted to know the age at which a group of African children learned to perform a particular task—so he wrote a grant proposal and went to Africa. When he arrived in the village some of the children could already perform the task, some children learnt while he was there and some hadn’t learnt when he left. In this case there is both left and right censoring. For some individuals we know an age by which they could perform the task (left censoring); for some we know an age at which they could not perform the task, and for the rest we know the age at which they learnt the task.

Survival Functions and Hazard Rates

Let the random variable T > 0 denote survival time. Let T have density f(t) and distribution

F (t) = Prob(T < t)

(Note that some books use F (t) = Prob(T > t) which makes sense but is not standard).

The survival function is

S(t) = 1 − F (t) = P (T > t)

and the hazard rate (function) is

λ(t) = \frac{f(t)}{1 − F (t)}

which can be interpreted as

λ(t) = \frac{Prob(t < T < d + Δt|T > t)}{Δt}

which can be read as

Prob{dying in(t, t + Δt)|lived to time t}

or the instantaneous force of mortality.

Now, notice that if we integrate λ(t) we get

∫_{0}^{t} λ(u)du = ∫_{0}^{t} \frac{f(u)}{1 − F (u)}du
= − log [1 − F (u)]_{0}^{t}
= − log (S(t))

Thus

S(t) = exp {− ∫_{0}^{t} λ(u)du}

These various forms, definitions and ideas are very important. If you don’t understand them now, stop me.

Types of Censoring

Let T_{1}, …, T_{n} be iid F.

• Type I Censoring. Let t_{c} be some preassigned fixed number, the fixed censoring time. Instead of seeing T_{1}, …, T_{n} we get to see

Y_{i} = min (T_{i}, t_{c})

so the density of Y might look like

• Type II Censoring. Let r < n be fixed. All we observe are T (1), …, T (r).

Both types of censoring are common—either we have only a certain amount of time in which to experiment (I) or we don’t know how long we need to go but know how what percentage of failures we need to wait for (II).

• III Random Censoring. Let C_{1}, …, C_{n} be iid G were C_{i} is the censoring time associated

with T_{i}. We observe

Y_{i} = min(T_{i}, C_{i})

and

δ_{i} = I(T i ≤ C_{i}) = \{ \begin{matrix}1 & if T_{i} not censored \\ 0 & if T_{i} is censored \end{matrix}

Many medical examples have people enter the study at different times by and the exit the study because of

1. loss of follow up (they move, etc.)

2. drop out (side effects, contrariness)

3. end of study.

The following picture shows what might happen in a study.

One might also imaging left censoring, where the start date is only know to be after a certain date.

Some distributions for Survival Analysis

Exponential This is the basic model with constant risk, λ(t) = λ > 0. Thus

∫_{0}^{t} λ(u)du = λ × t

and so

S(t) = exp ∫_{0}^{t} λ(u)du = e^{−λt}
F (t) = 1 − e^{−λt}
f(t) = λe^{−λt}

and thus E(T) = 1/λ and Var(T) = 1/λ^{2}. Here is a plot of the density.

A reasonable thing to do might be to plot log (Ŝ(t)) against t. It should be a straight line.

Gamma One possible generalization of the exponential.

f(t) = \frac{λ^{α}}{Γ(α)}t^{α−1}e^{λt} α > 0 λ > 0

F (t) = Incomplete Gamma Integral, no closed form
λ(t) = 1 − e^{−λt}

E(T) = α/λ and Var(T) = α/λ^{2}. Here is a plot of the density, for α = 0.5, 1, 2, 3 You really already know a lot about the Gamma distribution because

χ_{k}^{2} = Gamma(k/2, 1/2)

that is, α = k/2 and λ = 0.5.

Weibull Another generalization of the exponential.

S(t) = e^{−(λt)^{α}} α > 0 λ > 0

This can be forced to look a LOT like a gamma distribution. The Weibull is very important, but we will not dwell on it.

With complete (i.e., uncensored) data estimation of any of these distributions is relatively easy. Just write out the likelihood equations and differentiate. For some of the models there are even closed-form estimates. We will defer looking at the parametric ML estimates until we have looked at a non-parametric estimate of the survivor function.

Estimation of the Survivor Function

(Without parametric models).

Remember what the survivor function is, just the probability of surviving to be older than time t. A convenient and familiar object is the sample cumulative distribution function (cdf), written as

F̄_{n}(x) = \frac{# of obs ≤ x}{n}

Estimates of the sample survival function are obviously of the form, 1 − F̄_{n}(x). In the analysis of survival data it is often useful to summarize the survival experience in terms of F̄_{n}(x) or S̄_{n}(x). (Why? Well it might direct us to an appropriate parametric model, or it might be interesting in its own right). When the data is censored we need to think of something other than the simple cdf.

Life Tables The classical method of estimating S(t) is the so-called actuarial method which depends on the life table.

Let time be partitioned into a fixed sequence of intervals, I_{1}, …, I_{k} mostly of equal length. The intervals should be open at the left and closed at the right.

Consider the line

with

n_{i} = # alive at beginning of I_{i}
d_{i} = # died in I_{i}

l_{i} = # lost to follow-up in I_{i}
w_{i} = # withdrawn in I_{i}
p_{i} = P ( surviving through I_{i}| alive at start of I_{i}

As an example, look at the data below, where we would like to estimate S(5) say.

Years after diag Alive at start Died Lost to Followup Withdrawn
0–1  126  47  4  15
1–2  60  5  6  11
2–3  38  2  —  15
3–4  21  2  2  7
4–5  10  —  —  6

Reduced Sample Method To estimate S(τ_{k}) we only use those subjects who are at risk for the entire period (0, τ_{k}]. That is then

\begin{aligned}n & = & n_{1} − ∑_{i=1}^{k} l_{i} − ∑_{i=1}^{k} w_{i} \\ d & = & ∑_{i=1}^{k} d_{i} \\ and S(τ_{k}) & = & 1 − d/n \end{aligned}  \begin{aligned}n & = & 126 − 12 − 54 = 60 \\ d & = & 56 \\ & = & 1 − \frac{56}{60} = 0.078 \end{aligned}

What (if anything) is wrong with this? It throws away information, all the iformation in l_{i} and w_{i}, and it is thus an inefficient (and biased) estimate.

Actuarial Method The survival probability S(τ_{k}) can be broken into a product of probabilities,

S(τ_{k}) = P (T > τ_{k})
= P (T > τ_{1}) · P (T > τ_{1}|T > τ_{2}) ⋯ P (T > τ_{k}|T > τ_{k−1})
= p_{1} · p_{2} · p_{3} ⋯ p_{k}

where

p_{i} = P (T > τ_{i}|T > τ_{i−1})

The actuarial method estimates each p_{i} separately and multiplies to get S(τ_{k}). To estimate p_{i} one could use

1 − \frac{d_{i}}{n_{i}}

, but this ignores losses and withdrawals in interval I_{i}. A sensible(?) approach might be to assume that on average those that were lost or withdrew were at risk for 1/2 of the interval. Thus one might define the effective sample size as

n′_{i} = n_{i} − 1/2 · (l_{i} − w_{i})

with

q̂′ = \frac{d_{i}}{n′_{i}}
p̂_{i} = 1 − q_{i}
and Ŝ(τ_{k}) = ∏_{i=1}^{k} (p̂_{i})

For the above data, this yields,

Ŝ(5 years) = 0.44

(You work this out!) Variance of Ŝ(τ_{k}) By assuming that

n′_{i}p̂_{i} ≈ Binomial(n′_{i}, p_{i})

and some manipulations an estimate of the variance of the estimate can be found as

V̂ar(Ŝ(τ_{k})) = Ŝ^{2}(τ_{k}) ∑_{i=1}^{k} \frac{d_{i}}{n′_{i}(n′_{i} − d_{i})}

This is known as Greenwood’s formula (named after Major Greenwood).

Next topic, Kaplan-Meier Product Limit Estimate

The Kaplan-Meier Product Limit Estimate

This is very similar to the actuarial method we just saw, except that the intervals I_{i} are not fixed in advance. In fact the intervals are determined by when an event occurs. So a revised picture might be,

Recall that we observe Y_{1}, …, Y_{n} survival times, and δ_{1}, …, δ_{n} ( \begin{matrix}1 & censored \\ 0 & not censored \end{matrix} )

Order the Y_{i} into \begin{matrix}Y_{(1)}, …, Y_{(n)} \\ δ_{(1)}, …, δ_{(n)} \end{matrix} and carry only the δ_{i}

Let R(t) be the risk set at t and

n_{i} = # in R(Y_{(i)}) = # alive just before Y_{(i)}
d_{i} = #died at Y_{(i)}
p_{i} and q_{i} = same as before—conditional probability of surviving

Then

p̂_{i} = \{ \begin{matrix}1 − 1/n_{i} & if & δ_{(i)} = 1 \\ 1 & if & δ_{(i)} = 0 \end{matrix}

Thus, without ties, the PL estimate is

Ŝ(t) = ∏_{y(i)≤t} p̂_{i} = ∏ (1 − 1/n_{i})^{δ(i)}
= ∏ (1 − \frac{1}{(n − i + 1)})^{δ(i)}

Notes

1. Tied Uncensored observations—do the obvious thing which is p_{i} = 1 − \frac{d}{m} where d deaths occur.

2. Censored and uncensored tied data—consider the uncensored data to be just before the censored data. (Why in this order?)

Two Views of the PL estimate as an MLE

The PL estimate would not be nearly so popular if it could not be interpreted in other ways.

Let t_{1}, …, t_{k} be the observed failure times in a sample of size n_{0} with survivor function, S(t). Suppose d_{j} items fail at time t_{j} and m_{j} items are censored in [t_{j}, t_{j+1}) at known times t_{j_{l}}.

Let n_{j} = (m_{j} + d_{j}) + ⋯ + (m_{k} + d_{k}) be the number of items at risk at time just before t_{j}. The probability of failure at t_{j} is

S(t_{j}) − S(t_{j} + 0)

We assume that the contribution to the likelihood of a survival time censored at t_{j_{l}} is P (T > t_{j_{l}}) = S(t_{j_{l}} + 0).

Thus,

L = ∏_{j=0}^{k} {[S(t_{j}) − S(t_{j} + 0)]^{d_{j}} · ∏_{l=1}^{m_{j}} S(t_{j_{l}} + 0)}

Given the data, this can be viewed as a likelihood function on the space of all survivor functions, S(t). There are problems with this definition of the MLE, but the intuition is good.

Clearly to maximize L one needs and S which is discontinuous at the observed failure times—otherwise L is zero.

Further, L is maximized when S(t_{j_{l}} + 0) = S(t_{j} + 0)

An estimate Ŝ(t) is then a discrete survivor function with hazard components λ̂_{i} at t_{i}. Thus,

Ŝ(t_{j}) = ∏_{l=1}^{j−1}(1 − λ̂_{l})

and

Ŝ(t_{j} + 0) = ∏_{l=1}^{j}(1 − λ̂_{l})

where the λ are chosen to maximize the function

∏_{j=1}^{k} {[∏_{l=1}^{j−1}(1 − λ_{j})^{d_{j}} ] λ_{j}^{d_{j}} [∏_{l=1}^{j}(1 − λ_{l})^{m_{j}} ]}

= ∏_{j=1}^{k} λ_{j}^{d_{j}}(1 − λ_{j})^{n_{j}−d_{j}}

because

S(t_{j}) − S(t_{j} + 0) = ∏_{l=1}^{j−1} λ_{j}(1 − λ_{j})

Obviously the maximum is achieved when λ̂_{j} = \frac{d_{j}}{n_{j}} and the product limit estimate is

F̂(t) = ∏_{{j \atop t_{j} <t}} (\frac{n_{j} − d_{j}}{n_{j}} )

Interpretation We are in effect making the conditional probability of failure at each t_{j} agree exactly with the observed conditional relative frequency of failure at t_{j}.

Problem: If we assume that S(t) is continuous then Ŝ(t) is always bumpy, which might give you cause for concern. (This can all be fixed by appropriate technical mumbling).

Example

No Censoring Assume that n = 5 and x = death and o = censor with

\begin{aligned}& x & x & x & x & x \\ p_{i} & (1 − \frac{1}{5}) & (1 − \frac{1}{4}) & (1 − \frac{1}{3}) & (1 − \frac{1}{2}) & (1 − \frac{1}{1}) \end{aligned}

where p_{i} is the p(survive to age i|alive before timei)

\begin{aligned}S_{i} & = & \frac{4}{5} & \frac{4}{5} · \frac{3}{4} & \frac{4}{5} · \frac{3}{4} · \frac{2}{3} & … \\ & = & \frac{4}{5} & \frac{3}{5} & \frac{2}{5} & \frac{1}{5} & 0 \end{aligned}

Thus, without censoring the PL estimate is the obvious one

Some censored data Let n = 7 and take the data

\begin{aligned}& x & x & o & x & o & x \\ p_{i} & \frac{6}{7} & \frac{5}{6} & & \frac{3}{4} & & \frac{1}{2} \end{aligned}

The PL estimate is

\begin{aligned}S_{i} & = & \frac{6}{7} & \frac{6}{7} · \frac{5}{6} & \frac{6}{7} · \frac{5}{6} & \frac{6}{7} · \frac{5}{6} · \frac{3}{4} & … \\ & = & \frac{6}{7} & \frac{5}{7} & \frac{5}{7} & \frac{15}{28} & \frac{15}{28} & \frac{15}{56} & \frac{15}{56} \end{aligned}

How does one compare two survival curves?

That is a very good question. One answer is to say, just use the Splus function survdiff, which has a description of

Test if there is a difference between two or more survival curves using the G-rho family of tests. Or, test for a single curve against a known alternative.

But, to really understand what this is all about, we need a digression.

Digression on some Nonparametric Techniques

Wilcoxen or Mann-Whitney Test

Suppose we have samples from two distributions F and G and we wish to compare the distributions. Let X_{1}, … X_{m} and Y_{1}, … Y_{m} denote the random samples from the two distributions.

One test of the difference between the distributions proceeds as 1. Rank all the data—i.e., mix the two samples.

2. Let R(X_{i}) denote the ranks of the F sample in the combined data.

3. Form S = ∑_{i=1}^{m} R(X_{i}), the rank sum.

Now note that S will be large if the F-distribution has larger values than the G-distribution and vice-versa. It is easy to show that

E(S) = \frac{m(n + m + 1)}{2}

and

Var(S) = \frac{m(n + m + 1)n}{2}

(This comes from assuming that F = G in which case the ranks are just sampled form a discrete uniform distribution.

For small values of n and m there are tables of exact tail probabilities for S (or some linear transformation of it).

For large samples, just take a normal approximation

z = \frac{S − E(S)}{\sqrt{Var(S)}}

The test makes minimal assumptions and works quite well in certain situations? However, what is the test sensitive to? When normality can be assumed, the 2-sample t-test is a competitor and like most specific purpose tools does better than the Wilcoxon test.

You can find this test as wilcox.test in Splus.

Kolmogorov-Smirnov test

The previous test is not very sensitive to differences in the distribution shape, as opposed to solely location differences. A test which is sometimes useful here is the (Kolmogorov- )Smirnov test.

Let F̂ and Ĝ be estimates of two distribution functions. (For censored data one might think about using the Kaplan-Meier estimates). The exact form of the test depends on the alternate hypothesis (i.e., what do you want to test) so I will assume something like H_{1} : F ≠ G.

Now form T = max |F − G|, and look up the results in a table, OR just use the Splus function ks.gof.

The test is sensitive to all sorts of departures, particularly changes in tails and different variances.

End Digression

So how does one compare two Kaplan-Meier estimates? One of the things that survdiff does is a (version) of Gehan’s test, which itself is just a modification of the Wilcoxon test.

Instead of using S, the rank-sum, we need to develop another equivalent statistic. Let

U(X_{i}, Y_{j}) = \{ \begin{matrix}+1 & if & X_{i} > Y_{j} \\ 0 & if & X_{i} = Y_{j} \\ −1 & if & X_{i} < Y_{j} \end{matrix}

Let U = ∑_{i=1}^{m} ∑_{j=1}^{n} U_{i,j}, which can be read as the average number of times that X > Y. (Would ∑_{i=1}^{m} ∑_{j=1}^{n} |U_{i,j}| make sense? ).

It isn’t hard to see that

S = \frac{m(m + n + 1)}{2} = \frac{1}{2}U

thus E(U) = 0 and Var(U) = \frac{nm(n+m+1)}{3}.

For censored data, the Gehan test suggests defining U_{i,j} in a different form. Our basic data is (X_{i}, δ_{i}) and (Y_{i}, ϵ_{i}) where δ and ϵ are the indicators of censoring. The null hypothesis is F_{1} = F_{2} (distribution) and G_{1} = G_{2} (censoring), with alternate

interest in F_{1} ≠ F_{2}. (I am not sure what happens when the censoring is different). Now define

U_{i,j} = \{ \begin{matrix}+1 & if & we know that t_{i} > u_{j} i.e., (x_{i} > y_{j}, ϵ_{j} = 1) or (x_{i} = y_{j}, δ_{i} = 0, ϵ_{i} = 1) \\ 0 & otherwise \\ −1 & if & we know that t_{i} < u_{j} i.e., (x_{i} < y_{j}, δ_{j} = 1) or (x_{i} = y_{j}, δ_{i} = 1, ϵ_{i} = 0) \end{matrix}

and U = ∑_{i=1}^{m} ∑_{j=1}^{n} U_{i,j}

It is possible to show (via some permutation theory) that E(U) = 0 and Var(U) = a big mess, but there is a computational formula, and packages.

Digression on the Mantel-Haenszel Test

Suppose there are two populations and each individual in each population can have one of two characteristics. For example, two populations of cancer patients with different treatments who can either survive beyond a year or die within the year. We could write out the table of data as,

Dead Alive
Pop1 a  b  n_{1}
Pop2 c  d  n_{2}
m_{1}  m_{2} n

Let

p_{1} = P (Dead|Pop1)
p_{2} = P (Dead|Pop2)

A test of p_{1} = p_{2} can be based on (thing that we already know about homogeneity of proportions) or the statistic

χ^{2} = [\frac{p̂_{1} − p̂_{2}}{\sqrt{p̂(1 − p̂)(\frac{1}{n_{1}} + \frac{1}{n_{2}} )}}]^{2} = \frac{n(ad − bc)^{2}}{n_{1}n_{2}m_{1}m_{2}}

where p̂_{1} = \frac{a}{n_{1}} , p̂_{2} = \frac{c}{n_{2}} , and p̂ = \frac{m_{1}}{n}. The χ^{2} is approximately distributed as a χ_{1}^{2} random variable.

The intuition for this test comes (again) from the 2-sample t-test, but the test itself is derived from a hypergeometric distribution for the contingency table.

P (A = a) = \begin{matrix}\binom{n_{1}}{a} & \binom{n_{2}}{m_{1}−a} \\ ( & {n \atop m_{1}} & ) \end{matrix}

This test assume that the n and m margins are BOTH fixed in advance.

Under the hypergeometric distribution (which is our null)

E(A) = \frac{n_{1}m_{1}}{n}

and

Var(A) = \frac{n_{1}n_{2}m_{1}m_{2}}{n^{2}(n − 1)}

whence

ad − bc = n(a − E(A)

and

n_{1}n_{2}m_{1}m_{2} = n^{2}(n − 1)Var(A)

thus

χ^{2} = \frac{n}{n − 1} [\frac{a − E(A)}{\sqrt{Var(A)}}]^{2}

For the example above, the hypergeometric assumption is suspect. If I collected the data by finding m_{1} deaths and m_{2} survivors and then recorded their populations, then it is probably OK. On the other hand, if I just start out with two populations and watch for deaths a binomial (or product binomial) seem more reasonable. But the particular form of the χ^{2} test is (I think) just the same as the binomial one would have been, modulo a fudge about n or n − 1.

The real reason that the hypergeometric is used here is based on the statistical notion of ancilarity

Continue Digression, a sequence of 2 × 2 tables

The same setup as before, but now there are K hospitals. Because there may be differences between hospitals we don’t wish to combine all the 2 × 2 tables. Rather we would like to test

p_{1}^{1} = p_{2}^{1}, p_{1}^{2} = p_{2}^{2}, …, p_{1}^{k} = p_{2}^{k},

i.e., equality in each population, but not overall equality.

The data are now

Dead Alive Pop1 a_{1}  b_{1}  n_{1}1 Pop2 c_{1}  d_{1}  n_{1}2 m_{1}1 m_{1}2 n_{1}  …  Dead Alive Pop1 a_{k}  b_{k}  n_{k}1 Pop2 c_{k}  d_{k}  n_{k}2 m_{k}1 m_{k}2 n_{k}

The Mantel-Haenszel statistic for this problem is

MH = \frac{∑_{i=1}^{k}(a_{i} − E(A_{i}))}{\sqrt{∑_{i=1}^{k}(Var(A_{i}))}}

If the tables are independent then MH tends towards a N(0, 1) when either k is fixed and n_{i} → ∞ or k → ∞ and the tables are iid.

Of course we could use our contingency table techniques to answer these questions—but you end up with slightly different tests.

When used in survival analysis, neither of these conditions for asymptotic normality holds—but the result still does.

End Digression

Kaplan-Meier in Survival Analysis

In survival analysis this test is used in the following manner. Let \begin{matrix}X_{1}, …, X_{n} \\ δ_{1}, …, δ_{n} \end{matrix} be one

population and \begin{matrix}Y_{1}, …, Y_{n} \\ ϵ_{1}, …, ϵ_{n} \end{matrix} be the other. Consider Z_{(1)}, …, Z_{(n+m)} to be the combined

ordered populations. For each Z_{(i)} which corresponds to a death (in either population) construct a 2 × 2 table and compute the MH test for this sequence of tables. (E.G., at time (i) we have

Dead Alive
Pop1 a  b  n_{1}
Pop2 c  d  n_{2}

and at time (i+1) we have

Dead Alive
Pop1 a + 1 b − 1 n_{1}
Pop2 c  d  n_{2}

etc.,)

In other words, test F_{1} = F_{2} (which is the same as S_{1} = S_{2}) by using the MH test. (Independence of tables is violated here, but asymptotic normality of the test can still be pushed through).

What next

Before proceeding to other (maybe more complicated) matters, we need to ask if any of these ideas are useful or used in practice. The best I can do is present a section of the book “Biostatistics Casebook” by Miller, Efron, et. al. which shows a very simple yet convincing application.

Next topic, Spring Break

Some handwritten material was inserted

Some examples, from Venables and Ripley

Look at the leukaemia data, which is the survival times and white blood counts for a group of leukaemia patients. The data has 33 leukaemia patients, and variables.

wbc:   white blood count
ag:    a test result, present or absent
time:     survival time in weeks

The Sources are D.R. Cox and D. Oakes (1984) Analysis of Survival Data, Chapman & Hall, p. 9, who themselves took the data from P. Feigl & M. Zelen (1965) Estimation of exponential survival probabilities with concomitant information. Biometrics, 21, 826-838.

The data are,

> leuk
      wbc      ag time
 1   2300 present   65
 2    750 present  156
 3   4300 present  100
 4   2600 present  134
 5   6000 present   16
 6  10500 present  108
 7  10000 present  121
 8  17000 present    4
 9   5400 present   39
10   7000 present  143
11   9400 present   56
12  32000 present   26
13  35000 present   22
14 100000 present    1
15 100000 present    1
16  52000 present    5
17 100000 present   65
18   4400  absent   56
19   3000  absent   65
20   4000  absent   17
21   1500  absent    7
22   9000  absent   16
23   5300  absent   22
24  10000  absent    3
25  19000  absent    4
26  27000  absent    2
27  28000  absent    3
28  31000  absent    8

29  26000  absent    4
30  21000  absent    3
31  79000  absent   30
32 100000  absent    4
33 100000  absent   43

The first thing one might wish to do with these data is plot the two survival curves, via

> attach(leuk)
> plot(survfit(Surv(time) ~ ag), lty=c(2,3))

where the survfit function, omputes an estimate of a survival curve for censored data using either the Kaplan-Meier or the Fleming-Herrington method or computes the predicted survivor function for a Cox proportional hazards model. In our case the model is just the two levels of ag. The picture looks like,

There are clearly differences in the curves, but if they are significant or not is not obvious. See V & R (page 269) for some suggestions on how to test this.

There are now many ways we could proceed with these data, including

• Parametric models based on Gamma errors and the glm function,

• Parametric models, based on the exponential and Weibull distributions using survreg, and • The family of Cox models.

We will not do a lot of the theory, but will attempt to push the example through with just enough to make sense of what we are doing.

The next example

lwbc <- log(wbc)
par(mfrow=c(2,1))
plot(log(time) ~ ag + log(wbc))
options(contrasts=c("contr.treatment","contr.poly"))
leuk.glm <- glm(time ~ ag*lwbc, Gamma(log))
 print("## summary(leuk.glm, dispersion=1)")
print(summary(leuk.glm, dispersion=1)) # force exponential data summary
print("## anova(leuk.glm))")
print(anova(leuk.glm))
# drop interaction term
leuk.glm <- update(leuk.glm,  ~ .-ag:lwbc)
print("## summary(leuk.glm, dispersion=1))")
print(summary(leuk.glm, dispersion=1))
# compare with natural link
leuk.glmi <- glm(time ~ ag*wbc, Gamma(inverse))
print("## summary(leuk.glmi, dispersion=1))")
print(summary(leuk.glmi, dispersion=1))
print("## anova(leuk.glmi))")
print(anova(leuk.glmi))

[1] "## summary(leuk.glm, dispersion=1)"
Call: glm(formula = time ~ ag * lwbc, family = Gamma(log))
Deviance Residuals:
       Min       1Q     Median        3Q      Max
 -1.992227 -1.21155 -0.3270025 0.2160396 1.564344
Coefficients:
                 Value Std. Error   t value
(Intercept)  4.3291270  1.9661806  2.201795
         ag  4.1481711  2.5702903  1.613892
       lwbc -0.1525468  0.2026831 -0.752637
    ag:lwbc -0.3291857  0.2668885 -1.233420
(Dispersion Parameter for Gamma family taken to be 1 )
    Null Deviance: 58.13839 on 32 degrees of freedom
Residual Deviance: 38.55468 on 29 degrees of freedom

Number of Fisher Scoring Iterations: 6
Correlation of Coefficients:
        (Intercept)         ag       lwbc
     ag -0.7649644
   lwbc -0.9918835   0.7587556
ag:lwbc  0.7532658  -0.9907268 -0.7594297
[1] "## anova(leuk.glm))"
Analysis of Deviance Table
Gamma model
Response: time
Terms added sequentially (first to last)
        Df Deviance Resid. Df Resid. Dev
   NULL                    32   58.13839
     ag  1 11.94009        31   46.19829
   lwbc  1  5.87920        30   40.31909
ag:lwbc  1  1.76441        29   38.55468
[1] "## summary(leuk.glm, dispersion=1))"
Call: glm(formula = time ~ ag + lwbc, family = Gamma(log))
Deviance Residuals:
       Min        1Q     Median        3Q      Max
 -2.174855 -1.266119 -0.4246385 0.4961872 1.904819
Coefficients:
                 Value Std. Error   t value
(Intercept)  5.8134506  1.2931888  4.495438
         ag  1.0182978  0.3492227  2.915898
       lwbc -0.3042296  0.1318634 -2.307157
(Dispersion Parameter for Gamma family taken to be 1 )
    Null Deviance: 58.13839 on 32 degrees of freedom
Residual Deviance: 40.31909 on 30 degrees of freedom
Number of Fisher Scoring Iterations: 6
Correlation of Coefficients:
     (Intercept)         ag
  ag -0.2090766

lwbc -0.9811356   0.0720421
[1] "## summary(leuk.glmi, dispersion=1))"
Call: glm(formula = time ~ ag * wbc, family = Gamma(inverse))
Deviance Residuals:
       Min        1Q     Median        3Q      Max
 -1.890629 -1.299155 -0.4143898 0.3559901 1.909283
Coefficients:
                    Value   Std. Error     t value
(Intercept)  5.635111e-02 1.887885e-02  2.98488018
         ag -4.780969e-02 1.932473e-02 -2.47401586
        wbc -2.030302e-08 4.239835e-07 -0.04788634
     ag:wbc  5.926747e-07 5.232047e-07  1.13277789
(Dispersion Parameter for Gamma family taken to be 1 )
    Null Deviance: 58.13839 on 32 degrees of freedom
Residual Deviance: 39.73216 on 29 degrees of freedom
Number of Fisher Scoring Iterations: 6
Correlation of Coefficients:
       (Intercept)         ag        wbc
    ag -0.9769269
   wbc -0.6746535   0.6590872
ag:wbc  0.5467114  -0.6022141 -0.8103588
[1] "## anova(leuk.glmi))"
Analysis of Deviance Table
Gamma model
Response: time
Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev
  NULL                    32   58.13839
    ag  1 11.94009        31   46.19829
   wbc  1  5.31608        30   40.88221
ag:wbc  1  1.15005        29   39.73216

Logistic Regression—Again

The Challenger Data

You are probably all familiar with the Challenger Space Shuttle accident of January 20, 1986. I will leave out all the interesting substantive details about the accident and focus instead on the data about O-rings that was available at the time of the accident.

Well, I do need a little background, so here goes. The space shuttle uses two booster rockets and each rocket is made of several pieces whose joints are sealed by rubber O-rings. The O-rings are supposed to stop the hot gases that are inside the booster from escaping. Each booster contains 3 primary O-rings. There was data about the O-rings on 23 previous flights of space shuttles. Here is that data.

Temp Fail          Temp Fail
53   2             70   1
57   1             70   1
58   1             72   0
63   1             73   0
66   0             75   0
67   0             75   2
67   0             76   0
67   0             76   0
68   0             78   0
69   0             79   0
70   0             81   0
70   0

On the day of the fated flight, the forecast launch temperature was 31^{∘}F. The O-rings were known to be sensitive to temperature. The previous coldest launch temperature was 53^{∘}. There was discussion about the cold temperature before launch, and the engineers did look at the data, but only the data for the failures. They produced a plot that looked something like this.

> chal <- read.table("chal.dat", head=T)
> attach(chal, pos=1)
> ls()
[1] ".Last.value" "Fail"        "Temp"
> n <- rep(6, length(Fail))
> plot(Temp[Fail >0], Fail[Fail>0])

Based on the erroneous assumption that the flights with no failures did not add any information, they noted that there did not appear to be any relationship between failures and temperature.

Had they produced a plot of the full data, which follows, they might have behaved differently.

> plot(Temp, Fail)

In addition before each flight the O-rings were given a pressure test. The pressure at which that test was conducted increased over time, and it is possible that the higher pressure test could have caused more failures. That data, and a lot more discussion, is available in the paper Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure, by S.R Dalal, E.B. Fowlkes, and B. Hoadley, in V. 84, pp.945–957, of the Journal of the American Statistical Association (JASA). You should be able to read most of that article. It is well written and easy to follow.

How to look at these data

Obviously, we are most interested how probability of failure relates to temperature, and other quantities known before a flight. Given our data, it is not immediately obvious what response variable we should be looking at. For example, should we look at the number of failures of O-rings in a flight, or just the fact that an O-ring failed? These might lead to different results (but fortunately for us they can be looked at using the same analyses).

We have in the past looked at models for count data that were of the form

log m = some combination of factor effects

Although these loglinear models did not require discrete (categorical) predictors, the only examples we saw were of that form. The current data can also be modelled using the standard loglinear model, but an equivalent and more natural approach is to set up the model as a logistic regression.

In this case, allow p(t) to be the probability of failure given the temperature t. If one wanted to be really careful, one might write this as

p(F |t) = probability of failure given temperature was t

but the shorthand notation will probably suffice.

A first approach is to consider that the failures are a binomial process with the probability of failure dependent on t, but with the events themselves independent. Thus using the first form of the data (where we explicitly model the number of O-rings that failed) we are assuming that

• The 6 O-rings on any flight are fail independently (but with a common probability), and

• The 23 flights are independent.

The second assumption is more tenable than the first, but even the first seems quite reasonable. In order to fret about the first assumption, we will also analyze a version of the data that just records whether there was a failure or not.

The Logistic model

Clearly, we can not get away with modelling

p(t) = some arbitrary function of t

unless we are very careful about the arbitrary function. One would quickly find oneself in the position where the models produced negative probabilities, or probabilities greater than 1.

One might think of the model (based on our loglinear experience)

log p(t) = some arbitrary function of t

but it suffers from the same problems.

Lets go back to first principles and think about the data we are trying to model.

We believe (at least I am suggesting) that the data can be modelled as follows. Let n be the number of trials that are done for each run of the experiment (in the Challenger case this is either 6 if each O-ring is considered or 1 if just overall success/failure is considered) and let r be the number of failures in those n trials. Let R be the corresponding random variable. Then we believe that

P (R = r|t) = \binom{n}{r}p(t)^{r}(1 − p(t))^{n−r}
= \binom{n}{r}(\frac{p(t)}{1 − p(t)})^{r}(1 − p(t))^{n}
= Fn of r × Fn of p and r × Fn of only p

Likelihood theory (and common sense) suggests that the only thing of any importance here is the part that depends on both r and p. In other words, if you want to maximize the likelihood only the middle bit counts.

If we go one step further, and consider the log-likelihood

log \binom{n}{r} + r × log (\frac{p(t)}{1 − p(t)}) + n × log (1 − p(t))

then we see that modelling might be easy and natural (and the asymptotic statistics might work) if we form models that are linear in

log (\frac{p(t)}{1 − p(t)})

That is exactly what we will do.

The function

log (\frac{p}{1 − p})

is called the logistic transformation. If

x = log (\frac{p}{1 − p})

then the inverse function is

p = \frac{exp(x)}{1 + exp(x)}

We will (obviously, from the name) model the logit of the probabilities as a function of the predictors.

Back to the Challenger Data

How do we want to model the challenger data? We (with the current data) only have temperature as a predictor, so we are stuck with models of the form,

logit(p(t)) = some function oft

The obvious first step is to plot the data (how can we do that?), then fit some simple models, then think about transformations. In this example we are particularly interested in extrapolations (dangerous though that is!) to temperatures that are less than we have seen.

Plotting the data

This is harder than one might imagine. One (mostly) can not just plot the logits, because much of the data is at the extremes where p = 0 or (1 − p) = 0. Often the best we can do is the plots that we did earlier, perhaps adding some sort of smooth line to the data. (Next example, we will try that). We will be able to do better once we have a model and can consider residuals—but be careful which residuals we think about.

The model

The S-Plus logistic regression command needs to know both the number of successes and the number of failures, so we need to start out by creating such a variable,

> FN <- cbind(Fail, n-Fail)

Now, just use the glm command, with the family=binomial option, to fit the model, viz.

> mymod <- glm( FN ~ Temp, family=binomial)
> mymod <- glm( cbind(Fail, n-Fail) ~ Temp, family=binomial)

Either of the above forms will work, but the first is easier in the long run. Now look at the results of the fit,

> summary(mymod)
Call: glm(formula = FN ~ Temp, family = binomial)
Deviance Residuals:
        Min         1Q     Median          3Q      Max
 -0.9522773 -0.7829968 -0.5411832 -0.04378972 2.651492
Coefficients:
                 Value Std. Error   t value
(Intercept)  5.0848620  3.0477855  1.668379
       Temp -0.1155992  0.0469317 -2.463136
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 24.23036 on 22 degrees of freedom
Residual Deviance: 18.08633 on 21 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
     (Intercept)
Temp -0.9932182

> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: FN
Terms added sequentially (first to last)
     Df Deviance Resid. Df Resid. Dev
NULL                    22   24.23036
Temp  1 6.144035        21   18.08633

The model fits well (although, just like R^{2} is regression, one should be cautious about this, and the temperature coefficient is certainly significant.

At this point we need to

• Get predictions for the low-temperature extrapolations, and

• Evaluate the fit of the model

Extrapolating Lets work on that, in that order (which is really back-to-front, but useful for teaching purposes).

The fitted values (on the probability scale) are obtained via,

> fitted(mymod)
         1         2         3          4          5          6          7
 0.2607849 0.1817866 0.1652202 0.09994062 0.07278457 0.06535816 0.06535816
          8          9         10         11         12         13         14
 0.06535816 0.05864156 0.05257638 0.04710711 0.04710711 0.04710711 0.04710711
         15         16         17         18         19         20        21
 0.03775043 0.03376847 0.02698616 0.02698616 0.02411117 0.02411117 0.0192299
         22         23
 0.01716665 0.01367159

The predicted values (on the logit scale) are,

> predict(mymod)
         1        2         3         4         5         6         7
 -1.041893 -1.50429 -1.619889 -2.197885 -2.544682 -2.660281 -2.660281
         8        9        10        11        12        13        14
 -2.660281 -2.77588 -2.891479 -3.007079 -3.007079 -3.007079 -3.007079
        15        16        17        18        19        20        21
 -3.238277 -3.353876 -3.585074 -3.585074 -3.700673 -3.700673 -3.931872
        22        23
 -4.047471 -4.278669
> fitted(mymod) - unlogit(predict(mymod))
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

(unlogit is the following simple function),

> unlogit
function(x)
exp(x)/(1 + exp(x))

For communication purposes, it makes sense to work on the probability scale (and hence use fitted), but for all other purposes the predicted values are what we want. In particular, the predict function has an option that will give the Standard Errors of the predicted values (on the logit scale), viz.,

> predict(mymod,se=T)
$fit:
         1        2         3         4         5         6         7
 -1.041893 -1.50429 -1.619889 -2.197885 -2.544682 -2.660281 -2.660281
         8        9        10        11        12        13        14
 -2.660281 -2.77588 -2.891479 -3.007079 -3.007079 -3.007079 -3.007079
        15        16        17        18        19        20        21
 -3.238277 -3.353876 -3.585074 -3.585074 -3.700673 -3.700673 -3.931872
        22        23
 -4.047471 -4.278669
$se.fit:
         1         2        3         4         5         6         7
 0.7682786 0.5943294 0.556385 0.4298909 0.4298809 0.4441498 0.4441498
         8         9        10        11        12        13        14
 0.4441498 0.4647341 0.4908397 0.5216385 0.5216385 0.5216385 0.5216385
        15        16        17        18        19        20        21
 0.5942932 0.6348851 0.7222049 0.7222049 0.7682357 0.7682357 0.8637766
        22       23
 0.9129226 1.013299
$residual.scale:
[1] 1.189907
$df:
[1] 21

So, a 95% confidence interval for the first predicted value on the logit scale would be

−1.041893 ± 2 × 0.7682786

Transferred to the probability scale, this is an interval (0.070, 0.621), with a center at 0.260. But what about predictions nearer the take-off temperature?

 predict(mymod, data.frame(Temp=c(30,31,32,33)), se=T)
$fit:
        1        2        3       4
 1.616887 1.501288 1.385689 1.27009
$se.fit:
        1        2        3       4
 1.972255 1.917739 1.863302 1.80895

(That was pretty sneaky, but the predict function requires a data frame predictor values).

The bottom line is that at 31^{∘}F our model extrapolates to p = 0.83, with a 95% confidence interval of (0.089, 0.996), which is VERY wide. Even with this wide interval, the point estimate is something to worry about.

Evaluating the Model

> plot(mymod)

Notice that the residual plot has at least two recognizable components, but it couldn’t help but look like this. There is perhaps some suggestion that a function of temperature might be needed, like perhaps t^{2}, but don’t be fooled. The places where there are big residuals are for temperatures where there were both no failures, and failures. I tried a number of additional functions of temperature, and nothing much worked.

The Zero-One version of the model

Very quickly, lets look at the model when we just look at a failure or not.

> Bn <- cbind(sign(Fail), 1 - sign(Fail))
> Bn
   [,1] [,2]
 1    1    0
 2    1    0
 3    1    0
 4    1    0
 5    0    1
 6    0    1
 7    0    1
 8    0    1
 9    0    1
10    0    1
11    0    1
12    0    1
13    1    0
14    1    0
15    0    1
16    0    1
17    0    1
18    1    0
19    0    1
20    0    1
21    0    1
22    0    1
23    0    1

then

> mymod <- glm(Bn ~ Temp, family=binomial)
> summary(mymod)
Call: glm(formula = Bn ~ Temp, family = binomial)
Deviance Residuals:
       Min       1Q     Median        3Q      Max
 -1.061112 -0.76128 -0.3782822 0.4523828 2.217506

Coefficients:
                  Value Std. Error   t value
(Intercept)  15.0422911  7.3366528  2.050293
       Temp  -0.2321537  0.1076141 -2.157279
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 28.26715 on 22 degrees of freedom
Residual Deviance: 20.31519 on 21 degrees of freedom
Number of Fisher Scoring Iterations: 4
Correlation of Coefficients:
     (Intercept)
Temp -0.99716
> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: Bn
Terms added sequentially (first to last)
     Df Deviance Resid. Df Resid. Dev
NULL                    22   28.26715
Temp  1  7.95196        21   20.31519

So far things look very similar, including the residual plot (believe me).

Another Example, Stroke

Background

The data are from a prospective study whose purpose was to characterize the spontaneous recovery of language functioning in patients during the first 3 to 4 months following a stroke. Details of the design and results of this study may be found in Holland et al. (1985) and Holland, Greenhouse, Fromm, and Swindell (1989). Briefly, patients admitted to a hospital with the diagnosis of stroke and who met the study inclusion criteria were seen daily (6 days/week) for 15 minutes beginning 48 to 72 hours poststroke and continuing until hospital discharge. The daily visits consisted of semistructured conversational interactions and systemic scoring of the patient’s communication behavior. At hospital discharge, patients were formally tested using the Western Aphasia Battery (WAB). If the WAB indicated aphasic dysfunction (AQ < 93.8), the patient was seen and assessed at 1 month postdischarge. If this second WAB demonstrated continued difficulty the patient was assessed again at 2 months postdischarge. The data set to be discussed here consists of 61 patients who had unilateral (right- or left-hemisphere) stroke, ischemic or hemorrhagic etiology, and impaired language functioning at time of hospital discharge. We will use these data to investigate the correlates of language recovery during the two month period post-hospital discharge. Definitions and the codes for the variables of interest are presented in Table 2.

The multiple logistic regression model will allow us to investigate simultaneously the effect of a number of different variables on the likelihood of spontaneous language recovery. A feature of the logistic regression model is that it allows us to use explanatory variables that are either quantitative or discrete. Using the variables defined and coded in Table 2 for the stroke recovery study, the specification of the deterministic component of the multiple linear logistic regression model on the logit scale is

ln(\frac{π}{1−π}) = β_{0} + β_{1} · (SIDE) + β_{2} · (GENDER) + β_{3} · (AGE) + β_{4} · (LOSD) + β_{5} · (W AB1). (1)

Note that if we were to dichotomize length of stay, age and WAB1 then a multiple crossclassification table of these data would consist of 2^{6} = 64 cells or subgroups and that with 61 patients such a table would have numerous sparse or empty cells which would be a problem for a chi-square analysis. An advantage of a statistical model such as the one given in equation (1), however, is that the model provides a structure that smooths the fit of the probability of a response to all the cells in such a table.

Side  Gender  Age  LOSD  WAB1  Status
  0     1      74   17   92.5    1
  1     1      75   30   87.8    0
  0     1      72   17   65.9    0
  0     1      60   21   86.4    1
  0     0      69   11   90.0    1
  0     0      62   13   6.2     0

Table 1: Example of what the data file ’stroke.dat’ looks like

Interpretation of Regression Coefficients

The output below presents part of the S-Plus output listing the estimates of the regressions coefficients and their standard errors. We see that the coefficients for gender, length-of-stay, and the WAB score at discharge are all significantly related to the probability of spontaneous language recovery at the 0.05 level. The coefficient for length-of-stay is negative indicating that the longer patients are in the hospital the smaller the chance of spontaneous language recovery relative to patients who have shorter lengths-of-stay even after controlling for the effect of the WAB1 score. The coefficient for WAB1 is positive as we would expect. Males, after controlling for the other variables in the model, are exp(2.88) = 17.8 times as likely to recover spontaneous language function as females. A 95% confidence interval for β_{2}, the

Variable Description  Values/Codes
STATUS Status at last  0 = Not Recovered
observation  1 = Recovered
SIDE  Side of lesion  0 = Left

1 = Right

GENDER  0 = Female

1 = Male

AGE  Age at stroke  In years

onset

LOSD  Length of hospital In days

stay

WAB1  WAB at  Range: 0 - 100

hospital discharge

Table 2: Variables and Codes from the Stroke Recovery Study

adjusted log odds ratio of recovery for males relative to females, however, is large and ranges from 0.47 to 5.29. We also note that side of stroke and subject’s age do not have a statistically significant association with spontaneous language recovery after controlling for the effects of the other variables. Details about these tests can be found in Hosmer and Lemeshow (1989, Chapter 5).

> stroke <- read.table("stroke.dat", head=T)
> stroke[1:10,]
  Side Gender Age LOSD WAB1 Status
1    0      0  79   27 18.2      0
2    0      0  82   17 86.0      0
3    1      0  72   42 46.3      0
4    0      0  59   23 35.6      0
5    0      0  60   28  0.0      0
6    1      1  75   12 81.9      1
7    0      0  75   53 56.0      0
8    0      0  85   46 17.3      0
9    0      0  79   58 74.9      0
10    0      1  82   34 18.2      0

> mymod <- glm(Status ~ Side + Gender + Age + LOSD +WAB1, family=binomial,
+ data=stroke)
> summary(mymod)
Call: glm(formula = Status ~ Side + Gender + Age + LOSD + WAB1, family = binomial,
        data = stroke)
Deviance Residuals:
       Min         1Q       Median            3Q      Max
 -1.317524 -0.2999725 -0.003778562 -0.0002212266 2.862154
Coefficients:
                    Value Std. Error    t value
(Intercept) -11.072759221 6.28633143 -1.7614024
       Side  -1.369026127 1.10318267 -1.2409786
     Gender   2.876569185 1.25325078  2.2952862
        Age  -0.006938071 0.03062231 -0.2265691
       LOSD  -0.196137069 0.09595028 -2.0441532
       WAB1   0.164467235 0.07172681  2.2929674
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 65.71913 on 60 degrees of freedom
Residual Deviance: 28.9878 on 55 degrees of freedom
Number of Fisher Scoring Iterations: 7
Correlation of Coefficients:
       (Intercept)       Side     Gender        Age       LOSD
  Side  0.2553975
Gender -0.5007016  -0.2793548
   Age -0.3373404  -0.1265137  0.2006563
  LOSD  0.1199385   0.2046374 -0.5309336 -0.1894954
  WAB1 -0.9292116  -0.3055970  0.5278826  0.0616042 -0.2969773
> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: Status
Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev

  NULL                    60   65.71913
  Side  1  0.10011        59   65.61902
Gender  1  3.36010        58   62.25892
   Age  1  1.45449        57   60.80444
  LOSD  1 12.25425        56   48.55018
  WAB1  1 19.56238        55   28.98780
>

Analysis of the Fit of the Model

Figure 1 presents a residual analysis for this logistic regression.

One relatively easy method for the analysis of the fit of the logistic regression model is based on calculating the predicted probability of spontaneous language recovery based on the model in (1) compared to the observed probability of language recovery Specifically, substitute the estimated regression coefficients given in above into the multiple logistic regression equation in (1). For each patient in the study calculate his or her individual predicted probability of language recovery by calculating the linear combination of each patient’s individual set of explanatory variables with the estimated regression coefficients and use the inverse of the logit transformation to solve for the probability of language recovery, π̂. For example, consider a patient whose stroke was on the left side (SIDE=0), was a 61 year old male (GENDER=1, AGE=61), stayed in the hospital 19 days (LOSD=19), and whose WAB at discharge was 81.6 (WAB1=81.6). The estimated probability of spontaneous language recovery based on the model for such a patient is

π̂ = \frac{1}{1 + e^{−(−11.1−1.37·(0)+2.88·(1)−0.0069·(61)−0.196·(19)+0.165·(81.6))}} = 0.77.

This value can be interpreted to mean that 77 patients out of 100 with the above characteristics can be expected to recover normal language functioning spontaneously within the first 2 months post-hospital discharge. To study the goodness of fit of the model, we obtain the frequency distribution of the values of π̂, and determine, say, the quintiles of the distribution. The observed number of patients who had normal language functioning are tallied within each quintile, and the expected number of patients who had normal language functioning (predicted by the model) is determined by summing the π̂’s for all the patients in the quintile. Table 3 presents the expected and observed number of patients who had normal language functioning for each quintile of risk of language recovery. Comparisons of the expected number predicted by the model to the observed number suggests that this model fits the data reasonably well.

A test for whether there is a significant difference between the observed and expected numbers in Table reftab:fit can be constructed based on a Pearson chi-square test. This chi-square test would have degrees-of-freedom equal to the number of groups (e.g., 5 for the number of quintiles) minus 2 (Hosmer and Lemeshow, 1989, chapter 5). In our example, this chi-square statistics is equal to 3.60 with 3 degrees-of-freedom and has a non-significant p-value=0.69, indicating that the model adequately fits the data. A small p-value would indicate a significant lack of fit of the model.

Expected Observed
Quintile of risk Number Number
Highest  5  10.7  11
4  3.1  2
3  0.2  1
2  0.0  0
Lowest  1  0.0  0
Total  14.0  14

Table 3: Analysis of the fit of the logit model in Equation (1) Figure 1: Residual plots

More Logistic Regression

How to plot these data

Look at the stroke data. This is explicitly 0–1 data, and it is very hard to think about plotting it. A simple plot of, say, Status against WAB1, looks like the following, and it not all that informative.

> plot(WAB1, Status)

Even if I jitter the points a little, it is still not very easy to see any trends or patterns.

> plot(WAB1, Status+0.1*runif(length(Status)))

A somewhat more useful approach is to plot the data and add some sort of smooth curve that will allow you to at least look at some of the trends in the data. The particular approach I took with these data was to superimpose a (black-box) smooth curve generated by loess. The various pictures result from different smoothing parameters.

> scatter.smooth(WAB1, Status,span=0.2)
> title("Span = 0.2")
> scatter.smooth(WAB1, Status,span=0.35)

> title("Span = 0.35")
> scatter.smooth(WAB1, Status,span=0.5)
> title("Span = 0.5")
> scatter.smooth(WAB1, Status,span=0.65)
> title("Span = 0.65")

This is still not perfect, but you certainly get the sense that the probability goes up, probably more or less monotonely with WAB1.

Compare this to the similar picture(s) for Age.

> for (i in c(0.2, 0.35, .5, 0.65)){
+ scatter.smooth(Age, Status, span=i)
+ title(paste("Span = ", i))
+ }

If anything, there is a bit of evidence for some type of quadratic effect of Age, but the probable story is that Age is not that important, on its own. We should always remember that variables that do not seem to be important (from these marginal plots), may be quite important when paired with other variables (that is viewed conditionally).

With that hint about the Age variable, I re-ran the logistic regression model for these data and obtained,

> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age + Age^2, family=bin>
> summary(mymod)
Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2, family =
        binomial)
Deviance Residuals:
       Min         1Q        Median            3Q      Max
 -1.317304 -0.2637303 -0.0002506098 -6.014239e-07 2.770847
Coefficients:
                   Value   Std. Error   t value
(Intercept) -42.13652444 17.240544854 -2.444037
       Side  -1.70612965  1.320176507 -1.292350
     Gender   4.03384776  1.847166535  2.183803

       LOSD  -0.28388626  0.130491084 -2.175522
       WAB1   0.24891576  0.104667163  2.378165
        Age   0.92453559  0.421765743  2.192059
   I(Age^2)  -0.00791642  0.003578441 -2.212254
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 65.71913 on 60 degrees of freedom
Residual Deviance: 21.7621 on 54 degrees of freedom
> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: Status
Terms added sequentially (first to last)
         Df Deviance Resid. Df Resid. Dev
    NULL                    60   65.71913
    Side  1  0.10011        59   65.61902
  Gender  1  3.36010        58   62.25892
    LOSD  1 13.46451        57   48.79441
    WAB1  1 19.75552        56   29.03889
     Age  1  0.05110        55   28.98780
I(Age^2)  1  7.22570        54   21.76210

This is somewhat surprising (because I didn’t think there was an Age effect). Could we have seen the same trend from a residual analysis of these data?

Let us go back to the original model and look at residuals.

There is at least one large (Pearson) residual, and it is the same point as the large (deviance) residual. What about a plot of (deviance) residual against Age (when Age is in the model).

> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age, family=binomial)
> scatter.smooth(Age, (residuals(mymod)))
> identify(Age, (residuals(mymod)))
[1] 59
> scatter.smooth(Age[-59], (residuals(mymod))[-59])

We see that observation 59 is the one with the large residual. If I plot the data (the lower plot) without that observation, the quadratic trend seems to go away, at least a bit.

Observation 59 is now suspect, so let us remove that observation, and see what the resulting model looks like.

> WW <- rep(1, length(Side))
> WW[59] <- 0
> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age,
+ weight=WW, family=binomial)
Warning messages:
  linear convergence not obtained in  10  iterations. in: glm.fitter(x = X, y =
        Y, w = w, start = start, offset = offset, family = family, ....

Oh-oh, what does that mean? (Basically, the optimization algorithm is converging slowly, which is an indication that there are some fitted probabilities that are very close to 0 or 1). We can try a few more iterations, via.

> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age,
+ weight=WW, maxit=30,family=binomial)

Now, what does the summary and anova look like?

> summary(mymod)

Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age, family = binomial,
        weights = WW, maxit = 30)
Deviance Residuals:
      Min           1Q        Median            3Q      Max
 -1.21685 -0.006138687 -2.107342e-08 -2.107342e-08 2.120409
Coefficients:
                    Value  Std. Error     t value
(Intercept) -46.034665461 23.60416082 -1.95027757
       Side  -7.479606389  4.79366143 -1.56031178
     Gender   7.889726328  4.28891089  1.83956406
       LOSD  -0.677382165  0.43418335 -1.56012928
       WAB1   0.662724777  0.34770876  1.90597667
        Age  -0.001289562  0.04096666 -0.03147832
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 62.71879 on 60 degrees of freedom
Residual Deviance: 13.32727 on 54 degrees of freedom
Number of Fisher Scoring Iterations: 10
Correlation of Coefficients:
       (Intercept)       Side     Gender       LOSD       WAB1
  Side  0.8469550
Gender -0.9187997  -0.8860396
  LOSD  0.7944690   0.9174434 -0.8555061
  WAB1 -0.9831416  -0.9061715  0.9317554 -0.8758313
   Age -0.0387868   0.0615410  0.0243350 -0.0022418 -0.0596445
Warning messages:
  1 rows with zero weights not counted in: summary(mymod)

If I now start adding the quadratic term in Age, I get myself into all sorts of trouble.

> formula(mymod)
Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2
> anova(mymod)
Analysis of Deviance Table
Binomial model
Response: Status

Terms added sequentially (first to last)
         Df Deviance Resid. Df Resid. Dev
    NULL                    60   62.71879
    Side  2  0.62858        58   62.09021
  Gender  1  2.34693        57   59.74328
    LOSD  1 13.67051        56   46.07277
    WAB1  1 32.74451        55   13.32826
     Age  1  0.00099        54   13.32727
I(Age^2)  1 13.32727        53    0.00000
Warning messages:
1: linear convergence not obtained in  10  iterations. in: glm.fit(x, y, w,
        offset = offset, family = family, maxit = control$maxit,  ....
2: linear convergence not obtained in  10  iterations. in: glm.fit(x, y, w,
        offset = offset, family = family, maxit = control$maxit,  ....

Whoa. What is going on here? We suddenly have a model that fits perfectly?

The summary looks even weirder?

> summary(mymod)
Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2, family =
        binomial, weights = WW, maxit = 50)
Deviance Residuals:
         Min            1Q        Median            3Q          Max
 -3.9065e-05 -2.107342e-08 -2.107342e-08 -2.107342e-08 3.281089e-05
Coefficients:
                    Value  Std. Error       t value
(Intercept) -2858.6085408 2227320.001 -0.0012834297
       Side  -164.0618144  165491.082 -0.0009913635
     Gender   273.5752941  237195.052  0.0011533769
       LOSD   -16.2464405   13633.473 -0.0011916582
       WAB1    24.1122280   17153.825  0.0014056473
        Age    39.8782980   35880.527  0.0011114190
   I(Age^2)    -0.3549147     302.566 -0.0011730157
(Dispersion Parameter for Binomial family taken to be 1 )
    Null Deviance: 62.71879 on 60 degrees of freedom
Residual Deviance: 0 on 53 degrees of freedom
Number of Fisher Scoring Iterations: 29
Correlation of Coefficients:

         (Intercept)       Side     Gender       LOSD       WAB1        Age
    Side  0.3661039
  Gender -0.8397332  -0.5132934
    LOSD  0.7676113   0.5867586 -0.7405370
    WAB1 -0.9678428  -0.5023038  0.8384746 -0.7822973
     Age -0.9537356  -0.2802197  0.7986261 -0.8085521  0.8643352
I(Age^2)  0.9585630   0.3372172 -0.8129406  0.8428222 -0.8863173 -0.9963403
Warning messages:
1: fitted values close to 0 or 1 in: family$deviance(mu, y, w, residuals = T)
2: 1 rows with zero weights not counted in: summary(mymod)
>

On reflection, it is clear what is going on. Look at the first plot of Status against WAB1. Here we see, that absent the one observation in the middle (which is 59), we have an almost perfect discimination between the two levels of Status. That is the clue to what is “broken”. With the 7 parameters (anything like WAB1 squared would have done), we can find a function that perfectly discriminates between the two levels of status, and then maximum likelihood modelling breaks down. There is a lot more to be said here, I just don’t know what it is.

End of Stroke example

Influence of a point

There are good measures of influence for glm models, but they don’t seem to be implemented (i.e., I could not find them) in the standard packages. Rather than just give up, I quickly hacked together a “Poor Man’s” influence function. The basic idea of the following code is to fit the glm model with each data point in turn removed. The summary is then a plot of each of the coefficients, with as many points on the plot as there are data points. This plot is surprisingly effective and quickly shows us where the pathology in the data lies.

# Heavily annotated, for class
# Generate a vector of all ones, of length the number of observations.
W <-  rep(1,length(Side))
# Generate current weights.
WW <- W
# Fit the base model
mymod <- glm(Status ~ Gender + LOSD + Side + Age + WAB1, weight=WW,
   family=binomial)
# Save the coefficients, just in case
base.coef <- coefficients(mymod)
# Set up a matrix to hold the coefficients for each model with one
# point deleted.
infs <- matrix(0, nrow=length(Side), ncol=length(coefficients(mymod)))
# Loop through the data, setting one weight equal to 0.
# Save the coefficients in the infs matrix
for (i in 1:length(Side)){
        WW<-W; WW[i] <- 0;
        infs[i,] <- coefficients(update(mymod, . ~ .))
}
## Can't think of an easy, smart summary, so just plot the
## coefficients.
#Set up a 2 x 3 plotting frame
par(mfrow=c(2,3))
#Do the 6 plots.  After each plot, allow the user to
# identify the points which look interesting.
for(i in 1:6){
        plot(1:length(Side), infs[,i])
        title(names(coefficients(mymod))[i])

        identify(1:length(Side), infs[,i])
}

After running that code, and interacting with the various plots, you end up with the following plot.

There are other things we could have plotted (like following the deviance as the models changed), but this plot alone tells us much of the story.

More Examples

Here I will pull together a number of examples, drawn directly from the Book by Venables and Ripley.

First, the budworm data. The following data is the result of an experiment on the toxicity of the tobacco budworm to doses of particular drug, that the moths are beginning to show resistance to. Batches of 20 moths of each sex were exposed for 3 days to the chemical and the number of dead (or disabled) moths in each group were recorded. The data is,

Dose Response Gender
1       1       M
2       4       M
4       9       M
8       13      M
16      18      M
32      20      M
1       0       F
2       2       F
4       6       F
8       10      F
16      12      F
32      16      F

The doses were in micro-grams, and the dose level when up by powers of two. It is a judgement call as to whether you will try to model with dose as the response or log_{2}(dose. We will try both.

> bud <- read.table("budworm.dat", head=T)
> bud
   Dose Response Gender
 1    1        1      M
 2    2        4      M
 3    4        9      M
 4    8       13      M
 5   16       18      M
 6   32       20      M
 7    1        0      F
 8    2        2      F
 9    4        6      F
10    8       10      F
11   16       12      F
12   32       16      F
> attach(bud,pos=1)
> N <- 20

First, lets plot the data, and see if anything obvious shows up.

> plot( Dose, Response/N)
> points( Dose[Gender=="F"], (Response/N)[Gender=="F"], pch="F")
> points( Dose[Gender=="M"], (Response/N)[Gender=="M"], pch="M")

The same plot on the log(Dose) scale, is

> plot( Dose, Response/N, log="x")
> points( Dose[Gender=="F"], (Response/N)[Gender=="F"], pch="F")
> points( Dose[Gender=="M"], (Response/N)[Gender=="M"], pch="M")

Given that we have more than 0-1 data, it makes a lot of sense to at least try to plot the (observed) logits,

> par(mfrow=c(2,1))
> plot( Dose, logit(Response/N), log="x")
> title("Logit vs Log(Dose)")
> plot( Dose, logit(Response/N))
> title("Logit vs Dose")

This plot makes it quite clear, that the log scale has a much better chance of working. So, proceed to fit the model. We surely want a model with different slopes for the two

sexes

> ldose <- log(Dose)
> SF <- cbind(Response, N - Response)
> mymod <- glm(SF ~ Gender * ldose, family=binomial)
> summary(mymod)
 ..........
Coefficients:
                  Value Std. Error    t value
 (Intercept) -2.9935414  0.5525295 -5.4178852
      Gender  0.1749865  0.7781556  0.2248733
       ldose  1.3071341  0.2410133  5.4234939
Gender:ldose  0.5091459  0.3894474  1.3073547

    Null Deviance: 124.8756 on 11 degrees of freedom
Residual Deviance: 4.993727 on 8 degrees of freedom
> anova(mymod,test="Chisq")
Analysis of Deviance Table

Binomial model
Response: SF
Terms added sequentially (first to last)
             Df Deviance Resid. Df Resid. Dev   Pr(Chi)
        NULL                    11   124.8756
      Gender  1   6.0770        10   118.7986 0.0136955
       ldose  1 112.0415         9     6.7571 0.0000000
Gender:ldose  1   1.7633         8     4.9937 0.1842088
>

Now, what is going on here? Gender seems to have no effect (based on the t-value), but the deviance is quite large? Even more, our plots seem to imply what looks like a gender effect. The best way to find out is to distrust the parameter stuff (it can get easily confused), and try fitting a model without the interaction. This makes sense anyway, because there is little/no evidence of an interaction.

> anova(mymod <- glm(SF ~ Gender + ldose, family=binomial))
Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev
  NULL                    11   124.8756
Gender  1   6.0770        10   118.7986
 ldose  1 112.0415         9     6.7571
> summary(mymod)
Coefficients:
                Value Std. Error   t value
(Intercept) -3.473154  0.4682939 -7.416612
     Gender  1.100743  0.3557226  3.094385
      ldose  1.535336  0.1890119  8.122959
    Null Deviance: 124.8756 on 11 degrees of freedom
Residual Deviance: 6.757064 on 9 degrees of freedom

That is a little more like it, and seems to be a VERY good fitting model. We should at least look at some diagnostics, and some other plots, but unless they show something very odd, we are done.

> plot(mymod)

That looks very good, we should finish by plotting the data and fitted curves on the same graph, and that is left as an exercise for you (or look it up in the book). The book also goes on to talk about LD50, and how to estimate it

Time Series

I will be drawing my motivation and some of my lecture material from the book Time Series: A Biostatistical Introduction by Peter J. Diggle. I’ll also draw some material from the book The Analysis of Time Series: An Introduction by C. Chatfield, and various other sources as needed.

Many time series books draw their motivation from economic and finance, and much of the methodology is widely applied there. However the general ideas are useful in all areas of statistics from Biostatistics (the approach I will follow), to all sorts of physical phenomena (like ocean temperatures at different depths), to geology (where I first met this material— analyzing the core samples from exploratory drilling for nickel), to all sort of spatial problems (which need very similar mathematics).

The basic mathematical tools are some simple concepts about correlation, some ideas from other stochastic processes (things like Markov chains), and some Fourier theory. We will have to develop some mathematics, but we will try to motivate the work with examples.

What is a time series

In its simplest form a time-series is a set of data y_{t} : t = 1, …, n in which the subscript t indicates the time at which the datum is observed. That notation seems to imply that the data points must be equally spaced in time. While that is certainly the most common case, it is not necessary. A more honest notation might then be y(t_{i}) : i = 1, …, n. (Note that most software has difficulties with non-equally spaced data!)

Those two notations cover the most common situations, but there are a lot of bells and whistles to consider. The data may have a corresponding predictor series. For example we might monitor the progress, over time, of a fever in a baby. The predictor series might be the amount of medicine the child is given at varying times. Another common extension is to consider multivariate time-series where we measure multiple responses at the same time. To continue in the same vein, we might measure the temperature, pulse rate, and blood pressure of the baby.

What makes the data interesting (and difficult to analyze) is that it is NOT reasonable to assume that the observations are iid. At the very least we need to be willing to drop the independent part of the assumptions—while we might want to retain the identically distributed. In the above “examples” it makes little sense to assume that the baby’s temperature at one time point is independent of the temperature at the next time. If YOU were the parent, the whole point of keeping records would be to discover when the temperature is rising/falling, and you probably expect the observations to be correlated; which implies they are not independent.

Example 1

Levels of lutenizing hormone in blood samples. Blood samples are taken from a healthy woman every ten minutes over an eight hour period. The samples are assayed to determine the lutenizing hormone (LH) levels. The following data and plots show three series—one from the early follicular phase of the subject’s menstrual cycle, the others from the late follicular phase of two successive cycles. There are clear cycles. In principle, the phenomenon being investigated is happening in real time, but we can’t collect and assay the data that way. It makes a lot more sense to do equally spaced intervals. Here is the data.

   Lutenizing hormone.
   col 1 : serial number
   col 2 : Cycle 1, early follicular phase 1 (not discussed)
   col 3 : Cycle 1, late follicular phase
   col 4 : Cycle 2, early follicular phase
   col 5 : Cycle 2, late follicular phase
Serial Early1 Late1 Early2 Late2
 1  2.2  5.5  2.4  4.3
 2  2.2  4.5  2.4  4.6
 3  2.3  5.1  2.4  4.7
 4  2.0  5.5  2.2  4.1
 5  1.6  5.7  2.1  4.1
 6  1.4  5.1  1.5  5.2
 7  1.8  4.3  2.3  5.0
 8  2.2  4.8  2.3  4.4
 9  2.9  5.6  2.5  4.2
10  2.6  5.9  2.0  5.1
11  2.4  6.0  1.9  5.1
12  2.1  5.1  1.7  4.7
13  3.0  5.2  2.2  4.4
14  2.5  4.4  1.8  3.9
15  2.7  5.5  3.2  5.4
16  2.2  5.4  3.2  5.9
17  2.4  4.1  2.7  4.2
18  2.7  4.4  2.2  4.1
19  3.1  4.7  2.2  4.1
20  2.5  4.6  1.9  3.6
21  2.4  6.0  1.9  3.1
22  2.3  5.6  1.8  4.8
23  2.4  5.1  2.7  5.1
24  1.9  4.7  3.0  5.1
25  3.3  4.8  2.3  4.5
26  3.8  5.5  2.0  4.6
27  3.7  5.1  2.0  5.8
28  3.5  5.2  2.9  5.0
29  3.1  5.0  2.9  5.1
30  2.7  4.0  2.7  4.5
31  4.1  3.7  2.7  4.2

32  4.0  4.8  2.3  6.0
33  3.4  5.9  2.6  5.6
34  3.2  5.5  2.4  5.4
35  3.7  4.9  1.8  5.0
36  3.6  4.4  1.7  4.4
37  4.1  4.7  1.5  4.6
38  2.0  4.2  1.4  5.7
39  4.6  5.5  2.1  5.2
40  4.1  4.9  3.3  5.0
41  3.2  4.8  3.5  4.4
42  2.9  4.5  3.5  5.7
43  2.7  4.9  3.1  5.7
44  3.0  4.9  2.6  4.8
45  0.   4.5  2.1  3.4
46  0.   4.2  3.4  5.5
47  0.   4.9  3.0  5.5
48  0.   5.9  2.9  5.6

Here are some simple analyses in S.

> lh <- read.table("lh.dat", head=T)
> lh
 Serial Early1 Late1 Early2 Late2
1      1    2.2   5.5    2.4   4.3
2      2    2.2   4.5    2.4   4.6
3      3    2.3   5.1    2.4   4.7
. . . . .
> early <- ts(lh$Early2)
> late1 <- ts(lh$Late1)
> late2 <- ts(lh$Late2)
summary(lh)
    Serial          Early1          Late1           Early2
Min.   : 1.00   Min.   :0.000   Min.   :3.700   Min.   :1.40
1st Qu.:12.75   1st Qu.:2.200   1st Qu.:4.575   1st Qu.:2.00
Median :24.50   Median :2.650   Median :4.900   Median :2.30
Mean   :24.50   Mean   :2.594   Mean   :4.990   Mean   :2.40
3rd Qu.:36.25   3rd Qu.:3.225   3rd Qu.:5.500   3rd Qu.:2.75
Max.   :48.00   Max.   :4.600   Max.   :6.000   Max.   :3.50
    Late2
Min.   :3.100
1st Qu.:4.375
Median :4.800
Mean   :4.800
3rd Qu.:5.250
Max.   :6.000

I want to plot the data on the same scale, so I need the minimum (1.4) and maximum (6.0).

> par(mfrow=c(2,2))
> tsplot(early,ylim=c(1.4,6.0))
> tsplot(late1,ylim=c(1.4,6.0))
> tsplot(late2,ylim=c(1.4,6.0))

Notice the periodic nature of the plots. Is this real periodicity, or are we kidding ourselves. Could we predict into the future? The data shows a cyclic nature, but the frequency of the cycle is probably not known in advance. Other data shows similar trends that are known in advance. These are often called seasonal trends.

As a digression, when you have multiple plots per page how do you get to the next page? I wrote the following function, which seems to do the trick. There might be a better way—I don’t know it.

> newplot <- function() {
+ par(mfrow=par("mfrow"))}

Example 2

UK deaths from bronchitis, emphysema, and asthma. The table and figure show the monthly returns of deaths in the UK from bronchitis, emphysema, and asthma over the years 1974– 1979, inclusive. Males and females are shown separately.

Males

2134 1863 1877 1877 1492 1249 1280 1131 1209 1492
1621 1846 2103 2137 2153 1833 1403 1288 1186 1133
1053 1347 1545 2066 2020 2750 2283 1479 1189 1160
1113 970 999 1208 1467 2059 2240 1634 1722 1801
1246 1162 1087 1013 959 1179 1229 1655 2019 2284
1942 1423 1340 1187 1098 1004 970 1140 1110 1812
2263 1820 1846 1531 1215 1075 1056 975 940 1081
1294 1341

Females

901 689 827 677 522 406 441 393 387 582
578 666 830 752 785 664 467 438 421 412
343 440 531 771 767 1141 896 532 447 420
376 330 357 445 546 764 862 660 663 643
502 392 411 348 387 385 411 638 796 853
737 546 530 446 431 362 387 430 425 679
821 785 727 612 478 429 405 379 393 411
487 574

We read in the data and do a plot

> uk.female <- scan("uk-female.dat")
> uk.female <- matrix(uk.female,ncol=2,byrow=T)[,2]
> uk.female <- ts(uk.female, start=74.0,freq=12)
> uk.female
     Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
74:  901  689  827  677  522  406  441  393  387  582  578  666
75:  830  752  785  664  467  438  421  412  343  440  531  771
76:  767 1141  896  532  447  420  376  330  357  445  546  764
77:  862  660  663  643  502  392  411  348  387  385  411  638
78:  796  853  737  546  530  446  431  362  387  430  425  679
79:  821  785  727  612  478  429  405  379  393  411  487  574

(Similarly for the males).

Now notice that there is a strong seasonal nature to the plots; the winter of 75–76 was particularly severe; and the females have fewer deaths than the males.

At one level there is an obvious relationship between the mal and female returns, namely their common seasonal pattern. Is there more than that (i.e., once the seasonality is removed, is there still a relationship?)

We will look at other examples as we need them.

What are we trying to do?

One aim is to provide a concise description of a historic series. The description may be just a short summary, or some graphs, or a quite complicated model for the underlying stochastic process.

If we are more ambitious, we might try to forecast future values of the series. This is particularly important for economic series.

An important aim in many medical examples is detecting a change in a series. This is also often a consideration in quality control situations.

Some Notation

Random variables will be upper case letters, and realized values will be lower case, thus the random variable might be Y_{t} : t = 1, …, n and the observed values would be y_{t} : t = 1, …, n.

The underlying time scale might be continuous or discrete. (For example, my blood pressure could be monitored continuously, but the CLOSING value of the Dow Jones index is discrete). For most purposed the granularity of the time scale does not matter. However one must also worry about the whether the series itself is discrete-valued or continuousvalued. We will not (seriously) consider discrete-valued time series. There are entire books on Markov chains, point processes and the like.

Trend and Stationarity

The concepts of trend, serial dependence, and stationarity are central to our work in timeseries data.

Let Y (t) be a random function of time. That is at each time point t, Y (t) has some distribution.

The (generally non-random function), μ(t) = E[Y(t)] is called the mean-values function, or trend of the series.

Secondly, the autocovariance function

γ(s, t) = E[Y(t) − μ(t)Y(s) − μ(s)]

describes how the series at time t correlates with the series at time s. Clearly γ(t, t) is the variance at time t.

The series is considered to be stationary if the probabilistic structure of Y (t) is unaffected by time. Loosely speaking, the series looks about the same everywhere. A formal definition of strict stationarity is that for any positive integer m and times t_{i}, i = 1, …, m, the joint distribution of Y (t_{i} + s) : i = 1, …, m does not depend on s.

Strict stationarity is a very strong assumption that is hard to check and rarely needed. Much more common (and what we shall use) is the idea of second-order stationarity, namely: μ(t) must be a constant and γ(t, s) must depend only on |t − s|, not the particular values of t and s.

The simplest stationary random sequence is white noise This is a sequence of mutually independent random variables, each with mean zero and finite variance σ^{2}. For such a sequence the autocovariance function is

γ(t, s) = \{ \begin{matrix}σ^{2} & ift = s \\ 0 & ift ≠ s \end{matrix}

There are many versions of white noise, one standard one is where each of the rv’s a gaussian distribtion. Here are a few realizations of such a series.

> tsplot(ts(rnorm(200)))
> tsplot(ts(rnorm(200)))
> tsplot(ts(rnorm(200)))
> tsplot(ts(rnorm(200)))
> plot( 1:300/30, dgamma(1:300/30, 1.5))
> tsplot(ts(rgamma(200,1.5)))
> tsplot(ts(rgamma(200,1.5)))
> tsplot(ts(rgamma(200,1.5)))

224

Time Series

Recall that we had just looked at some examples of simple white noise series with several different distributions. Of what possible use is a white-noise sequence? Well, at the least it is a base sequence that we can use to compare with other series—series that we believe have significant signals.

A nonstationary example

While we will be assuming stationarity (a lot), it is worthwhile looking a simple time-series that is not stationary, but when generated looks like it has a strong trend. Consider the process

Y_{t} = Y_{t−1} + Z_{t} : t = 1, 2, …

where {Z_{t}} is just white noise, and Y_{0} is 0. What might this sequence look like (we’ll generate some examples, so hold on). What are the properties of this random sequence?

For starters, μ(t) = E[Y_{t}] can be easily derived as

μ_{t} = μ_{t−1} + E[Z_{t}] = μ_{t−1}

which is quickly seen to be zero. But a calculation that is just as simple shows that γt, t = tσ_{2}, which violates the stationarity assumption. (Digression, here is the calculation:

E((Y_{t})(Y_{t})) = E((Y_{t} + Z_{t})(Y_{t} + Z_{t}))
= E(2Y_{t}Z_{t} + (Y_{t})^{2} + (Z_{t})^{2})
= 0 + E((Y_{t})^{2}) + σ^{2}

So what about some pictures? Note all the different forms that a random walk (what this process is called) can take!

Simple Descriptive tools

Plots

We have already seen ample plots.

Smoothing

Often the signal in the data is hard to see because of all the underlying variability. It is quite useful to be able to see a (slightly) smoothed version of the data, to help us see the signal.

A very simple smoother is a moving average. This is defined as a new series

s_{t} = (y_{t−1} + y_{t} + y_{t+1})/3

(obviously we can average over a larger number of observations, and get “more” smoothing). Lets look at an example (based on the uk.female data) that we saw before.

> par(mfrow=c(4,1))
## Create a simple smoothing function
> mean3 <- function(y) { (lag(y,-1) + y + lag(y,1))/3}

> tsplot( uk.female, mean3(uk.female), type="pl")
> title("UK females, and smooth 3")
> tsplot( uk.female - mean3(uk.female))
> title("Residuals")
> mean5 <- function(y) { (lag(y,-2) + lag(y,-1) + lag(y) + lag(y,1)
+ + lag(y,2) )/5}
> tsplot( uk.female, mean5(uk.female), type="pl")
> title("UK females, smooth 5")
> tsplot( uk.female - mean5(uk.female))
> title("Residuals")

Notice that the longer the moving average, the more we have smoothed the data? What are the residuals telling us? Have we captured the bulk of the data? The residuals allow us to concentrate on the random fluctuations, partially ignoring the underlying signal.

The simple 3 and 5 point averages are just examples of what is called a moving average of a time series, defined as

s_{t} = ∑_{j=−p}^{p} w_{j}y_{t+j} t = p + 1, …, n − p

where p is a positive integer and the w_{j} are weights which sum to 1. Typically each w_{j} is positive. The order of the moving average is 2p + 1. There are many different choices for the order of smoothing and the exact weights.

A second type of smoothing is spline smoothing, where the series is approximated by a spline (piecwise polynomial) function. The details are left to the reader, but here are some examples.

> plot(1:length(uk.female), uk.female)
> jj <- smooth.spline( 1:length(uk.female), uk.female,spar=.0000001)
> lines(predict(jj))
> plot(1:length(uk.female), uk.female - predict(jj)$y)
> title("Residuals")
> plot(1:length(uk.female), uk.female)
> jj <- smooth.spline( 1:length(uk.female), uk.female,spar=.0001)
> lines(predict(jj))
> jj <- smooth.spline( 1:length(uk.female), uk.female,spar=.00001)
> lines(predict(jj))
> plot(1:length(uk.female), uk.female - predict(jj)$y)
> title("Residuals")

Both of these smoothing techniques are an attempt to decompose the series into (smooth) signal plus (random) noise.

Differencing

Differencing provides a simple approach to removing, rather than highlighting, trends in time-series data.

The first difference of a time series {y_{t}}, written as {Dy_{t}} is defined by the transformation

Dy_{t} = y_{t} − y_{t−1}

Higher order differences are defined by repeated use of the above definition. Thus the second order difference {D^{2}y_{t}} is

D^{2}y_{t} = D(Dy_{t}) = Dy_{t} − Dy_{t−1} = y_{t} − 2y_{t−1} + y_{t−2}

Why would one ever difference? Suppose that the series is deterministic, and a linear function of time, say y_{t} = α + βt, then

Dy_{t} = (α + βt) − {α + β(t − 1)} = β

a constant for all t. Similarly the k’th order difference will completely remove a k’th order polynomial effect.

Differencing also has the special property that it will remove the sort of nonstationarity that we generated earlier

Y_{t} = Y_{t−1} + Z_{t} : t = 1, 2, …

If we consider differences of this series, we end up with just white noise (Z_{t}) — generally a desirable thing.

If differencing is used to remove a trend in a time series, the usual approach is to construct {Dy_{t}}, {D^{2}y_{t}}, etc until visual inspection indicated the differenced series is stationary. Generally first or second order differences are enough. A simple calculation shows that a first order difference is the same as subtracting a moving average from the series. As an example, lets generate some non-stationary data (as above), and difference it, to see the effect.

> dd <- ts(cumsum(rnorm(200)))
> tsplot(dd)
> title("Nonstationary series")
> tsplot(diff(dd))
> title ("DIfferenced series")
> tsplot(dd<- ts(cumsum(rnorm(200))))
> title("Nonstationary series")
> tsplot(diff(dd))
> title ("DIfferenced series")

Autocovariance and autocorrelation function

The autocovariance function of a stationary random function {Y (t)} is

γ(k) = cov{Y(t), Y(t − k)}

This does not not depend on t because of the the assumption of stationarity. Since γ(0) is the variance, the autocorrelation function is

ρ(k) = γ(k)/γ(0)

An example Consider the first-order autoregressive process defined by

Y_{t} = αY_{t−1} + Z_{t}

where Z_{t} is white noise and −1 < α < 1. Some calculations (you do them!) show that

ρ_{k} = αρ_{k−1}

Inserting ρ_{0} = 1 gives ρ_{k} = α^{k} Lets look at some (generated) data, and correlations/covariances.

231

Autocorrelation Function and Periodogram

Last time we looked at the autocorrelation function (acf) which just gave the correletion of a series with itself at various lags. This time we will also look at the periodogram, an estimate of the spectrum of the time series.

The periodogram is a summary description based on a representation of an observed time series as a superposition of sinusoidal waves of various frequencies. The practical use of this is the empirical observation that many times series show cyclic fluctuations, but the underlying sin and cosine fucntions are not obvious.

To motivate the periodogram, consider a very simple model for a time series y_{t} exhibiting cyclic fluctuations with a known period, p. We assume that the the series is

y_{t} = αcos(ωt) + βsin(ωt) + z_{t}

where z_{t} is a white noise sequence, ω = 2π/p is the known frequency of the cyclic fluctuations and α, β are parameters that we shall estimate (by least squares). (Note the relationship between period p and frequency 2π/p).

After some calculations (…, I’ll ignore the math), one finds that the least squares estimates are

\begin{aligned}α̂ & = & 2 \{∑_{t=1}^{n} cos(tω)\} /n \\ α̂ & = & 2 \{∑_{t=1}^{n} sin(tω)\} /n \end{aligned}

Second (random) fact. The regression sum of squares (sum of square of the fitted values) is

n(α̂^{2} + β̂^{2})/2

This is proportional to the squared amplitude α̂^{2} + β̂^{2} of the fitted sinusoid.

That is a quite simplistic model. More realistic would be if we not only did not know the amplitude but also didn’t know the frequency of the series.

So consider a model of the form

y_{t} = ∑_{k=1}^{m} {α_{k}cos(ω_{k}t) + β_{k}sin(ω_{k}t)} + z_{t}

This is again an easy model to estimate, provided we pick the right set of frequencies. If we take ω_{k} = 2πk/n, with k < n/2 then some nice things happen. For starters this is the only set of frequencies that we really can estimate (smaller periods would jiggle more often than the data, longer periods would not see an entire cycle in the data!). This set of frequencies also represents a set of orthogonal functions. That is the correlation between any pair of these funcions (measured at the data points) is zero.

The nice thing that happens (because of the orthogonality) is that the estimates of α̂ and β̂ are the same as before, namely,

\begin{aligned}α̂_{k} & = & 2 \{∑_{t=1}^{n} cos(tω_{k})\} /n \\ α̂_{k} & = & 2 \{∑_{t=1}^{n} sin(tω_{k})\} /n \end{aligned}

Furthermore, α̂_{k}^{2} +β̂_{k}^{2} again gives the sums of squares associated with each of the m fitted sinusoids. Finally, (via some math, or common sense), if all the frequencies are fit, then the regression model fits exactly.

Some further mathematical manipulations show that

α̂^{2} + β̂^{2} = 2 [{∑ y_{t}cos(ωt)}^{2} + {∑ y_{t}sin(ωt)}^{2}] n

This is quite intuitive, the terms are just the inner product of the data with the (basis) function. This is what regression (and particularly ANOVA) was all about.

We can now define the ordinates of the periodogram function, I(ω) for 0 <= ω <= π, by

I(ω) = [{∑ y_{t}cos(ωt)}^{2} + {∑ y_{t}sin(ωt)}^{2}] n

.

Thus we can “parition” the total sum of squares for the data y_{t} as

∑_{t=1}^{n} = I(O) + 2 ∑_{j=1}^{m} I(2πj/n) + I(π)

where m is the largest integer less than n/2 and the final term (I(ω)) is absent if n is odd. The function represents that power (amount of energy) at each of the fundamental frequencies that can be estimated.

A plot of I(ω) against ω is called the periodogram of y_{t}.

Here is an example.

> par(mfrow=c(4,1))
> tsplot(uk.male)
> tsplot(ts(Mod(fft(uk.male))[2:36]))
> title("Really raw periodogram")
> spectrum(uk.male)
> spectrum(uk.male,spans=c(3,5))

Notice that the main power corresponds to a frequency of 12 months, but there is also significant power at some other frequencies (like a two-year cycle). The graphs are a timeseries plot of the raw data, a (raw, unsmoothed) periodogram, an estimate of the spectrum (on a 10*log 10) scale, and an even more smoothed spectral estimate.

Now remember the lutenizing hormone example. We will look at the late1 series from these data.

> tsplot(late1)
> tsplot(ts(Mod(fft(late1))[2:24]))
> title("Really raw periodogram")
> spectrum(late1)
> tsplot( ts(10^(spectrum(late1, plot=F)$spec/10)))
> title("Spectral estimate on natural scale")

235

Partial Autocorrelation Function

Recall, in the dim past, that we saw the autocorrelation function (acf) for an AR(1) model looked just like an exponential decline (or an oscillatory exponential decline). What does the ACF of an AR(2) or AR(3) process look like? The following graph shows plots of two AR processes,

y_{t} = .7y_{t−1} + .2y_{t−2} + z_{t}

and

y_{t} = .7y_{t−1} + z_{t}

and their autocorrelation functions. Notice that the ACF’s have long lags and it is not easy to diagnose anything.

> tsplot(dd <- arima.sim(200, model=list(ar=c(.7, .2))))
> acf(dd)
> tsplot(dd <- arima.sim(200, model=list(ar=c(.7))))
> acf(dd)
>

What about a similar set of plots for a MA(1) and MA(2) process?

> tsplot(dd <- arima.sim(200, model=list(ma=c(.6))))

> acf(dd)
> tsplot(dd <- arima.sim(200, model=list(ma=c(.6, .4))))
> acf(dd)
>

The two models that we fit are

y_{t} = .6z_{t−1} + z_{t}

and

y_{t} = .6z_{t−1} + .4z_{t−2} + z_{t}

Notice that the ACF tails off after lag 1 and 2. That is the ACF provides a diagnostic of the order of the MA process. A simple (too simple!) conclusion is that if the ACF is zero beyond a certain lag, then that is evidence for an MA process of the appropriate order.

That turns out to be a useful observation, but it doesn’t help us with AR processes. We need to search for a diagnostic that will help us determine the order of an AR process—but we don’t have to search all that far. What we need is the partial autocorrelation function (PACF). The conceptual idea (not computational, but that doesn’t matter) is that we take a time series and fit AR models of orders 1, 2, 3, and so on. For the PACF at lag 1 we use the coefficient α_{1} of the fitted AR(1) model. At lag 2 we use the coefficient α_{2} of the AR(2) model, and so on. Clearly (I hope it is clear) if we really have an AR(k) model then all of the coefficients at lags greater than k should be zero—or at least they should be as close to zero as sampling error allows. It works really well. To actually calculate the PACF (and plot it in S) requires the use of the acf function with an extra argument.

acf(dd,type="p")

I have produced four plots, two of AR(2) processes and two of AR(1) processes. Notice that these plots look like the ACF plots of MA processes.

> acf(arima.sim(200, model=list(ar=c(.7, .2))), type="p")
> acf(arima.sim(200, model=list(ar=c(.7, .2))), type="p")
> acf(arima.sim(200, model=list(ar=c(.7))), type="p")
> acf(arima.sim(200, model=list(ar=c(.7))), type="p")

Maybe the PACF of a MA process will look like the ACF of an AR process? This sort of duality would be quite useful. Lets look and see. The simulations were generated by

> acf( arima.sim(200, model=list(ma=c(-.9))), type="p")
> acf( arima.sim(200, model=list(ma=c(.9))), type="p")
> acf( arima.sim(200, model=list(ma=c(.9, .6))), type="p")
> acf( arima.sim(200, model=list(ma=c(.9, .6))), type="p")

The following plots really do look like the ACF of AR’s. There is a mathematical duality between ACF’s and PACF’s, but the important thing to note is that each of them is quite good at helping to diagnose the order of a pure MA or a pure AR process.

Examples of ARIMA fitting

The purpose of these notes is to take a set of data (river data from the trinity river) and attempt a first pass at ARIMA modelling. The tools we will use are arima.mle (to do the mle estimation), arima.diag (to do the diagnostic checking), acf (to get ACF and PACF plots), diff (to do differences of the series), and tsplot(to plot things). In order to help a little, I’ll include the help pages for three of these commands here

Estimate    Autocovariance,    Autocorrelation     or     Partial
Autocorrelation

DESCRIPTION:
Estimates and displays autocovariance, autocorrelation  or
partial autocorrelation functions.
USAGE:
acf(x, lag.max=<<see below>>, type="correlation", plot=T)
REQUIRED ARGUMENTS:
x:     a univariate or multivariate time series or  a  vector  or
matrix.  Missing  values are allowed only at the beginning
or end of  series.  If  a  matrix,  rows  are  treated  as
timepoints and columns as univariate series.
OPTIONAL ARGUMENTS:
lag.max:    the maximum number of lags at which to  estimate  the
autocovariance.   If  this is not supplied, it is a number
proportional to the logarithm of the length of the series.
type:      a  character  string:  "covariance"  to  estimate  the
autocovariance    function,    "correlation"    for    the
autocorrelation function, or  "partial",  if  the  partial
autocorrelation  function is desired.  The start of one of
the strings will suffice.
plot:     If TRUE, the autocovariance or autocorrelation function
between  pairs  of univariate series will be plotted in an
array of at most 5 by  5  plots  per  page.   If  type  is
"correlation"  or  "partial"  approximate  95%  confidence
limits are drawn on the plots.
VALUE:
a list with the following components:
acf:   a three-dimensional array containing the autocovariance or
autocorrelation  function  estimates.   acf[i,j,k]  is the
covariance (or correlation) between  the  j-th  series  at
time t and the k-th series at time t+1-i.
lag:   an array the same shape as acf  containing  the  lags  (as

fractions   of  the  sampling  period)  at  which  acf  is
calculated.  If j > k  and  i  >  1,  then  lag[i,j,k]  is
negative.
n.used:   the number of observations in which no  missing  values
occur.
type:     a character string indicating  the  type  of  function,
"covariance", "correlation" or "partial".
series:   the name of x, including transformations.
SIDE EFFECTS:
if plot is TRUE, a high density  plot  (type="h")  of  the
resulting  function  will  be  produced  including the 95%
confidence limits about zero.
DETAILS:
The autocovariance function is estimated  by  summing  the
lagged  products and dividing by the length of the series.
For the  autocorrelation  function,  all  covariances  are
further divided by the geometric mean of the corresponding
variances.  For the partial autocorrelation function,  the
Levinson-Durbin recursion is used to fit AR(p) models to x
successively for p = 1, ..., lag.max,  and  from  the  AR-
coefficients  the  partial autocorrelation function values
are derived.
BACKGROUND:
Both  the  autocorrelation  function  (acf)  and   partial
autocorrelation  function (pacf) are useful tools in ARIMA
model   identification.    The   autocovariance   function
describes  the  serial  dependence  structure  of  a  time
series.  The  time  series  is  generally  assumed  to  be
stationary before these tools are used for Box and Jenkins
modeling purposes.
REFERENCES:
The Time Series chapter of the S-PLUS User's Manual.
SEE ALSO:
acf.plot, ar, arima.mle.
EXAMPLES:
acf.lynx <- acf(lynx, 36, "correlation")
pacf.lynx <- acf(lynx, 36, "partial")

ARIMA Modeling via Gaussian Maximum Likelihood

DESCRIPTION:
Returns a  list  representing  a  univariate  ARIMA  model
estimated by Gaussian maximum likelihood.
USAGE:
arima.mle(x, model, n.cond=<<see below>>, xreg=NULL,  ...)
REQUIRED ARGUMENTS:
x:     a univariate time series  or  a  vector.   Missing  values
(NAs) are allowed.
model:    a list specifying an ARIMA model.  If  a  simple  ARIMA
model  is specified (i.e., only one component), then model
should be a list with component names from: order, period,
ar,  ma,  ndiff,  ar.opt,  ma.trans,  and ma.opt.  See the
DETAILS section for further details.
OPTIONAL ARGUMENTS:
n.cond:   the number of observations on which  to  condition  the
likelihood.   n.cond  must  be at least P+D (the default),
where  P  and  D  are   the   orders   of   the   expanded
autoregressive and differencing polynomials.
xreg:     a univariate or multivariate time series or a vector or
matrix with univariate time series per column.  These will
be used as additive regression variables.
...:   additional arguments can be passed to nlmin.
VALUE:
a list representing the result of the fitted model:
model:     the  same  as  the  input  model  with  the  estimated
coefficients substituted in the components ar and ma.
var.coef:   the variance-covariance matrix for the autoregressive
and moving average coefficients.
loglik:   -2 times the log likelihood  of  the  model  (up  to  a
constant factor).
aic:   Akaike's information criteria,  which  is  loglik  plus  2
times the number of parameters fit.
sigma2:   the estimated innovations variance.
n.used:     the  number  of  observations  used  to  compute  the
likelihood.
n.cond:   the number of observations on which the  likelihood  is
conditioned.
reg.coef:   the estimated regression  coefficients  (returned  if
xreg is provided).
converged:      if  the  optimizer  has  apparently  successfully

converged  to a minimum, then converged is TRUE; otherwise
converged is FALSE.
conv.type:       a  character  string  describing  the  type   of
convergence.
series:   a character string giving  the  name  of  x,  including
transformations.
reg.series:    a  character  string  giving  the  name  of  xreg,
including   transformations  (only  returned  if  xreg  is
provided).
DETAILS:
SPECIFYING A MODEL.  The components of the model  argument
should  be  some or all of:  order, period, ar, ma, ndiff,
ar.opt, ma.trans, and ma.opt.  order is a vector of length
3 specifying the order of an ARIMA(p,d,q); the elements of
the vector are (p, d, q), where p and q are the  order  of
the  autoregressive  and moving average operators and d is
the number of differences.  period specifies the period of
the ARIMA operators; e.g., a seasonal operator for monthly
data would have a period of 12 (the default is 1).  ar and
ma  are  vectors  of initial values for autoregressive and
moving average coefficients respectively (if not provided,
the  initial values to the optimizer are set to 0).  ndiff
is  an  alternative  way  to   specify   the   number   of
differences.   A model can be specified using either order
or ar, ma, and ndiff, or both (but they must agree if both
are  specified).  The model doesn't explicitly allow for a
non-zero mean of the  series,  though  the  estimation  is
invariant to the mean of the data for many of the possible
models; if you estimate an AR model, you must subtract the
mean  yourself  or  use  the xreg argument in order to get
sensible estimates.
If ma.trans is TRUE  (the  default),  the  moving  average
coefficients  will  be  transformed before passing them to
the  optimizer,  ensuring  invertibility  of  the   model.
ar.opt  and  ma.opt  are logical vectors of length p and q
respectively.  If  the  ith  element  is  TRUE,  then  the
optimizer  will  optimize  over  the ith autoregressive or
moving average coefficient.  By default, the  vectors  are
TRUE.   This  option is useful to fit models in which some
low order coefficients are set to 0.
If a multiplicative ARIMA model is specified  (i.e.,  more
than  one component), then model is a list of lists.  Each
list represents a component of  the  multiplicative  ARIMA
model, and is specified in the same manner as above (i.e.,

the components of each list  are  ar,  ma,  ndiff,  order,
period, ma.trans, ar.opt, ma.opt).
THE ALGORITHM.  The likelihood is  conditioned  on  n.cond
observations.   If  no  missing  values  are  present, the
likelihood   is   computed   by   forming   the   Choleski
decomposition of the covariance matrix of the process (see
Ansley, 1979).  If missing values are  present,  then  the
likelihood  is computed using the Kalman filter applied to
a state space representation of the likelihood.  The state
space  representation is the same as that used by Kohn and
Ansley (1986).  The  method  of  initializing  the  Kalman
filter  recursions  is  that  given  by  Bell  and Hillmer
(1987).  The regression parameters are concentrated out of
the  likelihood,  as  in  Kohn  and  Ansley  (1985).   The
autoregressive   and   moving   average   parameters   are
transformed to ensure stationarity and invertibility as in
Jones (1980).
WARNING:
Unlike the ar function, the mean of the series will not be
estimated  by  arima.mle unless you use the xreg argument.
arima.mle assumes a zero mean series.
REFERENCES:
Ansley,  C.  F.  (1979).   An  algorithm  for  the   exact
likelihood   of   a  mixed  autoregressive-moving  average
process.  Biometrika 66, 59-65.
Bell, W. and Hillmer, S. (1987).  Initializing the  Kalman
filter   in   the  nonstationary  case.   Research  Report
CENSUS/SRC/RR-87/33, Statistical Research Division, Bureau
of the Census, Washington, DC, 20233.
Jones, R. H. (1980).  Maximum likelihood fitting  of  ARMA
models   to   time   series   with  missing  observations.
Technometrics 22, 389-395.
Kohn, R. and Ansley, C. F. (1985).   Efficient  estimation
and   prediction   in   time   series  regression  models.
Biometrika 72, 694-697.
Kohn,  R.  and  Ansley,   C.   F.   (1986).    Estimation,
prediction,   and  interpolation  for  ARIMA  models  with
missing  data.   Journal  of  the   American   Statistical
Association 81, 751-761.

The Time Series chapter of the S-PLUS User's Manual.
SEE ALSO:
ar, arima.diag, arima.filt, arima.forecast, nlmin, sabl.
EXAMPLES:
# Simulate an MA(2) series and fit the series using Gaussian
# Maximum Likelihood
ma <- arima.sim(model=list(ma=c(-.5, -.25)))
ma.fit <- arima.mle(ma, model=list(ma=c(-.5, -.25)))
# Fit a Box-Jenkins (0,1,1)x(0,1,1)12 Airline model to the ship data
# Use zeros as the starting values for the optimizer
al.mod <- list(list(order=c(0,1,1)), list(order=c(0,1,1), period=12))
fit <- arima.mle(ship, model=al.mod)

Compute Diagnostics for ARIMA Model

DESCRIPTION:
Computes diagnostics for an ARIMA model.  The  diagnostics
include the autocorrelation function of the residuals, the
standardized residuals, and the  portmanteau  goodness  of
fit test statistic.
USAGE:
arima.diag(z, acf.resid=T, gof.lag=10, lag.max=<<see below>>, resid=F,
           std.resid=T, plot=T, ...)
REQUIRED ARGUMENTS:
z:     a list like the output from arima.mle.
OPTIONAL ARGUMENTS:
acf.resid:     logical flag: if TRUE, the autocorrelation of  the
residuals will be returned.
lag.max:    the maximum number of lags at which to  estimate  the
autocovariance.   If  this  is  not  supplied,  it  is the
maximum  between  gof.lag  plus  the   number   of   model
parameters  and  a number proportional to the logarithm of
the length of the series.
gof.lag:    if gof.lag>0, then gof.lag plus the number  of  model
parameters  is the number of lags to use for computing the
Portmanteau goodness of fit statistic.  If gof.lag=0, then
the statistic will not be computed

resid:    logical flag: if  TRUE,  then  the  residuals  will  be
returned.
std.resid:     logical  flag:  if  TRUE,  then  the  standardized
residuals will be returned.
plot:     logical flag: if TRUE, the diagnostics will be  plotted
using the function arima.diag.plot
...:    additional  arguments  may  be  passed  to  the  function
arima.diag.plot
VALUE:
a list (which is returned invisibly when  plot=TRUE)  with
the following elements:
acf.list:   a list representing the autocorrelation  function  of
the residuals.  See acf for details.
gof:   a  list  representing  the  Portmanteau  goodness  of  fit
statistics  computed  for  a  range of lags.  The list has
four elements: lag, statistic,  df,  p.value.   lag  is  a
vector   of  the  number  of  lags  used  to  compute  the
statistics.   statistic  is  the  vector   of   statistics
corresponding  to  each  lag  used.   df  is the number of
degrees of freedom the test statistics have under the null
hypothesis that the model is correct.  p.value is a vector
of the p-values for the  statistics  using  a  Chi-Squared
distribution with the appropriate degrees of freedom.
resid:    the residuals or innovations for the process.
std.resid:     the standardized  residuals.   The  residuals  are
standardized  to  have  unit variance under the assumption
that the model is correct and the process is Gaussian.
series:   the name of x, including transformations.
SIDE EFFECTS:
if plot is TRUE, the diagnostics will be plotted using the
function arima.diag.plot
DETAILS:
The residuals (both standardized  and  raw)  are  computed
using   the   function  arima.filt.   The  autocorrelation
function of the residuals is computed using  the  function
acf.   The  portmanteau test statistic is derived from the
autocorrelation function of the residuals  (see  the  Time
Series chapter in the S-PLUS User's Manual for details).
REFERENCES:
Box, G. E. P. and Jenkins,  G.  M.  (1976).   Time  Series
Analysis:  Forecasting  and Control.  Holden-Day, Oakland,
Calif.  Chapter 8.

       The Time Series chapter of the S-PLUS User's Manual.
SEE ALSO:
       arima.diag.plot arima.filt, arima.mle.
EXAMPLES:
       # compute and plot diagnostics for simulated AR(1) series with mean 5
       x <- arima.sim(model=list(ar=.9)) + 5
       xreg <- rep(1,100)
       fit <- arima.mle(x,model=list(ar=.9), xreg=xreg)
       diag <- arima.diag(fit)
       lynx.arma11 <- arima.mle(lynx, model=list(ar=0, ma=0))
       arima.diag(lynx.arma11)

The data we with to look at have the following first few lines

TRINITY RIVER AT LEWISTON, CALIFORNIA, OCTOBER 1912 - SEPTEMBER 1960
83.42 100.71 86.5 106.15 107.77 78.79
85.67 110.47 107.49 76.07 90.48 105.36
74.58 99.29 89.93 113.55 101.06 59.08
97.34 105.89 108.67 83.38 82.43 104.72

We read in the data as (with skip to skip the first line of rubbish).

trinity <- as.ts( scan("trinity.dat", skip=1))

The first thing we want to do is a number of plots of the data, and acf, pacf plots. Notice that the seasonality is easier to see in the partial data (second plot) and that the 12 month effects show up strongly in the ACF and PACF

> tsplot(trinity)
> title("tsplot(trinity)")
> tsplot(trinity[1:96])
> title("tsplot(trinity[1:96])")
> acf(trinity)
> acf(trinity, type="p")

We have a lot of reason to think about 12th order difference of the data, so lets try that.

> delme <- diff(trinity, lag=12)
> tsplot(delme)
> title("tsplot(delme)")
> tsplot(delme[1:96])
> title("tsplot(delme[1:96])")
> acf(delme)
> acf(delme, type="p")

Those plots now look quite good. One might temporarily guess that a MA(2) model might do the trick, or maybe an ARMA(1,1). Lets try those and see where we get. First fit the model that has 12th order differences, and is otherwise an MA(2) model.

>  mymod <- arima.mle(trinity, model=list( list(ndiff=1, period=12),
+ list(order=c(0,0,2))))
> mymod$loglik
[1] 4271.951

Also try the same except an ARIMA(1,0,1) model.

>  mymod2 <- arima.mle(trinity, model=list( list(ndiff=1, period=12),
+ list(order=c(1,0,1))))
 > mymod2$loglik
[1] 4256.779

Just for kicks we’ll try the raw AR(1) and raw MA(1) and see how they do.

> arima.mle(trinity, model=list( list(ndiff=1, period=12),
+ list(order=c(1,0,0))))$loglik
[1] 4307.937

> arima.mle(trinity, model=list( list(ndiff=1, period=12),
+ list(order=c(0,0,1))))$loglik
[1] 4305.263

Either of these is demonstrably works (likelihood ratio test) than the combined models. Lets quickly look at the diagnostics for the two candidate models.

> arima.diag(mymod)
> arima.diag(mymod2)

I could not not figure out how to automatically get the pacfs, so here is the manual approach to that

> resid2 <- arima.diag(mymod2,std=T, plot=F)$std.resid
> resid <- arima.diag(mymod,std=T, plot=F)$std.resid
> acf(resid,type="p")
> acf(resid2,type="p")

The plots are very similar, but I (again) prefer the second one. The summary of our fit is

$model:
$model[[1]]:
$model[[1]]$ndiff:
[1] 1
$model[[1]]$period:
[1] 12
$model[[1]]$order:
[1] 0 1 0

$model[[2]]:
$model[[2]]$order:
[1] 1 0 1
$model[[2]]$ar:
[1] 0.5485699
$model[[2]]$ndiff:
[1] 0
$model[[2]]$ma:
[1] 0.8451971

251 252