13.3 – Test assumption of normality

Introduction

I’ve commented numerous times that your conclusions from statistical inference are only as good as the validity of making the and applying the correct procedures. This implies that we know the assumptions that go into the various statistical tests, and where possible, we critically test the assumptions. From time to time then I will provide you with “tests of assumptions.”

Here’s one. The assumption of normality, that your data were sampled from a population with a normal distribution for the variable of interest, is key and there are a number of ways to test this assumption. Hopefully as you read through the next section you can extend the logic to any distribution; if the data are presumed to come from a binomial, or a Poisson, or a uniform distribution, then the logic of goodness of fit tests would apply.

How to test normality assumption

It’s not a statistical test per se, but the best is to simply plot (histogram) the data. You can do these graphics by hand, or, install the data mining package rattle which will generate nice plots useful for diagnosing normality issues.

Note 1: rattle (R Analytic Tool To Learn Easily) is a great “data-mining” package. We introduced the package in Chapter 3 — Exploring data.

The rattle histogram plot (Fig 1) superimposes a normal curve over the data, which allows you to “eyeball” the data.

First, the eye test. I used the R-package rattle for this on a data set of comet tail lengths of rat lung cells exposed to different amounts of copper in growth media (scroll to bottom of page or click here to get the data).

In addition to the histogram (Fig 1 top image), I plotted the cumulative function (Fig 1 bottom image). In short, if the data set were normal, then the cumulative frequency plot should look like a straight line.

rattle normal

cumulative frequency

Figure 1. Rattle descriptive graphics on Comet Copper dataset. Dotted line (top image) and red line (bottom image) follow the combined observations regardless of treatment.

So, just looking at the data set, we don’t see clear evidence for a normal-like data set. The top image (Fig 1) looks stacked to the left and the cumulative plot (bottom image) is bumped in the middle, not falling on a straight line. We’ll need to investigate the assumption of normality more for this data set. We’ll begin by discussing some hypothetical points first, then return to the data set.

Goodness of fit tests of normality hypothesis

While graphics and the “eyeball test” are very helpful, you should understand that whether or not your data fits a normal distribution, that’s a testable hypothesis. The null hypothesis is that your sample distribution fits are normal distribution. In general terms, these “fit” hypotheses are viewed as “goodness of fit” tests. Often times, the test is some variation of a \chi^{2} problem: you have your observed (sample distribution) and you compare it to some theoretical expected value (eg, in this situation, the normal distribution). If your data fit a normal curve, then the test statistics will be close to zero.

We have discussed before that the data should be from a normal distributed population. To the extent this is true, then we can trust our statistical tests are performing the way they are expected to do. We can test the assumption of normality by comparing our sample distribution against a theoretical distribution, the normal curve. I’ve shown several graphs in the past that “looked normal”. What are the alternatives for unimodal (single peak) distributions (Fig 2)?

Kurtosis describes the shape of the distribution, whether it is stacked up in the middle (leptokurtosis), or more spread out and flattened (platykurtosis).

Skewness describes differences from symmetry about the middle. Left skew the tail of the distribution extends to the left, ie, smaller values.

“Leptokurtosis”

Leptokurtosis

“Platykurtosis”

Platykurtosis

Negative skew, left skewed

Negative skew, left skewed

Positive skew, right skewed

Positive skew, right skewed

Figure 2. Graphs describing different distributions. From top to bottom: Leptokurtosis, platykurtosis, negative skew, positive skew.

The easiest procedures for goodness of fit tests of normality are based on the \chi^{2} distribution and yield a “goodness of fit” test for normal distribution. We discussed the \chi^{2} distribution in Chapter 6.9, and used the test in Chapter 9.1.

    \begin{align*} \chi^2 = \sum_{i = 1}^{k}\frac{\left ( O_{i}-E{i} \right )^2}{E_{i}} \end{align*}

where O refers to the observed data (what we’ve got) and E refers to the expected (eg, data from a normal curve with same mean and standard deviation as our sample).

To illustrate, I simulated a data set in R.

Rcmdr: Distributions → Continuous distributions → Normal distribution → Sample from normal distribution

I created 100 values, stored them in column 1, and set the population mean = 125 and population standard deviation = 10. And therefore the population standard error of the mean was 1.0.

The resulting MIN/MAX was from 99.558 to 146.16; the sample mean was 124.59 with a sample standard deviation of 9.9164. And therefore the sample standard error of the mean was 0.9916.

Question: After completing the steps above, your data will be slightly different from mine… Why?

But getting back to the main concept, Does our data agree with a normal curve? We have discussed how to construct histograms and frequency distributions.

Let’s try six categories (why six? we discussed this when we talked about histograms).

All chi-square tests are based on categorical data, so we use the counts per category to get our data. Group the data, then count the number of OBSERVED in each group. To get the EXPECTED values, use the Z-score (normal deviate) with population mean and standard deviation as entered above.

Table 1. Tabulated values for test of normality.

