More tests

\(\chi^2\) Goodness of Fit Test

As you may remember, in the survey we asked students to complete in welcome week, one of the questions concerned the month in which you were born. You can download the data from https://uoepsy.github.io/data/surveydata_allcourse.csv.

We need to make sure we know how many observations we have. While there are currently 393 observations in our sample, we may have missing data in our variable of interest!

We can use filter to remove any which are missing (or NAs).
Take a bit of time to look at what each bit of this does. To help figure it out, try seeing what is.na(c(NA,1,2)) does. Then try !is.na(c(NA,1,2)).

survey_data <- read_csv("https://uoepsy.github.io/data/surveydata_allcourse.csv") 
survey_data <-
  survey_data %>% filter(!is.na(birthmonth))

Question A1

What is your intuition about the distribution of all students’ birth-months?
Do you think they will be spread uniformly across all months of the year (like a fair 12-sided dice), or do you think people are more likely to be born in certain months more than others?

Plot the distribution and get an initial idea of how things are looking.

Tip: You can do this quickly with barplot() and table(), or you could create try using ggplot() and looking into geom_bar().

(We’ve made the solutions visible immediately for this question, because they contain some useful little bits of code for re-ordering the months)

Solution

The quick and dirty way:

barplot(table(survey_data$birthmonth))

A ggplot option:

ggplot(data = survey_data, aes(x = birthmonth)) +
    geom_bar() +
    labs(x = "- Birth Month -")

To order the months correctly, we need to make the “birthmonth” variable a factor (a categorical variable), and order the levels:

# we can make use of month.abb which is the abbreviated months in R
month.abb

##  [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"

# but we want them all lower case:
tolower(month.abb)

##  [1] "jan" "feb" "mar" "apr" "may" "jun" "jul" "aug" "sep" "oct" "nov" "dec"

#This is one way to do it
survey_data$birthmonth <- factor(survey_data$birthmonth, levels = tolower(month.abb), ordered=T)

# and this is another (the tidyverse way):
survey_data <-
  survey_data %>%
    mutate(
      birthmonth = fct_relevel(factor(birthmonth), tolower(month.abb))
    )

Now when we plot it, they’re in the right order!

ggplot(data = survey_data, aes(x = birthmonth)) +
    geom_bar() +
    labs(x = "- Birth Month -")

Question A2

We’re going to perform a statistical test to assess the extent to which our data conforms to the hypothesis that people are no more likely to be born on one month than another.

Under this hypothesis, what would be the proportional breakdown of observed births in each of the months?

Solution

If people are no more likely to be born in one month than another, then we would expect the same proportion of observed births in each month.
There are 12 months, so we would expect \(\frac{1}{12}\) observations in each month.

We can write these as: \[ \begin{align} & p_{jan} = 1/12 \\ & p_{feb} = 1/12 \\ & ... \\ & p_{dec} = 1/12 \\ \end{align} \]

Question A3

Given that there are there are 390 people who entered their birth-month, how many observations would we expect to find with a birthday in January? And in February? … and so on?

Solution

We would expect \(\frac{1}{12} \times\) 390 = 32.5 observations in each month.

The test procedure we are about to apply is called a chi-square goodness-of-fit test.

It applies to a categorical variable, and the null hypothesis asserts specific values for the population proportion in each category. The alternative hypothesis simply states that at least one of the population proportions is not as specified in the null hypothesis.

As always, the test statistic measures how far the observed sample results deviate from what is expected if the null hypothesis is true. With a chi-square test, you construct the test statistic by comparing the observed sample counts in each category to the expected counts under the null hypothesis.

The test-statistic (denoted \(\chi^2\), spelled chi-square, pronounced “kai-square”) is obtained by adding up the standardized squared deviations in each category:
\[ \chi^2 = \sum_{i} \frac{(\text{Observed}_i - \text{Expected}_i)^2}{\text{Expected}_i} \]

Question A4

The code below creates counts for each month:

survey_data %>%
    group_by(birthmonth) %>%
    summarise(
        observed = n()
    )

(A shortcut for this would be survey_data %>% count(birthmonth))

Add to the code above to create columns showing:

• the expected counts \(E_i\)
• observed-expected (\(O_i - E_i\))
• the squared differences \((O_i - E_i)^2\)
• the standardised square differences \(\frac{(O_i - E_i)^2}{E_i}\)

Then calculate the \(\chi^2\) statistic (the sum of the standardised squared differences).
If your observed counts matched the expected counts perfectly, what would the \(\chi^2\) statistic be?

