Department of Statistics, Virginia Tech

This homework is due on **Thursday, October 26th at 2pm** (the start of class). Please turn in all your work. This homework primarily covers \(r \times c\) contingency tables, median test, goodness-of-fit test, contingency coefficients and Cochran’s test.

*Three professors are teaching large classes in introductory statistics. At the end of the semester, they compare grades to see if there are significant differences in their grading policies. Are these differences significant?*

Prof | A | B | C | D | F | WP | WF |
---|---|---|---|---|---|---|---|

Smith |
12 | 45 | 49 | 6 | 13 | 18 | 2 |

Jones |
10 | 32 | 43 | 18 | 4 | 12 | 6 |

White |
15 | 19 | 32 | 20 | 6 | 9 | 7 |

*(3 pts) Which test are you using?*

Pearson’s Chi-square test for difference in probabilities.

*(7 pts) Perform the appropriate test “by hand”.*

The hypotheses are

\[ \begin{aligned} \mathcal{H}_0 &: \; \mbox{All of the probabilities for the same grade are equal to each other} \\ \mathcal{H}_1 &: \; \mbox{All of the probabilities for the same grade are not equal to each other} \end{aligned} \]

First, lets collect the essence of the table in R.

```
O <- rbind(smith=c(12, 45, 49, 6, 13, 18, 2),
jones=c(10, 32, 43, 18, 4, 12, 6),
white=c(15, 19, 32, 20, 6, 9, 7))
colnames(O) <- c("A", "B", "C", "D", "F", "WP", "WF")
```

Now we can calculate column and row sums, etc.

```
n <- rowSums(O)
C <- colSums(O)
N <- sum(n)
```

And we can flesh out the table with those values.

```
tab <- rbind(O, C)
tab <- cbind(tab, c(n, N))
rownames(tab)[nrow(tab)] <- colnames(tab)[ncol(tab)] <- "tot"
tab
```

```
## A B C D F WP WF tot
## smith 12 45 49 6 13 18 2 145
## jones 10 32 43 18 4 12 6 125
## white 15 19 32 20 6 9 7 108
## tot 37 96 124 44 23 39 15 378
```

Now, calculate the expected counts under the null hypothesis.

`E <- outer(n, C/N)`

And then the test statistic.

```
t <- sum((O - E)^2/E)
t
```

`## [1] 28.91509`

Finally, the \(p\)-value is

```
df <- (length(n)-1)*(length(C)-1)
df
```

`## [1] 12`

`pchisq(t, df, lower.tail=FALSE)`

`## [1] 0.004055965`

- An easy reject of \(\mathcal{H}_0\) at the 5% level. We conclude that there are significant differences in their grading policies.

*(5 pts) Perform the test in (ii) “using a software library”.*

`chisq.test(O, correct=FALSE)`

```
## Warning in chisq.test(O, correct = FALSE): Chi-squared approximation may be
## incorrect
```

```
##
## Pearson's Chi-squared test
##
## data: O
## X-squared = 28.915, df = 12, p-value = 0.004056
```

- Same result; never-mind the warning.

*A random sample of 30 graduating university seniors was categorized by college and religious preference as follows. Is there a relationship between college and religious preference?*

Arts & Sciences | Business | Engineering | Other |
---|---|---|---|

Catholic | Protestant | Catholic | Protestant |

Catholic | Catholic | Protestant | Catholic |

Jewish | Jewish | Protestant | Catholic |

Catholic | Jewish | Protestant | Other |

Protestant | Protestant | Other | Catholic |

Protestant | Protestant | Other | |

Other | Other | Catholic | |

Protestant | Jewish | ||

Other | Other |

*(3 pts) Which test are you using?*

Chi-squared test for independence.

*(7 pts) Perform the appropriate test “by hand”.*

Here we will test the following hypotheses.

\[ \begin{aligned} \mathcal{H}_0: & \; \mbox{Religion is independent of college affiliation} \\ \mathcal{H}_1: & \; \mbox{Religion depends, in part, on college affiliation} \end{aligned} \]

First we have to build the contingency table, which I will do below directly in R.