Number of observations Weight Normal deviate (Z) Expected Proportion Expected number (Obs – Exp)2 / Exp
105 3 less than or equal to -2 0.0228 2.28 0.227368421
105 < 115 17 between -1 & -2 = 0.1587 – 0.0228 0.1359 13.59 0.855636497
115 < 125 34 between -0 & -1 = 0.5 – 0.1587 0.3413 34.13 0.000495166
125 < 135 30 between +0 & +1 = 0.5 – 0.1587 0.3413 34.13 0.499762672
135 < 145 15 between +1 & +2 = 0.9772 – 0.8413 0.1359 13.59 0.146291391
> 145 1 greater than or equal to +2 = 1 – 0.9772 0.0228 2.28 0.718596491
\chi^2 = 2.44815064

Then obtain the critical value of the \chi^{2} with df = 6 – 1 = 5 (see Appendix, Table of chi-square critical values , critical \chi^{2} = 11.1 with df = 5).

Thus, we would not reject the null hypothesis and would proceed with the assumption that our data could have come from a normally distributed population.

This would be an OK test, but different approaches, although based on a chi-square-like goodness of fit, have been introduced and are generally preferred. We have just shown how one could use the chi-square goodness of fit approach to testing whether your data fit a normal distribution. A number of modified tests based on this procedure have been proposed; we have already introduced and used one (Wilks-Shapiro), which is easily accessed in R.

Rcmdr: Summaries → Wilks-Shapiro

Like Wilks-Shapiro, another common “goodness-of-fit” test of normality is the Anderson-Darling test. This test is now included with Rcmdr, but it’s also available in the package nortest. After the package is installed and you run the library (you should by now be able to do this!), then at the R prompt (>) type:

ad.test(dataset$variable)

replacing “dataset$variable” with the name of your data set and variable name.

In the context of goodness of fit, a perfect fit means the data are exactly distributed as a normal curve, and the test statistics would be zero. Differences away from normality increase the value of the test statistic.

How do these tests perform on data?

The histogram of our simulated normal data of 100 observations with mean = 125 and standard deviation = 10 (Fig 3).

Histogram of simulated normal dataset, μ = 125, σ = 10.

Figure 3. Histogram of simulated normal dataset, μ = 125, σ = 10.

and here’s the cumulative plot (Fig 4).

Cumulated frequency of simulated normal dataset, μ = 125, σ = 10

Figure 4. Cumulated frequency of simulated normal dataset, μ = 125, σ = 10.

Results from Anderson-Darling test were A = 0.2491, p-value = 0.7412, where A is the Anderson-Darling test statistic.

Results of the Shapiro-Wilks test on the same data. W = 0.9927, p-value = 0.8716, where W is the Shapiro-Wilks test statistic. We would not reject the null hypothesis in either case because p > 0.05.

Example

Histogram of a data set, highly skewed to the right. 90 observations in three groups of 30 each, with mean = 0 and standard deviation = 1 (Fig 5).

Histogram of simulated normal dataset, μ = 0, σ = 1.

Figure 5. Histogram of simulated normal dataset, μ = 0, σ = 1.

and the cumulative frequency plot (Fig 6).

Cumulated frequency of simulated normal dataset, μ = 0, σ = 1.

Figure 6. Cumulated frequency of simulated normal dataset, μ = 0, σ = 1.

Results from Anderson-Darling test were A = 9.0662, p-value < 2.2e-16. Results of the Shapiro-Wilks test on the same data. W = 0.6192, p-value = 6.248e-14. Therefore, we would reject the null hypothesis because p < 0.05.

The Shapiro-Wilk test in Rcmdr

Let’s go back to our data set and try tests of normality on the entire data set, ie, not by treatment groups.

Rcmdr: Statistics → Summaries → Test of normality…

normalityTest(~CometTail, test="shapiro.test", data=CometCopper)

 Shapiro-Wilk normality test

data: CometCopper$CometTail
W = 0.91662, p-value = 0.006038

End of R output

Another test, built on the basic idea of a chi-square goodness of fit is the Anderson Darling test. Some statisticians prefer this test and it is one built-in to some commercial statistical packages (eg, Minitab). To obtain the Anderson-Darling test in R, you need to install a package. After installing nortest package, run the AD test at the command prompt.

require(nortest)
ad.test(CometCopper$CometTail)

Anderson-Darling normality test

data: CometCopper$CometTail
A = 1.0833, p-value = 0.006787

End of R output.

Note 2: Because you’ve installed Rcmdr the Anderson-Darling test is available (since version 2.6) without installing the nortest package. You can call Anderson-Darling via the menu system, Statistics → Summaries → Test of normality … or type in the script the following command

normalityTest(~CometTail, test="ad.test", data=CometCopper)

In contrast, Shapiro-Wilks test is part of the base R package, so can be called directly

shapiro.test(CometCopper$CometTail)

P-values are both much less than 0.05, so we would reject the assumption of normality for this data set.

Which test of normality?

Why show you two tests for normality, the Shapiro-Wilks and Anderson-Darling? The simple answer is that both are good as general tests of normality, both are widely used in scientific papers, so just pick one and go with it as your general test of normality.