Solution

# the nrow() function tells us how many rows there are in the dataset!
totaln <- nrow(survey_data)

chi_table <- 
    survey_data %>%
    group_by(birthmonth) %>%
    summarise(
        observed = n(),
        expected = totaln/12,
        diff = observed-expected,
        sq_diff = diff^2,
        std_sq_diff = sq_diff / expected
    )
chi_table

## # A tibble: 12 x 6
##    birthmonth observed expected  diff sq_diff std_sq_diff
##    <ord>         <int>    <dbl> <dbl>   <dbl>       <dbl>
##  1 jan              36     32.5   3.5   12.2      0.377  
##  2 feb              30     32.5  -2.5    6.25     0.192  
##  3 mar              31     32.5  -1.5    2.25     0.0692 
##  4 apr              23     32.5  -9.5   90.2      2.78   
##  5 may              32     32.5  -0.5    0.25     0.00769
##  6 jun              31     32.5  -1.5    2.25     0.0692 
##  7 jul              34     32.5   1.5    2.25     0.0692 
##  8 aug              32     32.5  -0.5    0.25     0.00769
##  9 sep              36     32.5   3.5   12.2      0.377  
## 10 oct              38     32.5   5.5   30.2      0.931  
## 11 nov              36     32.5   3.5   12.2      0.377  
## 12 dec              31     32.5  -1.5    2.25     0.0692

And we can calculate our \(\chi^2\) test statistic by simply summing the values in the last column we created:

sum(chi_table$std_sq_diff)

## [1] 5.323077

If all our observed counts are equal to our expected counts, then the diff column above will be all \(0\), and \(0^2=0\), and \(\frac{0}{E_i}\) will be \(0\). So \(\chi^2\) will be \(0\).

Question A5

You can see the distribution of \(\chi^2\) distributions with different degrees of freedom below.

Figure 1: Chi-Square Distributions

We can find out the proportion of the distribution which falls to either side of a given value of \(\chi^2\) using pchisq(). We need to give it our calculated \(\chi^2\) statistic, our degrees of freedom (df), which is equal to the number of categories minus 1. We also need to specify whether we want the proportion to the left (lower.tail=TRUE) or to the right (lower.tail=FALSE).

1. Using pchisq(), calculate the probability of observing a \(\chi^2\) statistic as least as extreme as the one we have calculated.
   
2. Check that these results match with those provided by R’s built-in function: chisq.test(table(survey_data$birthmonth)).

Solution

sum(chi_table$std_sq_diff)

## [1] 5.323077

pchisq(sum(chi_table$std_sq_diff), df = 11, lower.tail = FALSE)

## [1] 0.9145164

chisq.test(table(survey_data$birthmonth))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(survey_data$birthmonth)
## X-squared = 5.3231, df = 11, p-value = 0.9145

Question A6

Which months of year had the highest contributions to the chi-square test statistic?

Hint: think about your standardised squared deviations.

Solution

Standardized squared deviations

One possible way to answer this question is to look at the individual contribution of each category to the \(\chi^2\) statistic. We computed these values in Question A4.

chi_table %>%
  select(birthmonth, std_sq_diff)

## # A tibble: 12 x 2
##    birthmonth std_sq_diff
##    <ord>            <dbl>
##  1 jan            0.377  
##  2 feb            0.192  
##  3 mar            0.0692 
##  4 apr            2.78   
##  5 may            0.00769
##  6 jun            0.0692 
##  7 jul            0.0692 
##  8 aug            0.00769
##  9 sep            0.377  
## 10 oct            0.931  
## 11 nov            0.377  
## 12 dec            0.0692

From the barplot we created earlier on, we can see which months make up higher/lower proportions than expected:

ggplot(chi_table, aes(x = birthmonth, y = observed/nrow(survey_data))) +
  geom_col(fill = 'lightblue') +
  geom_hline(yintercept = 1/12, color = 'red') +
  theme_classic(base_size = 15)

Pearson residuals

Equivalently, you could answer by looking at Pearson residuals:

chisq.test(table(survey_data$birthmonth))$residuals

## 
##        jan        feb        mar        apr        may        jun        jul 
##  0.6139406 -0.4385290 -0.2631174 -1.6664102 -0.0877058 -0.2631174  0.2631174 
##        aug        sep        oct        nov        dec 
## -0.0877058  0.6139406  0.9647638  0.6139406 -0.2631174

The greatest absolute values are for apr and oct, showing that for these months the deviations from expected to observed were the greatest.

