Two-Way Tables (Reading: Faraway (2006), section 4.1~4.2)

Two-by-Two Tables.

The data shown below were collected as part of a quality improvement study at a semiconductor factory.

A sample of wafers was drawn and cross-classified according to whether a particle was found on the die that produced the wafer and whether the wafer was good or bad. First, let us set up the data in a convenient form for analysis:

> y <- c(320,14,80,36)
> particle <- gl(2,1,4,labels=c("no","yes"))
> quality <- gl(2,2,labels=c("good","bad"))
> wafer <- data.frame(y,particle,quality)
> wafer

We will need the data in this form with one observation per line for our model fitting, but usually we prefer to observe its table form:

> (ov <- xtabs(y ~ quality+particle))

The data might have arisen under several possible sampling schemes:

1. observe the manufacturing process for a certain period of time and observed 450 wafers Þ all the numbers are random Þ we could use a Poisson model

2. decide to sample 450 wafers Þ 450 is a fixed number and other numbers are random Þ we could use a multinomial model

3. select 400 wafers without particles and 50 wafers with particles Þ 400, 50, and 450 are fixed numbers and other numbers are random Þ we could use a binomial model

4. select 400 wafers without particles and 50 wafers with particles that also included, by design, 334 good wafers and 116 bad ones Þ 400, 50, 334, 116, and 450 are fixed numbers and other numbers are random Þ we could use a hypergeometric model

The main question of interest concerning these data is whether the presence of particles on the wafer affects the quality outcome. We shall see in the following that all four sampling schemes lead to exactly the same conclusion.

Sampling Scheme 1: Poisson Model.
Under sampling scheme 1, it would be natural to view these outcomes, Y_ij, 1£i£I, 1£j£J, occurring at different rates, lp_ij,  where l is an unknown value of a size variable and p_ij is the probability that a randomly selected subject falls in the ij^th cell. We can form a Poisson model for these rates, i.e., Y_ij~Poisson(lp_ij). Suppose that we fit a Poisson GLM with only main effect terms:

> modl <- glm(y ~ particle+quality, poisson)
> summary(modl)

Notice that:

• the null model, which suggests all four outcomes occur at the same rate, does not fit because the deviance of 474.10 is very large for 3 degrees of freedom

• the fitted model has a residual deviance of 54.03, which is clearly an improvement over the null model. An alternative goodness-of-fit measure is the Pearson X²:

> sum(residuals(modl,type="pearson")^2)

[1] 62.81231

• we might also want to test the significance of the individual predictors. The null of the tests corresponds to the hypothesis of equivalent marginal proportions, i.e., p₁₊=p₂₊=...=p_I+. or p₊₁=p₊₂=...=p_+J. For the table, we could use the z-values in the output, but it is better to use the difference-in-deviance test:

> drop1(modl,test="Chi")

Single term deletions

 

Model:

y ~ particle + quality

         Df Deviance    AIC    LRT   Pr(Chi)   

<none>         54.03  83.77                    

particle  1   363.91 391.66 309.88 < 2.2e-16 ***

quality   1   164.22 191.96 110.18 < 2.2e-16 ***

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We see that both predictors are significant relative to the fitted model. It suggests that the marginal proportions are not equivalent for both predictors.

• by examining the coefficients, we see that wafers without particles occur at a significantly higher rate than wafer with particles because the estimated coefficient of particleyes is negative. Similarly, we see that good-quality wafers occur at a significantly higher rate than bad-quality wafers because the estimated coefficient of qualitybad is negative

• the coefficients of a model containing only main effect terms are closely related to the marginal totals because the MLE of the rates m=(lp_ij) must satisfy X^Ty=X^T(m-hat), where X^Ty is, in this example:

  > (t(model.matrix(modl)) %*% y)[,]

  (Intercept) particleyes  qualitybad

          450          50         116

  So we see that the fitted values, m-hat, must be a function of marginal totals. The coding we use will determine exactly how the coefficients relate to the marginal totals. For example, consider the particle factor, which is coded as 0 and 1 for No Particle and Particle categories, respectively. Its fitted coefficient is:

  > coef(modl)[2]

  particleyes

    -2.079442

  > (coefparticle <- exp(c(0, coef(modl)[2])))

              particleyes

        1.000       0.125

  > coefparticle/sum(coefparticle)

              particleyes

    0.8888889   0.1111111

  We see that these are just the marginal proportions for the no particle wafers and particle wafers respectively:

  > c(400, 50)/450

  [1] 0.8888889 0.1111111