```
O <- rbind(catholic=c(3, 1, 1, 4),
jewish=c(1, 2, 0, 1),
protestant=c(3, 3, 3, 1),
other=c(2, 1, 1, 3))
colnames(O) <- c("art & science", "business", "engineering", "other")
R <- rowSums(O)
C <- colSums(O)
N <- sum(R)
tab <- rbind(O, C)
tab <- cbind(tab, c(R, N))
rownames(tab)[nrow(tab)] <- colnames(tab)[ncol(tab)] <- "total"
tab
```

```
## art & science business engineering other total
## catholic 3 1 1 4 9
## jewish 1 2 0 1 4
## protestant 3 3 3 1 10
## other 2 1 1 3 7
## total 9 7 5 9 30
```

Now we can calculate the expected counts.

`E <- outer(R, C/N)`

Then the observed value of the \(T\) statistic can be calculated as follows.

`t <- sum((O-E)^2/E)`

And the \(p\)-value can be calculated as.

```
df <- (length(R)-1)*(length(C) - 1)
pchisq(t, df, lower.tail=FALSE)
```

`## [1] 0.6781851`

- We fail to reject \(\mathcal{H}_0\) at the 5% level and conclude that religious preference and college are independent.

*(5 pts) Perform the test in (ii) “using a software library”.*

`chisq.test(O, correct=FALSE)`

```
## Warning in chisq.test(O, correct = FALSE): Chi-squared approximation may be
## incorrect
```

```
##
## Pearson's Chi-squared test
##
## data: O
## X-squared = 6.6048, df = 9, p-value = 0.6782
```

- Same result; never-mind the warning.

*An experiment was conducted to determine the compression strength (lb) for different types of boxes. We are interested in determining if the median is the same for all types of boxes.*

Box 1 |
655.5 | 788.3 | 734.3 | 721.4 |

Box 2 |
789.2 | 772.5 | 786.9 | 686.1 |

Box 3 |
737.1 | 539.0 | 696.3 | 671.7 |

Box 4 |
535.1 | 628.7 | 542.4 | 559.0 |

*(10 pts) Perform the test “by hand” including pairwise comparisons if appropriate.*

We will be performing the median test with the following hypotheses.

\[ \begin{aligned} \mathcal{H}_0 &: \; \mbox{All boxes have the same median.} \\ \mathcal{H}_1 &: \; \mbox{At least two boxes have different medians.} \end{aligned} \]

Lets create the lists in R.

```
boxes <- list(box1=c(655.5, 788.3, 734.3, 721.4),
box2=c(789.2, 772.5, 786.9, 686.1),
box3=c(737.1, 539.0, 696.3, 671.7),
box4=c(535.1, 628.7, 542.4, 559.0))
```

First we extract the lengths of the lists

```
n <- sapply(boxes, length)
n
```

```
## box1 box2 box3 box4
## 4 4 4 4
```

and the grand median.

```
grandmed <- median(unlist(boxes))
grandmed
```

`## [1] 691.2`

With that we can create the table. First lets calculate the “greater” row and then the “lesser” row.

```
Og <- sapply(boxes, function(x) { sum(x > grandmed) } )
Oleq <- n - Og
O <- rbind(Og, Oleq)
rownames(O) <- c("> Med", "<= Med")
O
```

```
## box1 box2 box3 box4
## > Med 3 3 2 0
## <= Med 1 1 2 4
```

Then we can calculate the row sums.

```
a <- sum(Og)
b <- sum(Oleq)
N <- sum(n)
```

Now we can dress up the table for inspection, although we don’t need it for anything specifically.

```
tab <- rbind(O, n)
tab <- cbind(tab, c(a, b, N))
rownames(tab)[nrow(tab)] <- colnames(tab)[ncol(tab)] <- "tot"
tab
```

```
## box1 box2 box3 box4 tot
## > Med 3 3 2 0 8
## <= Med 1 1 2 4 8
## tot 4 4 4 4 16
```

The test statistic may be calculated as follows.

```
t <- (N^2 / (a*b)) * sum((Og - n*a/N)^2/n)
t
```

`## [1] 6`

Finally, the \(p\)-value is …

```
df <- length(boxes) - 1
pchisq(t, df, lower.tail=FALSE)
```

`## [1] 0.1116102`

- We fail to reject the null hypothesis at the 5% level, and conclude that all boxes have the same median.

*(5 pts) Perform the test “using a software library”.*