Question A7

According to the internet (that reliable source of information!), 76% of people in the world have brown eyes, 10% have blue, 5% hazel, 5% amber, 2% green, 1% grey, and 1% have some other eye colouring (red/violet/heterochromia).

Perform a \(\chi^2\) goodness of fit test to assess the extent to which our sample of students conform to this theorised distribution of eye-colours.

• Hint: chisq.test(..., p = c(?,?,?,...) ).

Solution

Let’s look at the observed counts:

table(survey_data$eyecolour)

## 
## amber  blue brown green  grey hazel other 
##     3    92   176    52     8    43    15

Our theoretical probabilities of different eye colours must match the order in the table which we give chisq.test(). They must also always sum to 1.

chisq.test(table(survey_data$eyecolour), p = c(.05,.1,.76,.02,.01,.05,.01))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(survey_data$eyecolour)
## X-squared = 450.74, df = 6, p-value < 2.2e-16

Question A8

What are the observed proportions of our sample with each eyecolour?
Can you figure out how to use the prop.table() function?

Solution

From the help documentation (?prop.table()), we see that we can pass prop.table() the argument x, which needs to be a table.

prop.table(table(survey_data$eyecolour))*100

## 
##      amber       blue      brown      green       grey      hazel      other 
##  0.7712082 23.6503856 45.2442159 13.3676093  2.0565553 11.0539846  3.8560411

barplot(prop.table(table(survey_data$eyecolour))*100)

\(\chi^2\) Test of Independence

Research Question: Can telling a joke affect whether or not a waiter in a coffee bar receives a tip from a customer?
A study published in the Journal of Applied Social Psychology¹ investigated this question at a coffee bar of a famous seaside resort on the west Atlantic coast of France. The waiter randomly assigned coffee-ordering customers to one of three groups. When receiving the bill, one group also received a card telling a joke, another group received a card containing an advertisement for a local restaurant, and a third group received no card at all.

The data are available at https://uoepsy.github.io/data/TipJoke.csv.
The dataset contains the variables:

• Card: None, Joke, Ad.
• Tip: 1 = The customer left a tip, 0 = The customer did not leave tip.

Question B1

Produce a plot and a table to display the relationship between whether or not the customer left a tip, and what (if any) card they received alongside the bill.

Solution

tipjoke <- read_csv('https://uoepsy.github.io/data/TipJoke.csv')

table(tipjoke$Card, tipjoke$Tip)

##       
##         0  1
##   Ad   60 14
##   Joke 42 30
##   None 49 16

plot(table(tipjoke$Card, tipjoke$Tip))

Question B2

What would you expect the cell counts to look like if there were no relationship between what the waiter left and whether or not the customer tipped?

Solution

In total, 60 customers tipped (14+30+16), and 151 did not. So overall, 0.28 (\(\frac{60}{(60+151)}\)) of customers tip.

74 customers got an Ad card, 72 customers got a Joke, and 65 got None. If this were independent of whether or not they left a tip, we would expect equal proportions of tippers in each group.
So we would expect 0.28 of each group to leave a tip.

You can think about observed vs expected by looking at the two-way table along with the marginal row and column totals given:

	0	1
Ad			74
Joke			72
None			65
	151	60	211

For a given cell of the table we can calculate the expected count as \(\text{row total} \times \frac{\text{column total}}{\text{samplesize}}\):

Expected:

	0	1
Ad	52.96	21.04	74
Joke	51.53	20.47	72
None	46.52	18.48	65
	151.00	60.00	211

(If you’re wondering how we do this in R, we saw in the lectures briefly a complicated bit of code using %o% which could do this for us):

t <- tipjoke %>%
  select(Card, Tip) %>% table()

e <- rowSums(t) %o% colSums(t) / sum(t)
e

##             0        1
## Ad   52.95735 21.04265
## Joke 51.52607 20.47393
## None 46.51659 18.48341

Question B3

Just like we gave the chisq.test() function a table of observed frequencies when we conducted a goodness of fit test in earlier exercises, we can give it a two-way table of observed frequencies to conduct a test of independence.

Try it now.

Why is the degrees of freedom 2 (df = 2)?
Degrees of freedom for a \(\chi^2\) test of independence is calculated as:
\[ \begin{align} & df = (r - 1)(c - 1) \\ & \text{Where:} \\ & r = \text{number of rows} \\ & c = \text{number of columns} \\ \end{align} \] Why? Well, remember that the degrees of freedom is the number of values that are free to vary as we estimate parameters. In a \(3 \times 2\) table like the one we have for Cards \(\times\) Tips, the degrees of freedom is the number of cells in the table that can vary before we can simply calculate the values of the other cells (where we’re constrained by the need to sum to our row/column totals).