• this model has a deviance of 54.03 on one degree of freedom and so does not fit the data. A Poisson GLM with only main effect terms corresponds to independent contingency tables, i.e., p_ij=p_i+p_+j. Because the model does not fit, we reject the hypothesis that good- or bad-quality outcomes would occur independently of whether a particle was found on the wafer.

• the addition of an interaction term would saturate the model and so would have zero deviance and zero degrees of freedom. So, an hypothesis comparing the models with and without interaction would use a test statistic of 54.03 on one degree of freedom, which is exactly the goodness-of-fit test. The hypothesis of no interaction would be rejected --- the necessity to fit a Poisson GLM with interaction indicates that the two predictors are associated, i.e., dependent. 

• for each cells, the fitted values under the Poisson GLM with only main effect terms are:

  > xtabs(fitted(modl) ~ quality+particle)

         particle

  quality        no       yes

     good 296.88889  37.11111

     bad  103.11111  12.88889

  Notice that these fitted values form an independent contingency table because its odds-ratio is one:

  > (296.88889*12.88889)/(103.11111*37.11111)

  [1] 1.000000

Sampling Scheme 2: Multinomial Model.
Suppose that we assume that the total sample size was fixed at 450 and that the frequency of the four possible outcomes was recorded. In this circumstances, it is natural to use a multinomial distribution to model the response. There is a close connection between multinomial and Poisson models as follows. Let Y_ij be the independent Poisson random variables with means lp_ij, 1£i£I, 1£j£J, then the joint distribution of (Y₁₁, Y₁₂, ..., Y_IJ|S_ijY_ij=y₊₊) is multinomial(y₊₊, p₁₁, p₁₂, ..., p_IJ). This essentially means that we can model multinomial GLM's using a Poisson GLM. The inferences for p_ij in the multinomial model would coincide with that in the Poisson GLM. For example, let us compute the MLE of p_i+ and p_+j under the multinomial model, which are S_jY_ij/y₊₊ and S_iY_ij/y₊₊, respectively:  

One alternative to the deviance is the Pearson X² statistic S_iS_j[(O_ij-E_ij)²/E_ij]:

which is also the same as that in the Poisson model with only main effect terms. We can apply Yates' continuity correction on the Pearson X² statistic, which gives superior results for small samples. This correction is implemented in:

The X² statistic with continuity correction is 60.1239 for one degree of freedom, which clearly rejects the null, i.e., independence hypothesis.

Larger Two-Way Table.
Snee (1974) presents the data on 592 students cross-classified by hair and eye color. Let us read the data into R and take a look of it:

> haireye <- read.table("haireye.txt")
> haireye

The data is more conveniently displayed using:

> (ct <- xtabs(y ~ hair + eye, haireye))

(Q: what information have you found by merely looking at the table?)

We can execute the usual Pearson's X² test for independence as:

> summary(ct)

where

• we see that hair and eye color are clearly not independent because an X² of 138.29 on 9 degrees of freedom is way too large

• 9=(4-1)(4-1) is the degrees of freedom for the interaction between the two categorical variables.

Now, let us fit a Poisson GLM containing only main effect terms:

> modc <- glm(y ~ hair + eye, family=poisson, haireye)
> summary(modc)

Notice that:

• most of the levels of hair and eye color show up as significantly different from the reference levels of black hair and green eyes, which indicates that there are higher numbers of people with some hair colors than others and some eye colors than others, i.e., the marginal proportions of hair and eye color are not equivalent

• the deviance of 146.44 on 9 degrees of freedom shows that they are clearly dependent

• we can use the Pearson residuals under this Poisson model to calculate a Pearson X²:

> sum(residuals(modc,type="pearson")^2)

[1] 138.2898

We get the same value as the previous X².

One option for displaying contingency table data is a graphical tool, dotchart:

> dotchart(ct)

An alternative graphical tool is the mosaic plot, which divides the plot region according to the frequency of each level in a recursive manner:

> mosaicplot(ct,color=T,main="")

In the plot,

• the area is first divided according to the frequency of hair color

• within each hair color, the area is then divided according to the frequency of eye color

• we can easily read the information of p_i+ and p_j|i from the plot

• a different plot could be constructed by reversing the order of hair and eye in the xtabs command above

• we can now readily see the frequency of various outcomes. For example, brown hair and brown eyes is the most common combination while black hair and green eyes is the least common

(Q: what patterns will you see in the dotchart and mosaic plots if the two variables are independent? what if they have equivalent marginal proportions? Generate some tables with the two properties and draw the plots.)