# 2.5 - Examples in SAS/R: Dice Rolls & Tomato

Printer-friendly version

## Example - Dice Rolls

Suppose that we roll dice 30 times and observe the following table showing the number of times each face ends up on top.

 Face Count 1 3 2 7 3 5 4 10 5 2 6 3 Total 30

We want to test the null hypothesis that the dice is fair. Under this hypothesis, X ∼ Mult(n = 30, π0) where π0 = (1/6 , 1/6 , 1/6 , 1/6 , 1/6 , 1/6 ).

This is our assumed model, and under this H0, the expected counts are Ej = 30*1/6= 5 for each cell j. One way to calculate the difference between what we are assuming and what we are observing is to evaluate the the goodness-of-fit statistics as discussed in  the previous section.

\begin{eqnarray*}
X^2 &= & \dfrac{(3-5)^2}{5}+\dfrac{(7-5)^2}{5}+\dfrac{(5-5)^2}{5}\\
& & +\dfrac{(10-5)^2}{5}+\dfrac{(2-5)^2}{5}+\dfrac{(3-5)^2}{5}\\
&=& 9.2
\end{eqnarray*}

\begin{eqnarray*}
G^2 &=& 2\left(3\text{log}\dfrac{3}{5}+7\text{log}\dfrac{7}{5}+5\text{log}\dfrac{5}{5}\right.\\
& & \left.+ 10\text{log}\dfrac{10}{5}+2\text{log}\dfrac{2}{5}+3\text{log}\dfrac{3}{5}\right)\\
&=& 8.8
\end{eqnarray*}

Note that numerical values of Χ2 and G2 can be different. We can calculate the p-values for these statistics in order to help determine how well this model fits. The p-values are P25 ≥ 9.2) = .10 and P25 ≥ 8.8) = .12. Given these p-values, with the critical value or Type I error of α=0.05, we fail to reject the null hypothesis. But rather than relying on p-values, we should also think about the actual values of X2 and G2statistics. Small values imply a good fit here, i.e., the observed values are close to what you had assumed. Large values are going to imply the opposite. In this case, the fair-dice model doesn't fit the data very well, but the fit isn't bad enough to conclude beyond a reasonable doubt that the dice is unfair.

Next, we show how to do this in SAS and R; you can run either.

The following SAS code, dice_rolls.sas, will perform the goodness-of-fit test for the example above.

Here is the output generated by this program, and the headings are self-explanatory:

Notice that this SAS code only computes the Pearson chi-square statistic and not the deviance statistic.

Here is a short video describing the above output.
 Can you identify the relevant statistics and the p-value in the output? What does the column labeled "Percent" represent?

The following R code, dice_rolls.R will perform the same analysis as what we just did in SAS. You can copy the entire code in the R window and your output will be saved into two files, dice_rolls.out and dice_rolls_Results. If you are new to R, however, I suggest you run this line by line to get more familiar with the commands; you can simply copy and past into the R-terminal. For other options, see a very basic intro to R.

 dice_rolls.R dice_rolls.out (part of the output)

Here are a few additional comments regarding the above R code and its output. If you run Inspect, the data there will be labeled as "ex7" where in the updated code the label is "dice", but it's the same data. The last part of the R code does the same analysis, but gives an output that looks more like what SAS produces; we did this to show a bit more programming in R and how you can play with creating different outputs.

IMPORTANT! The Pearson chi-squared statistic is X2 = 9.2 (p-value=0.101), the deviance statistic (from the R output) is G2=8.78 (p-value=0.118), and they both follow the chi-squared sampling distribution with df=5. It is not uncommon to observe the discrepancy in the values between these two statistics especially for the smaller sample sizes. Notice, however, that in this case they do lead to the same conclusion --- there is a moderate evidence that the dice is fair.

 Can you identify the relevant statistics and the p-value in the output? What does the column labeled "Percentage" in dice_rolls.out represent?

### Example - Tomato Phenotypes

Tall cut-leaf tomatoes were crossed with dwarf potato-leaf tomatoes, and n = 1611 offspring were classified by their phenotypes.

 Phenotypes Count tall cut-leaf 926 tall potato-leaf 288 dwarf cut-leaf 293 dwarf potato-leaf 104

Genetic theory says that the four phenotypes should occur with relative frequencies 9 : 3 : 3 : 1, and thus are not all equally as likely to be observed. The dwarf potato-leaf is less likely to observed than the others. Do the observed data support this theory?

Under the null hypothesis, the probabilities are

π1 = 9/16 , π2 = π3 = 3/16 , π4 = 1/16 ,

and the expected frequencies are

E1 = 1611(9/16) = 906.2,
E
2 = E3 = 1611(3/16) = 302.1, and
E
4 = 1611(1/16) = 100.7.

We calculate the fit statistics and find that X2 = 1.47 and G2 = 1.48, which are nearly identical. The p-values based on the χ2 distribution with 3 d.f. are approximately equal to 0.69. Therefore, we fail to reject the null hypothesis and the data are consistent with the genetic theory. Again, the small values of these statistics imply that the fit between the data and the proposed model is good!

Here is how to do the computations in R using the following code :

This has step-by-step calculations and using of built-in chisq.test() with producing some nice output including Pearson and deviance residuals.

Here is the above analysis done in SAS:

 Do you recall what the residuals are from linear regression? How would you define them in this context? What do they tell you about the tomato example?