Solution

chisq.test(table(tipjoke$Card, tipjoke$Tip))

## 
##  Pearson's Chi-squared test
## 
## data:  table(tipjoke$Card, tipjoke$Tip)
## X-squared = 9.9533, df = 2, p-value = 0.006897

More RMarkdown

Question

Can you create an RMarkdown document which:

1. Reads in the https://uoepsy.github.io/data/TipJoke.csv data.
2. Conducts and reports a \(\chi^2\) test of independence examining whether telling a joke affect whether or not a waiter in a coffee bar receives a tip from a customer.
3. Successfully compiles (“knits”) into an .html file.

Revisiting NHST

Now that you have performed a number of different types of statistical hypothesis test, it is worth revisiting the general concept in order to consolidate what you’ve been doing.

• Step 1. We have been starting by considering what a given statistic is likely to be if a given hypothesis (the null) were true.
  • For the \(t\)-tests, if the null hypothesis is true (there is no difference between group means/between our observed mean and some value), then our \(t\)-statistics (if we could do our study loads of times) will mainly fall around 0, and follow a \(t\)-distribution. The precise \(t\)-distribution depends on the degrees of freedom, which in turn depends on how much data we have.
    
  • For the \(\chi^2\) tests, if the null hypothesis is true and there is no difference between the observed and expected frequencies, then our \(\chi^2\)-statistics will follow the \(\chi^2\) distribution (i.e., with 2 categories, most of them will be between 0 and 2, with fewer falling >2, see the yellow line in Figure 1).
    
• Step 2. We calculate our statistic from our observed data.
• Step 3. We ask what the probability is of getting a statistic at least as extreme as we get from Step 2, assuming the null hypothesis we stated in Step 1.

If you’re finding the programming easy, but the statistical concepts difficult

Another way which might help to think about this is that if we can make a computer do something over and over again, we can do stats! You may already be familiar with this idea from exercises with the function replicate()!

1. make the computer generate random data, based on some null hypothesis. Do it lots of times.

2. what proportion of the simulations produce results similar to the observed data (i.e., as extreme or more extreme)? This is \(p\). The only difference between this and “statistics” is that we calculate \(p\) using math, rather than having to generate random data.

Statistical vs Practical Significance

Question C1

An agricultural company is testing out a new fertiliser they have developed to improve tomato growth.

They plant 1000 seeds (taken from the same tomato plant) in the same compost and place them in positions with the same amount of sunlight. 500 of the plants receive 100ml of water daily, and the other 500 receive a 100ml of the fertiliser mixed with water. After 100 days, they measure the height of all the tomato plants (in cm).

You can find the data at https://uoepsy.github.io/data/tomatogrowth.csv.

Conduct the appropriate test to determine whether the fertiliser provides a statistically significant improvement to tomato plant growth.

Solution

tomato <- read_csv("https://uoepsy.github.io/data/tomatogrowth.csv")
summary(tomato)

##      height         group          
##  Min.   : 84.8   Length:1000       
##  1st Qu.:109.3   Class :character  
##  Median :116.0   Mode  :character  
##  Mean   :115.8                     
##  3rd Qu.:123.1                     
##  Max.   :144.8

Our alternative hypothesis is that the difference in means \((treatment - control)\) is greater than 0 (i.e., it improves growth). The t.test() function will use alphabetical ordering of the group variable, so if we say alternative="less" then it is the direction we want \((control - treatment < 0)\):

t.test(tomato$height ~ tomato$group, alternative = "less")

## 
##  Welch Two Sample t-test
## 
## data:  tomato$height by tomato$group
## t = -2.0085, df = 997.97, p-value = 0.02243
## alternative hypothesis: true difference in means between group control and group treatment is less than 0
## 95 percent confidence interval:
##        -Inf -0.2296311
## sample estimates:
##   mean in group control mean in group treatment 
##                115.1955                116.4692

Or, we can set the order by making it a factor and setting the levels:

tomato %>% mutate(
  group = factor(group, levels=c("treatment","control"))
) -> tomato

t.test(tomato$height ~ tomato$group, alternative = "greater")