Here we will borrow the library we created in class.

```
median.test <- function(x)
{
## summarize lengths of populations
n <- sapply(x, length)
## calculate grand mean and build observation table
grandmed <- median(unlist(x))
Og <- sapply(x, function(x) { sum(x > grandmed) } )
Oleq <- n - Og
O <- rbind(Og, Oleq)
## calculate totals columns of table
a <- sum(Og)
b <- sum(Oleq)
N <- sum(n)
## build the table for visualization
tab <- rbind(Og, Oleq, n)
tab <- cbind(tab, c(a, b, N))
rownames(tab)[nrow(tab)] <- colnames(tab)[ncol(tab)] <- "tot"
## calculate the test statistic, null distribution and p-value
t <- (N^2 / (a*b)) * sum((Og - n*a/N)^2/n)
df <- length(x) - 1
phi <- pchisq(t, df, lower.tail=FALSE)
## return info
return(list(grand.med=grandmed, table=tab, t.stat=t, p.value=phi))
}
```

Here it is on our `boxes`

data.

`median.test(boxes)`

```
## $grand.med
## [1] 691.2
##
## $table
## box1 box2 box3 box4 tot
## Og 3 3 2 0 8
## Oleq 1 1 2 4 8
## tot 4 4 4 4 16
##
## $t.stat
## [1] 6
##
## $p.value
## [1] 0.1116102
```

- Same result; excellent!

*A random sample of 30 stocks was selected from each of the three major U.S exchanges, and their performance over the previous year was monitored. The median performance for all 90 stocks was noted and the following table constructed:*

Exchange | Above Median |
---|---|

New York | 18 |

American | 17 |

NASDAQ | 10 |

*Was there a significant difference in the performance of stocks on the three exchanges during the previous year?*

*(8 pts) Perform the test “by hand”.*

Again, we are performing the median test with the same hypotheses as above.

This is a little easier because we have some of the table already populated for us.

```
n <- rep(30, 3)
Og <- c(18, 17, 10)
Oleq <- n - Og
O <- rbind(Og, Oleq)
colnames(O) <- c("NY", "American", "NASDAQ")
rownames(O) <- c("> Med", "<= Med")
O
```

```
## NY American NASDAQ
## > Med 18 17 10
## <= Med 12 13 20
```

Now we can dress up the table for inspection, although we don’t need it for anything specifically.

```
tab <- rbind(O, n)
tab <- cbind(tab, c(a, b, N))
rownames(tab)[nrow(tab)] <- colnames(tab)[ncol(tab)] <- "tot"
tab
```

```
## NY American NASDAQ tot
## > Med 18 17 10 8
## <= Med 12 13 20 8
## tot 30 30 30 16
```

The test statistic may be calculated as follows.

```
t <- (N^2 / (a*b)) * sum((Og - n*a/N)^2/n)
t
```

`## [1] 5.066667`

Finally, the \(p\)-value is …

```
df <- length(n) - 1
pchisq(t, df, lower.tail=FALSE)
```

`## [1] 0.07939393`

- Fail to reject the null at the 5% level; we conclude that all exchanges have the same median.

*(7 pts) Perform the test “using a software library”.*

We have two choices here: (a) we could try to use our library function above, but we’ll need to fabricate some original data with an appropriate list of values; or (b), we can modify the function to work only with the counts, rather than directly with the list.

For (a) consider the following.

```
grandmed <- 0
exchange <- list(NY=c(rep(1,Og[1]), rep(-1,Oleq[1])),
American=c(rep(1,Og[2]), rep(-1,Oleq[2])),
NASDAQ=c(rep(1,Og[3]), rep(-1,Oleq[3])))
median.test(exchange)
```

```
## $grand.med
## [1] 0
##
## $table
## NY American NASDAQ tot
## Og 18 17 10 45
## Oleq 12 13 20 45
## tot 30 30 30 90
##
## $t.stat
## [1] 5.066667
##
## $p.value
## [1] 0.07939393
```

- Bingo.

Now for (b), consider the following new library function.