The more complicated answer is that each is designed to be sensitive to different kinds of departure from normality. By some measures, the Shapiro-Wilks test is somewhat better (ie, more statistical power to test the null hypothesis), than other tests, but this is not something you want to get into as a beginner. So, I show both of them to you so that you are at least introduced to the concept that there are often more than one way to test a hypothesis. The bottom-line is that plots may be best!

Questions

  1. Work describe in this chapter involves statistical tests of the assumption of normality. It is just as important, maybe more so, to also apply graphics to take advantage of our built-in pattern recognition functions. What graphic techniques, besides histogram, should be used to view the distribution of the data?
  2. In R, what command would you use so that you can call the variable name, CometTail, directly instead of having to refer to the variable as CometCopper$CometTail?
  3. Why are Anderson-Darling, Shapiro-Wilks and other related tests referred to as “goodness of fit” tests? You may wish to review discussion in Chapter 9.1.
  4. The example tests presented for the Comet Copper data set were conducted on the whole set, not by treatment groups. Re-run tests of normality via Rcmdr, but this time, select the By groups option and select Treatment.

Quiz Chapter 13.3

Test assumption of normality

Data set used in this page

Treatment CometTail
Control 17.86
Control 16.52
Control 14.93
Control 14.03
Control 13.33
Control 8.81
Control 14.70
Control 9.26
Control 21.78
Control 6.18
Control 9.20
Control 5.54
Control 6.72
Control 2.63
Control 7.19
Control 5.39
Control 11.29
Control 15.44
Control 17.86
Contro l4.25
Copper 53.21
Copper 38.93
Copper 18.93
Copper 30.00
Copper 28.93
Copper 15.36
Copper 17.86
Copper 17.50
Copper 21.07
Copper 29.29
Copper 28.21
Copper 16.79
Copper 21.07
Copper 37.50
Copper 38.22
Copper 17.86
Copper 29.64
Copper 11.07
Copper 35.00
Copper 49.29

Chapter 13 contents

10.1 – Compare two independent sample means

Introduction

We introduced the concept of comparing a sample statistic (mean) against a population parameter (Chapter 6.7, Normal deviate) or one-sample t-test against a specified mean (eg, from published data or from theory, Chapter 8.5).

Consider now a basic experimental design, the randomized control trial, or RCT (Fig 1), introduced in Chapter 2.4.

simple RCT

Figure 1. A two group Randomized Control Trial.

Subjects randomly selected from population of interest, then again — random assignment — once recruited into one of two treatment groups. Importantly, subjects belong to one treatment arm only: no subject simultaneously receives the treatment and the control. This is in contrast to the paired design, in which subjects receive both treatments (see Chapter 10.3).

In inferential statistics about an experiment, we are more likely trying to test if sample means are different. For example

  • two species grown in common garden, do they differ in growth rate?
  • human subjects given a new treatment have better outcomes compare to those receiving a control treatment (eg, placebo).

The equivalent null hypothesis is that two samples are pulled from the same population. We write the null hypothesis as

    \begin{align*} H_{0}\ :\ \bar{X}_{1}=\bar{X}_{2} \end{align*}

and the corresponding alternative hypothesis, HA, then must be

    \begin{align*} H_{A}:\bar{X}_{1}\neq \bar{X}_{2} \end{align*}

Question: Is this a one tailed or two-tailed hypothesis?

Answer: Two-tailed (review Chapter 8.4)

Note 1: In this day and age, there’s really no compelling reason to learn the t-test. First, it is just a special case of the one-way ANOVA, therefore, it’s a special case of the general linear model. Struggling to learn R commands? Well, one solution would be just to learn the general linear model approach — the R function lm() (OK — don’t get too excited — lm() has many options and details). Second, few experiments or observational studies are likely to have only two groups; thus, the temptation to carry out a series of t-tests, taking all groups two at a time, or “pairwise,” while tempting, actually violates a whole bunch of basic statistical rules (discussed in Chapter 12.1). It will also make statisticians go crazy when they see it. That said, if your experiment has but two groups, then by all means, the t-test is a choice. The t-test is also a statistical test that you have likely already used before so we present the discussion here to build on what you may already have learned. We also present the independent t-test as a vehicle.

Worked example

We introduce the two-sample t-test, or better, the independent sample t-test.

    \begin{align*} t=\frac{\bar{X}_{1}-\bar{X}_{2}}{s_{\bar{X}_{1}-\bar{X}_{2}}} \end{align*}

where the numerator is the difference between the two sample means and the denominator is the standard error of the differences between the two groups standard errors. The formula for

    \begin{align*} s_{\bar{X}_{1}-\bar{X}_{2}}=\sqrt{ \frac{s_1^2}{n_{1}}+\frac{s_2^2}{n_{2}}} \end{align*}

The choice of independent sample over two-sample is best because it emphasizes that the two groups (the two samples), must be comprise of independent sampling units. This is a pretty straight-forward requirement; you have randomly assigned twenty individuals to two groups, a control group (n = 10) and a treatment group (n = 10). Individuals are either in the control group or they are in the treatment group — they cannot simultaneously appear in both groups.