## 
##  Welch Two Sample t-test
## 
## data:  tomato$height by tomato$group
## t = 2.0085, df = 997.97, p-value = 0.02243
## alternative hypothesis: true difference in means between group treatment and group control is greater than 0
## 95 percent confidence interval:
##  0.2296311       Inf
## sample estimates:
## mean in group treatment   mean in group control 
##                116.4692                115.1955

Question C2

Figure 2: Relationship between tomato plant height and crop yield

Consider the above plot. For every 5cm taller a tomato plant is, it tends to provide 1 more tomato.
Now consider the agricultural company’s situation in the previous question. Given that the fertiliser is comparitively pricey for them to manufacture, is it worth putting into production?

Solution

While the fertiliser does improve plant growth to a statistically significant (at \(\alpha=0.05\)) degree, the improvement is minimal. The difference in means is only 1.2737, which would result in about 1/4 of a tomato per plant.

Furthermore, if we take a look at the confidence interval provided by the t.test() function, we can see that a plausible value for the true difference in means is 0.23cm, which is tiny!

Things to think about

The above exercises are aimed to demonstrate that whether or not our p-value is below some set criteria (e.g., .05, .01, .001) is only a small part of the picture. There are many things which are good to remember about p-values:

1. With a big enough sample size, even a tiny tiny effect is detectable at <.05. For example, you might be interested in testing if the difference in population means across two groups is 0 (\(\mu_1 - \mu_2 = 0\). Your calculated sample difference could be \(\bar{x}_1 - \bar{x}_2 = 0.00002\) but with a very small p-value of 0.00000001. This would tell you that there is strong evidence that the observed difference in means (0.00002) is significantly different from 0. However, the practical difference, that is - the magnitude of the distance between 0.00002 and 0 - is negligible and of pretty much no interest to practitioners. This is the idea we saw in the tomato-plant example.

2. The criteria (\(\alpha\)) which we set (at .05, .01, etc.), is arbitrary.
   Two things need to be kept in mind: there is the true status of the world (which is unknown to us) and the collected data (which are available and reveal the truth only in part).
   An observed p-value smaller than the chosen alpha does not imply the true presence of an effect. The observed difference might be due to sampling variability.
   

   

   Figure 3: Two possible samples (blue dots) drawn from two populations with same mean. On the left, the selected sample shows a big difference. On the right, the sample shows no difference. Samples such as that on the left are very unlikely to happen (e.g., 5% of the time). It is for these unlikely samples that we would reject the null hypothesis incorrectly 5% of the time.

3. Even if a null hypothesis about the population is actually true, then 5% (if \(\alpha\) = 0.05) of the test-statistics computed on different samples from that population would result in a p-value <.05. If you were to obtain 100 random samples from that population, five out of the 100 p-values are likely to be <.05 even if the null hypothesis about the population was actually true.

4. If you have a single dataset, and you perform several tests of hypotheses on those data, each test comes with a probability of incorrectly rejecting the null (making a type I error) of 5%. Hence, considering the entire family of tests computed, your overall type I error probability will be larger than 5%. In simple words, this means that if you perform enough tests on the same data, you’re almost sure to reject one of the null hypotheses by mistake. This concept is known as multiple comparisons.

Further reading

There are many different competing approaches to doing statistical analyses.
In this course we are learning about what is known as the frequentist framework. Roughly speaking, this is where probabilities are defined as “long-run frequencies” (i.e., the probability of \(x\) happening over many many trials²). Even within the frequentist approach, there are different views as to how to how this definition of probability is best utilised.

The following links provide some introductory readings to some of the different schools of thought:

1. Perezgonzalez, J. D. (2015). Fisher, Neyman-Pearson or NHST? A tutorial for teaching data testing. Frontiers in Psychology, 6, 223.

2. Calin-Jageman, R. J., & Cumming, G. (2019). The new statistics for better science: ask how much, how uncertain, and what else is known. The American Statistician, 73(sup1), 271-280.

3. The correctly-used p value needs an effect size and CI - don’t worry too much about the background of this blog, but it offers some useful visualisations to show how important it is to remember about the uncertainty in our estimates.

---

1. Gueaguen, N. (2002). The Effects of a Joke on Tipping When It Is Delivered at the Same Time as the Bill. Journal of Applied Social Psychology, 32(9), 1955-1963.↩︎

2. For those of you who are interested in what alternative definitions there are, do a google search for “frequentist vs bayesian.” Be prepared that this will open a big can of worms!↩︎

---

This workbook was written by Josiah King, Umberto Noe, and Martin Corley, and is licensed under a Creative Commons Attribution 4.0 International License.