```
median.test.counts <- function(O)
{
## summarize lengths of populations
n <- colSums(O)
## calculate totals columns of table
a <- sum(Og)
b <- sum(Oleq)
N <- sum(n)
## build the table for visualization
tab <- rbind(Og, Oleq, n)
tab <- cbind(tab, c(a, b, N))
rownames(tab)[nrow(tab)] <- colnames(tab)[ncol(tab)] <- "tot"
## calculate the test statistic, null distribution and p-value
t <- (N^2 / (a*b)) * sum((Og - n*a/N)^2/n)
df <- ncol(O) - 1
phi <- pchisq(t, df, lower.tail=FALSE)
## return info
return(list(table=tab, t.stat=t, p.value=phi))
}
```

Lets try it.

`median.test.counts(O)`

```
## $table
## NY American NASDAQ tot
## Og 18 17 10 45
## Oleq 12 13 20 45
## tot 30 30 30 90
##
## $t.stat
## [1] 5.066667
##
## $p.value
## [1] 0.07939393
```

- Perfect.

*A certain type of insect that is found in lakes in the southwestern United States is studied to see if the chromosomal structure is significantly different among states. The number of insects of various chromosomal types is recorded as follows:*

Type | Texas | New Mexico | Arizona | California |
---|---|---|---|---|

A |
54 | 72 | 83 | 96 |

B |
20 | 6 | 18 | 6 |

C |
17 | 3 | 12 | 0 |

D |
0 | 12 | 14 | 1 |

E |
0 | 10 | 0 | 0 |

*(10 pts) Compute the following “by hand”:**(3 pts) \(T\).*

To compute \(T\) we will need to complete the totals column in the table.

First, lets tabulate the data values in R.

```
O <- rbind(A=c(54, 72, 83, 96),
B=c(20, 6, 18, 6),
C=c(17, 3, 12, 0),
D=c(0, 12, 14, 1),
E=c(0, 10, 0, 0))
colnames(O) <- c("Texas", "New Mexico", "Arizona", "California")
```

Now we can compute the row an column totals.

```
r <- rowSums(O)
c <- colSums(O)
N <- sum(r)
```

Then we can fill out the table.

```
tab <- cbind(O, r)
tab <- rbind(tab, c(c, N))
rownames(tab)[nrow(tab)] <- colnames(tab)[ncol(tab)] <- "tot"
tab
```

```
## Texas New Mexico Arizona California tot
## A 54 72 83 96 305
## B 20 6 18 6 50
## C 17 3 12 0 32
## D 0 12 14 1 27
## E 0 10 0 0 10
## tot 91 103 127 103 424
```

The expected counts are …

`E <- outer(r, c/N)`

… so the observed value of the test statistic is

```
t <- sum((O-E)^2/E)
t
```

`## [1] 100.913`

*(2 pts) \(R_1\).*

```
r1 <- t/(N*(min(length(r),length(c))-1))
r1
```

`## [1] 0.07933409`

*(1 pts) Cramer’s contingency coefficient.*

```
Ccc <- sqrt(r1)
Ccc
```

`## [1] 0.2816631`

*(2 pts) Pearson’s contingency coefficient.*

```
r2 <- sqrt(t/(N+t))
r2
```

`## [1] 0.4384598`

*(2 pts) Pearson’s Mean-Square contingency coefficient.*

```
r3 <- t/N
r3
```

`## [1] 0.2380023`

*(5 pts) Find the available coefficients “using a software library”.*

`library(vcd)`

`## Loading required package: grid`

`assocstats(O)`

```
## X^2 df P(> X^2)
## Likelihood Ratio 108.40 12 0.0000e+00
## Pearson 100.91 12 3.3307e-16
##
## Phi-Coefficient : NA
## Contingency Coeff.: 0.438
## Cramer's V : 0.282
```

A die was cast 600 times with the following results:

Occurrence |
1 | 2 | 3 | 4 | 5 | 6 |

Frequency |
87 | 96 | 108 | 89 | 122 | 98 |

Is the die balanced?

*(5 pts) Perform the test by hand.*

First, lets code the data in R.

```
O <- c(87, 96, 108, 89, 122, 98)
N <- sum(O)
N
```

`## [1] 600`

A balanced die would have a probability of \(1/6\) for each number on each of the six faces. Therefore, in a total number of counts of, 600, we would expect the following counts.

```
pstar <- rep(1/6, 6)
E <- N * pstar
E
```

`## [1] 100 100 100 100 100 100`