We will work our way through this test by example. For starters, let’s use the same lizard dataset (see Example data set, below), four body mass recordings (grams) each for house geckos (Hemidactylus frenatus, Fig 2) and the Carolina anole (Anolis carolensis, Fig 3), two of many lizard species introduced to Hawaiʻi.

Male Hemidactylus frenatus, central Oahu, M. Dohm 2018

Figure 2. Male Hemidactylus frenatus, central Oahu, M. Dohm.

Male Anolis carolinensis, ʻAkaka Falls, Big Island of Hawai`i, M. Dohm 2018

Figure 3. Male Anolis carolinensis, ʻAkaka Falls, Big Island of Hawaiʻi, M. Dohm.

Example data set

Geckos: 3.186, 2.427, 4.031, 1.995

Anoles: 5.515, 5.659, 6.739, 3.184

Question: How would you go about creating a data frame with the values in long form (stacked worksheet), including a label variable and the body mass?

Note 2: The independent sample t-test in Rcmdr requires a stacked worksheet and not in unstacked worksheet two columns implied by the above pseudo-code. If you need help with worksheet format, then see Part07 in Mike’s Workbook for Biostatistics.

Answer: At the R prompt, type

Geckos <- c(3.186, 2.427, 4.031, 1.995); Anoles = c(5.515, 5.659, 6.739, 3.184) #create two vectors
bmass <- c(Geckos, Anoles)    #combine the two vectors into a single vector holding all of the body mass records
species <- c("gecko", "gecko", "gecko", "gecko", "anole", "anole", "anole", "anole")   #create your label variable
lizards <- data.frame(species, bmass)  #create your data frame
lizards   #print your data frame  species bmass
1   gecko 3.186
2   gecko 2.427
3   gecko 4.031
4   gecko 1.995
5   anole 5.515
6   anole 5.659
7   anole 6.739
8   anole 3.184

Note also that you can enter data into the Data editor by creating the data frame first then adding values. To edit the data frame “lizards” type fix(lizards) at the R prompt, then close the data frame when you have added or changed values as needed.

As always, begin with an exploration of the data, including a graph (Fig 4).

Box plot of lizard body mass.

Figure 4. Box plot of lizard body mass.

We can see already that there’s greater spread of data for the Anoles compared to the Geckos, but the median values differ. Small sample sizes can be a problem for analyses as we can only have reduced confidence in our conclusions. However, we press on for the sake of demonstration.

Let’s test the null hypothesis, H0, ie, the two species of lizards have the same mean body mass.

Rcmdr: Statistics → Means → Independent-samples t-test…

In this next image I posted the Rcmdr menu popup for the Independent Samples t-test (Fig 5). Later versions of Rcmdr split the settings for this command into two tabs; the first tab allows for the selection of the variables and setting the hypotheses whereas the second tab, labeled Options, permits additional choices. The default selections need your attention: to actually conduct the t-test you need to answer “No” to the question, “Assume equal variance?” (Fig 6).

Screenshot Rcmdr t-test menu

Figure 5. Rcmdr menu for Independent sample t-test.

Select the Options tab (Fig 6) to select null hypothesis and to select the t-test and not the Welch-test (which is the default, ie, No to the prompt “Assume equal variances?”).

Screenshot Rcmdr t-test options

Figure 6. Rcmdr Options menu for Independent sample t-test.

Let’s look at the results and break down the parts of the test.

t.test(Body.mass~Lizard, alternative='two.sided', conf.level=.95, 
+ var.equal=TRUE, data=lizards)

Two Sample t-test

data: Body.mass by Lizard 
t = 2.7117, df = 6, p-value = 0.03503
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
0.2308685 4.4981315 
sample estimates:
mean in group Anolis mean in group Gecko 
5.27425 2.90975

Consider the R session output above and answer the following questions.

Questions for the worked example

  1. Which lizard group had the greater mean value, Anolis or Gecko?
  2. What are the assumptions necessary for you to use the independent sample t-test?
  3. What does “two-sided” mean?
  4. What was the null hypothesis?
  5. Was this a one-tailed or two-tailed test of the null hypothesis?
  6. What is the value of the test statistic?
  7. How many degrees of freedom?
  8. What is the critical value for this test?
  9. What is the value of the lower limit of the 95% confidence interval?
  10. What is the value of the lower limit of the 99% confidence interval?
  11. True or False. If the null hypothesis is accepted, then zero is a value included in the 95% confidence interval.
  12. Do you accept the null hypothesis? Explain your selection.

Try another example

DNA damage, changes in the chemical structure of nucleotide bases or breakage of the DNA chains, occurs in cells under many circumstances. The comet assay, or single-cell gel electrophoresis, is one method for visualizing and measuring DNA strand breaks in cells. Exposed cells are mixed with a low-melting temperature agarose and placed onto a microscope slide. The cells are then lysed with an alkaline detergent and high salts. When current is applied across the slide, undamaged DNA remains in the nucleus, whereas damaged DNA extends towards the anode to form a comet-like tail, with imaging assisted by including a fluorescent dye like Sybr-Green. Examples of comets are shown below (Fig 7).

example comet assay

Figure 7. Comet examples. A Intact cell, no DNA damage, B Cell with some DNA damage, a slight tail to the right is evident, C  Cell with significant DNA damage, a large tail is evident. M. Dohm.

In an experiment, immortalized lung epithelial cells were exposed to dilute copper solutions for 30 minutes then washed with PBS. The comet assay was applied to these cells and for comparison, to cells without copper exposure but otherwise treated the same way (controls). The data are available at the bottom of this page (scroll down or click here).

Again, you should begin all analyses with an exploration of the data, including a graph (Fig 8).

box plot comet tails

Figure 8. Boxplot of comet tail lengths for cells with and without (control) exposure to copper in the cell medium for 30 minutes.

Let’s look at the R output for the t-test analysis.

    Two Sample t-test

data:  CometTail by Treatment
t = -5.8502, df = 38, p-value = 9.139e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -22.39865 -10.88213
sample estimates:
mean in group Control  mean in group Copper 
             11.14533              27.78571

Consider the R session output above and answer the following questions.

Questions for Comet assay data set

  1. Which cell group had the greater mean value, Copper-exposed or Control-exposed cells?
  2. What are the assumptions necessary for you to use the independent sample t-test?
  3. What does “two-sided” mean?
  4. What was the null hypothesis?
  5. Was this a one-tailed or two-tailed test of the null hypothesis?
  6. What is the value of the test statistic?
  7. How many degrees of freedom?
  8. What is the critical value for this test?
  9. What is the value of the lower limit of the 95% confidence interval?
  10. What is the value of the lower limit of the 99% confidence interval?
  11. True or False. If the null hypothesis is accepted, then zero is a value included in the 95% confidence interval.
  12. Do you accept the null hypothesis? Explain your selection.

T test from summary statistics

In some cases you may only have to summary statistics for data, eg, the means and the standard deviations. We can use the equations of the t test to write a simple formula, where the user provides the known means, standard deviations, and sample size. For example, create a simple function with readline for user input.

myTtest <- function() {
   mnx <- as.numeric(readline(prompt="Enter mean of x: "))
   stdevx <- as.numeric(readline(prompt="Enter sd of x: "))
   nx <- as.numeric(readline(prompt="Enter n of x: "))
   mny <- as.numeric(readline(prompt="Enter mean of y: "))
   stdevy <- as.numeric(readline("Enter sd of y: "))
   ny <- as.numeric(readline(prompt="Enter n of y: "))
   myTvalue <- abs(((mnx-mny)-0)/sqrt(((stdevx^2)/nx)+(stdevy^2)/ny))
   myDF <- as.integer(nx+ny-2)
   myPvalue <- pt(myTvalue,myDF,lower.tail=FALSE)*2
   myResults <- c(myTvalue, myDF, myPvalue)
   report <- c("T-test: ", "df: ", "two-tailed p-value: ")
   cat(sprintf("%s %3.3f, ", report, myResults))
}

then run the function by typing myTest() at the R prompt and entering the means, standard deviations, and sample size when prompted.

myTtest()
Enter mean of x: 2.91
Enter sd of x: .895
Enter n of x: 4
Enter mean of y: 5.27
Enter sd of y: 1.497
Enter n of y: 4
T-test: 2.706, df: 6.000, two-tailed p-value: 0.035

P-values from confidence intervals

While we expect certain reporting criteria for published, it is not uncommon that one or more elements are missing. For example, while means plus or minus standard deviations for two groups may be reported, the confidence interval of the difference may not be reported, or, an inexact p-value is reported like “< 0.05,” but we “need” an exact p-value for our meta-analysis.  A little math, and these missing statistics can be calculated.

We can use the lizard study as an example; while all of the expected elements were reported, including the p-value (0.03503), let’s say the p value was reported as “statistically significant at p < 0.05.”

We need

lower and upper confidence intervals: 0.2309 and 4.498, respectively.

critical value from t distribution, df = 6 (2 tailed, therefore 0.05/2). A little R code we have:

critVal <- qt(c(.05), df=6, lower.tail=FALSE)=2.447

the standard error, which can be calculated as SEliz = (upper − lower)/(2×critVal)=1.097968

and then the p-value, exp(-.717*testLiz-0.416*testLiz^2) = 0.03101

which is pretty close to the result from R (p = 0.03503). The difference is likely due to the small sample size.

Questions

  1. Don’t forget to work through the Questions for the Comet tail data set (scroll up or click here).
  2. Microsoft Excel, LibreOffice Calc, and Google sheets spreadsheet software all include t-test functions and return the p-value. Consider two variables big (100, 110, 120, 100, 110, 210, 200) and small (0,1,1,2,0,1,0). (Note — these two groups are obviously very different, calculating a t-test on their difference is silly, just for this question.)  If formatting is set to the default two decimal places for Number cell category, the p-value will return as “0.00.” How should you report the p-value in this case?

Quiz Chapter 10.1

Compare two independent sample means

Comet assay data set

Treatment CometTail
Control 17.856139
Control 16.52125
Control 14.925449
Control 14.029174
Control 13.332945
Control 8.811185
Control 14.701654
Control 9.261025
Control 21.779311
Control 6.180284
Control 9.201752
Control 5.54472
Control 6.717885
Control 2.625092
Control 7.191583
Control 5.392866
Control 11.284813
Control 15.441254
Control 17.857176
Control 4.250956
Copper 53.214287
Copper 38.92857
Copper 18.928572
Copper 30
Copper 28.928572
Copper 15.357142
Copper 17.857143
Copper 17.5
Copper 21.071428
Copper 29.285715
Copper 28.214285
Copper 16.785715
Copper 21.071428
Copper 37.5
Copper 38.214287
Copper 17.857143
Copper 29.642857
Copper 11.071428
Copper 35
Copper 49.285713

Chapter 10 contents

6.9 – Chi-square distribution

Introduction

As noted earlier, the normal deviate or Z score can be viewed as randomly sampled from the standard normal distribution. The chi-square distribution describes the probability distribution of the squared standardized normal deviates with degrees of freedom, df, equal to the number of samples taken. (The number of independent pieces of information needed to calculate the estimate, see Ch. 8.) We will use the chi-square distribution to test statistical significance of categorical variables in goodness of fit tests and contingency table problems.
The equation of the chi-square is

    \begin{align*} \chi^2 = \sum_{i = 1}^{k}\frac{\left ( f_{i}-\hat{f}_{i} \right )^2}{\hat{f}_{i}} \end{align*}

where k is the number of groups or categories, from 1 to k, and fi is the observed frequency and fi “hat” is the expected frequency for the kth category. We call the result of this calculation the chi-square test statistic. We evaluate how often that value or greater of a test statistic will occur by applying the chi-square distribution function. Graphs to show chi-square distribution for degrees of freedom equal to 1 – 5, 10, and 20 (Fig 1).

GIF, chi-square distribution, df = 1 - 5, 10, and 20

Figure 1. Animated GIF of plots of chi-square distribution over range of degrees of freedom.

Note that the distribution is asymmetric, particularly at low degrees of freedom. Thus tests using the chi-square are one-tailed (Fig 2).

upper-tail chi-square distribution

Figure 2. The test of the chi-square is typically one-tailed. In this case, probability of values greater than the critical value.

By convention in the Null Hypothesis Significance Testing protocol (NHST), we compare the test statistic to a critical value. The critical value is defined as the value of the test statistic — the cutoff boundary between statistical significance and insignificance — that occurs at the Type I error rate, which is typically set to 5%. The interpretation of the result is as follows: after calculating a test statistic, we can judge significance of the results relative to the null hypothesis expectation. If our test statistics is greater than the critical value, then the p-value of our results are less than 5% (R will report an exact p-value for the test statistic).  You are not expected to be able to follow this logic just yet — rather, we teach it now as a sort of mechanical understanding to develop in the NHST tradition. The justification for this approach to testing of statistical significance is developed in Chapter 8. A portion of the critical values of the chi-square distribution are shown in Figure 3.

Portion of chi-square table

Figure 3. Portion of the table of some critical values of chi-square distribution, one tailed (right-tailed or “upper” portion of distribution).

See Appendix for a complete chi-square table.

The chi-squared test is a one-tailed test.

The result of the calculation of the test statistic is always non-negative because it is calculated by squaring the differences between the observed and expected frequencies.  As a consequence, the test is always “right-tailed.”

The chi-square test does not provide information about the direction of an association, nor the strength of the association, only whether or not the variables are independent or associated. We use “correlation,” (eg, Phi coefficient or Cramer’s V, see Chapter 9.1 – Chi-square test: Goodness of fit and Chapter 9.2 – Chi-square contingency tables), to quantify the strength and direction of the association. While statistical significance is important, examining the direction and magnitude, the effect size, provides actionable insight and context for making informed decisions. For example, a new drug tested in a clinical trial may show an association between subject improvement (statistical significance), but yet have a small effect.

Example

Professor Hermon Bumpus of Brown University in Providence, Rhode Island, received 136 House Sparrows (Passer domesticus) after a severe winter storm 1 February 1898. The birds were collected from the ground; 72 of the birds survived, 64 did not (Table 1). Bumpus made several measures of morphology on the birds and the data set has served as a classical example of Natural Selection (Chicago Field Museum). We’ll look at this data set when we introduce Linear Regression.

Table 1. Survival statistics of Bumpus House sparrows

Yes No
Female 21 28
Male 51 36

Was there a survival difference between male and female House Sparrows? This is a classic contingency table analysis, something we will at length in Chapter 9. For now, we report the Chi-square test statistic for this test was 3.1264 and the test had one degree of freedom. What is the critical value of the chi-square distribution at 5% and one degree of freedom?. Typically we would simply use R to look this up

qchisq(c(0.05), df=1, lower.tail=FALSE)

But we can also get this from the table of critical values (Fig 4). Simply select the row based on the degrees of freedom for the test then scan to the column with the appropriate significance level, again, typically 5% (0.05).

Figure 4. Portion of the chi-square distribution which shows how to find critical value of the chi-square distribution. 

Figure 4. Portion of the chi-square distribution which shows how to find critical value of the chi-square distribution. 

For 1 degree of freedom at 5% significance, the critical value is 3.841. Back to our hypothesis: Did male and female survival differ in the Bumpus data set? Following the NHST logic, if the test statistic value (eg, 3.1264) is greater than the critical value (3.841), then we would reject the null hypothesis. For this example, we would conclude no statistical difference between male and female survival because the test statistic was smaller than the critical value. How likely are these results due to chance? That’s where the p-value comes in. Our test statistic value falls between 5% and 10% (2.706 < 3.1264 < 3.841). In order to get the actual p-value of our test statistic we would need to use R.

R code

Given a chi-square test statistic you can use R to calculate the probability of that value against the null hypothesis. At the R prompt

pchisq(c(3.1264), df=1, lower.tail=FALSE)

And R output

[1] 0.07703368

Because we are using R Commander, simply select the command by following the menu options.

Rcmdr: Distributions → Continuous distributions → Chi-squared distribution → Chi-squared probabilities …

Enter the chi-square value and degrees of freedom (Fig 5).

Figure 4. Screenshot of input box in Rcmdr for Chi-square probability values.

Figure 5. Screenshot of input box in Rcmdr for Chi-square probability values.

Questions

  1. What happens to the shape of the chi-square distribution as degrees of freedom are increased from 1 to 5 to 20 to 100?

Be able to answer these questions using the Chi-square table, Appendix 20.2, or using Rcmdr

  1. For probability α = 5%, what is the critical value of the chi-square distribution (upper tail)?
  2. The value of the chi-square test statistic is given as 12. With 3 degrees of freedom, what is the approximate probability of this value, or greater from the chi-square distribution?

Quiz Chapter 6.9

Chi-square distribution


Chapter 6 contents

6.7 – Normal distribution and the normal deviate (Z)

Introduction

In chapter 3.3 we introduced the normal distribution and the Z score, aka normal deviate, as part of a discussion about how some knowledge about characteristics of the dispersion of our data sampled from a population could be used to calculate how many samples we need (the empirical rule). We introduced Chebyshev’s inequality as a general approach to this problem, where little is known about the distribution of the population and contrasted it with the Z score, for cases where the distribution is known to be Gaussian or the normal distribution. The normal distribution is one of the most important distributions in classical statistics. All normal distributions are bell-shaped and symmetric about the mean. To describe a normal distribution only two parameters are needed: the population mean, μ, and the population standard deviation, σ. The normal distribution with mean equal to zero and standard deviation equal to one is called the standard normal, or Z distribution. With use of the Z score any normal distribution can be quickly converted to the standard normal distribution.

Proportions of a Normal Distribution

This concept will become increasingly important for the many statistical tests we will learn over the next few weeks. What is the proportion of the populations that is greater than some specific value? Below (Fig 1), again, I have generated a large data set, now with population mean µ = 5 and σ = 2. The red line corresponds to the equation of the normal curve using our values of µ = 5 and σ = 2.

Frequency of observations expected to be greater than 7 (red bars) from a large population with mean µ = 5 and σ = 2.

Figure 1. Frequency of observations expected to be greater than 7 (red bars) from a large population with mean µ = 5 and σ = 2.

R code for Fig 1.

set.seed(123)
myData <- rnorm(1000, mean = 5, sd = 2)
h <- hist(myData, breaks = 20, plot = FALSE) 
threshold <- 7
# Bars with midpoints greater than threshold will be "red", 
# and others will be "white."
bar_colors <- ifelse(h$mids > threshold, "red", "white")
plot(h, col = bar_colors, xlab = "Value", ylab = "Frequency")

Note that this is a crucial step! We assume that our sample distribution is really a sample from a population density (= “area under the curve”) function (= “an equation”) for a normal random (= “population”) variable.

Once I (you) make this assumption, then we have powerful and easy to use tools at our command to answer questions like

Question: What proportion of the population is greater than 7 (Fig 1, colored in red)?

This gets to the heart of the often-asked question, How many samples should I measure? If we know something about the mean and the variability, then we can predict how many samples will be of a particular kind. Let’s solve the problem.

The Z score

We could use the formula for the normal curve (and a lot of repetitions), but fortunately, some folks have provided tables that short-cut this procedure. R and other programs also can find these numbers because the formulas are “built in” to the base packages. First, let’s introduce a simple formula that lets us standardize our population numbers so that we can use established tables of probabilities for the normal distribution.

Below, we will see how to use Rcmdr for these kinds of problems.

However, it’s one of the basic tasks in statistics that you should be able to do by hand. We’ll use the Z score as a way to take advantage of known properties of the standard normal curve.

    \begin{align*} Z = \frac{\left ( X_{i} - \mu\right )}{\sigma} \end{align*}

Z = 1 (with the mean and SD). Z (say “Z-score”) is called the normal deviate (aka “standard normal score”; it is also called the “Z-score”); it gives us a shortcut for finding the proportion of data greater than 7 in this case).

We use the normal deviate to do a couple of things; one use is to standardize a sample of observations so that they have a mean of zero and a standard deviation of one (the Z distribution). The data would then said to have been normalized.

The second use is to make predictions about how often a particular observation is likely to be encountered. As you can imagine, this last use is very helpful for designing an experiment — if we need to see a specified difference, we can conduct a pilot study (or refer to the literature) to determine a mean and level of variability for our observation of choice, plug these back into the normal equation and predict how likely we can expect to see a particular difference. In other words, this is one way to answer that question — how many observations need I make for my experiment to be valid?

Table of normal distribution

A portion of the table of the normal curve is provided at our web site and in your workbook. For our discussions, here’s another copy to look at (Fig 2).

portion of the table of normal distribution

Figure 2. Portion of the table of the normal distribution. Only values equal to or greater than Z = 0 are visible.

See Table 1 in the Appendix for a full version of the normal table.

We read values of Z from the first column and the first row. For Z = 0.23 we would scan the top row, scoot over to fourth column, then trace to where the row and column intersect (Fig 3); the frequency of occurrence of values at Z = 0.23 is 0.409046 or 40.9% (Fig 3).

Portion of normal table, z = 0.23 highlighted in red

Figure 3. Highlight Z = 0.23, frequency is 0.409046.

Z on the standard normal table is going to range between -4 and +4, with Z = 0 corresponding to 0.500. The Normal table values are symmetrical about the mean of zero.

What to make of the values of Z, from -4, -3, … +2, +3, up to +4 and beyond?  These are the standard deviations! Recall that using the Z score you corrected to mean of zero (got it!), and a standard deviation of one! A Z of 2 is twice the standard deviation; a Z = 3 is therefore three times the standard deviation, and so forth. The distribution is symmetrical: you get the same frequency for negative as for positive values. So on the “X” axis on a standard normal distribution, we have units of standard deviation plus (greater) or minus (less) than the mean. In Figure 4, area less than -1 standard deviations is highlighted (Fig 4).

Plot of standard normal distribution; area less than -1 σ.

Figure 4. Plot of standard normal distribution; area less than -1 σ.

R code for Fig 4.

set.seed(123)
myData2 <- rnorm(100000, mean = 0, sd = 1)
threshold <- -1
bar_colors <- ifelse(h$mids < threshold, "red", "white")
# Additional code used, same as Figure 1.

Question. How many multiples of standard deviations would you have for a Z score of Z = 1.75?

Answer = 1.75 times

Examples

See Table 1 in the Appendix for a full version of the normal table as you read this section.

What proportion of the data set will have values greater than 7? After applying our Z score equation, I get Z = 1.0 and for a Z of 1.0 = 0.1587 or 15.87% of the observations are greater than 7.

What proportion of the data set will have values less than -7? After applying our Z score equation, I get Z = -1.0. Taking advantage of symmetry argument, I just take my Z = -1.0 and make it positive — instead of values smaller than -1.0 we now have values greater than +1.0. And for a Z of 1.0 = 0.1587 or 15.87% of the observations are greater than 7, which means that 15.87% will be -7 or smaller.

What proportion of the data set will have values greater than 8? Again, apply the Z score equation I get Z of 1.5 = 0.0668 or 6.68% of the observations are greater than Z.

What proportion of the population is between 5 and 7 ? (colored in red). Draw the problem (example shown in Fig 5).

proportion of the population is between 5 and 7.

Figure 5. Proportion of the population is between 5 and 7.

R code for Fig 5.

# Used myData from Fig 1 code.

highlight_start <- 5
highlight_end <- 7
bar_colors <- rep("white", length(h$counts)) 
bar_colors[h$breaks[-length(h$breaks)] >= highlight_start & h$breaks[-1] <= highlight_end] <- "red"
# Additional code used, same as Figure 1.

Note 1: h$breaks[-length(h$breaks)] gives the left edges of the bins, and h$breaks[-1] gives the right edges.

Worked problem

1 – (proportion beyond 7) – (proportion less than 5)

1 – (0.1587) – (proportion less than 5)

And the proportion less than 5?

Use Z-score equation again. Now we find that

Z = 0

And look up Z of 0 from the table  = 0.5 proportion or 50.0%

Therefore, the proportion between 5 & 7 equals

1 -0.1587 – 0.50 = 0.3413

Answer = 34.13 % of the observations are between 5 and 7 when the mean µ = 5 and s = 2.

Questions

Several questions were asked in the text.

1. Proportion of observations greater than 7  from a large population with mean µ = 5 and σ = 2.

2. How many multiples of standard deviations would you have for a Z score of Z = 1.75?

3.

Quiz Chapter 6.7

Normal distribution questions


Chapter 6 contents