To settle the discrepancy we can perform a Goodness-of-fit test, with the null hypothesis that \(p_i* = 1/6\) for \(i=1,\dots, 6\), versus the alternative that at least two of those probabilities are different.

The test statistic is calculated as follows.

```
t <- sum((O-E)^2/E)
t
```

`## [1] 8.58`

And the \(p\)-value is

`pchisq(t, length(O)-1, lower.tail=FALSE)`

`## [1] 0.1270355`

- Therefore we fail to reject the null hypothesis at the 5% level and conclude that the die is fair.

*(5 pts) Perform the test using software.*

`chisq.test(O, p=pstar)`

```
##
## Chi-squared test for given probabilities
##
## data: O
## X-squared = 8.58, df = 5, p-value = 0.127
```

- Same result.

*In an attempt to compare the relative power of three statistical tests, 100 sets of artificial data were generated using a computer. On each set of data the three statistical tests were used with \(\alpha=0.05\), and the decision to accept or reject the null hypothesis was recored. The results are as follows:*

Test 1 | Test 2 | Test 3 | Numbers of Sets of Data |
---|---|---|---|

Accept | Accept | Accept | 26 |

Accept | Accept | Reject | 6 |

Accept | Reject | Accept | 12 |

Reject | Accept | Accept | 4 |

Reject | Reject | Accept | 18 |

Reject | Accept | Reject | 5 |

Accept | Reject | Reject | 7 |

Reject | Reject | Reject | 22 |

*Is there a difference in the power of the three tests when applied to populations from which the simulated data were obtained?*

*(3 pts) Which test are you using?*

Cochran’s test.

*(7 pts) Perform the appropriate test “by hand”.*

The trick here is to expand out each binary pattern above the requestie number of times, as indicated in the final column. The version below codes “Accept” as 1 and “Reject” as 0.

```
m1 <- matrix(rep(c(1,1,1),26), nrow=26, byrow=TRUE)
m2 <- matrix(rep(c(1,1,0),6), nrow=6, byrow=TRUE)
m3 <- matrix(rep(c(1,0,1),12), nrow=12, byrow=TRUE)
m4 <- matrix(rep(c(0,1,1),4), nrow=4, byrow=TRUE)
m5 <- matrix(rep(c(0,0,1),18), nrow=18, byrow=TRUE)
m6 <- matrix(rep(c(0,1,0),5), nrow=5, byrow=TRUE)
m7 <- matrix(rep(c(1,0,0),7), nrow=7, byrow=TRUE)
m8 <- matrix(rep(c(0,0,0),22), nrow=22, byrow=TRUE)
X <- rbind(m1, m2, m3, m4, m5, m6, m7, m8)
```

Now we can calculate the row and column sums.

```
R <- rowSums(X)
C <- colSums(X)
N <- sum(R)
```

I’ll forgo presenting the table because it is too big.

The observed value of the test statistic is

```
t <- length(C)*(length(C)-1)*(sum((C - N/length(C))^2))/sum(R*(length(C)- R))
t
```

`## [1] 10.42308`

Then we may calculate the \(p\)-value as follows.

`pchisq(t, length(C)-1, lower.tail=FALSE)`

`## [1] 0.005453278`

- Reject \(\mathcal{H}_0\) at the 5% level and conclude that the three test do not have the same statistical power.

- (5 pts) Perform the test in (ii) “using a software library”.

The library we built in class is pasted below

```
cochran.test <- function(X)
{
## build out table
R <- rowSums(X)
C <- colSums(X)
N <- sum(R)
tab <- rbind(X, C)
tab <- cbind(tab, c(R, N))
rownames(tab)[nrow(tab)] <- colnames(tab)[ncol(tab)] <- "tot"
## calculate test statistic
t <- length(C)*(length(C)-1)*(sum((C - N/length(C))^2))/sum(R*(length(C)- R))
## calculate the p-value
phi <- pchisq(t, length(C)-1, lower.tail=FALSE)
return(list(table=tab, t.stat=t, p.value=phi))
}
```

Here is the result on this data, showing just the test statistic and \(p\)-value, so that big table doesn’t get printed to the screen.

```
out <- cochran.test(X)
out$t.stat
```

`## [1] 10.42308`

`out$p.value`

`## [1] 0.005453278`

- Bingo!