Fitting probability models to frequency data (Part I)

• Working on grading Exam #1!
• Lab #5 will be posted tonight and due next Wednesday
  • Data Manipulation with dplyr

\[ \hat{p} = \frac{\mathrm{Number \ of \ successes \ in \ sample}}{\mathrm{Total \ sample \ size}} = \frac{1}{n}\hat{X} \]

• \( \hat{X} \) is a binomial random variable
• \( \hat{p} \) is thus a scaled binomial random variable
• Standard deviation of \( \hat{X} \) is \( \sigma_{\hat{X}} = \sqrt{np(1-p)} \)
• From rules for scaled random variables, the standard deviation of \( \hat{p} \) is thus
  \[ \sigma_{\hat{p}} = \frac{1}{n}\sqrt{np(1-p)} = \sqrt{\frac{p(1-p)}{n}} \]

This is the standard error for \( \hat{p} \)!!!

Definition: The binomial test uses data to test whether a population proportion (\( p \)) matches a null expectation (\( p_{0} \)) for the proportion.

Definition: The null hypothesis \( H_{0} \) and alternative hypothesis \( H_{A} \) for a binomial test are given by:

    \( H_{0} \): Relative frequency of successes in population is \( p_{0} \).
    \( H_{A} \): Relative frequency of successes in population is not \( p_{0} \).

Do people typically use a particular ear preferentially when listening to strangers? Marzoli and Tomassi (2009) had a researcher approach and speak to strangers in a noisy nightclub. An observer scored whether the person approached turned either the left or right ear toward the questioner. Of 25 participants, 19 turned the right ear toward the questioner and 6 offered the left ear. Is this evidence of population difference from 50% for each ear?

Discuss: What is the null and alternative hypotheses?

Answer:
\[ \begin{array}{ll} H_{0}\,: & p = 0.5 \\ H_{A}\,: & p \neq 0.5 \end{array} \]

Do people typically use a particular ear preferentially when listening to strangers? Marzoli and Tomassi (2009) had a researcher approach and speak to strangers in a noisy nightclub. An observer scored whether the person approached turned either the left or right ear toward the questioner. Of 25 participants, 19 turned the right ear toward the questioner and 6 offered the left ear. Is this evidence of population difference from 50% for each ear?

Discuss: What is the observed value of the test statistic?

Answer: Number of right ears is 19 (\( \hat{X}=19 \)).

Discuss: Under the null hypothesis, calculate the probability of getting exactly 19 right ears and six left ears.

(prob <- dbinom(x = 19, size = 25, prob = 0.5))

[1] 0.005277991

\[ \mathrm{Pr[19]} = \left(\begin{array}{c}{25 \\ 19}\end{array}\right)0.5^{19}0.5^{6} = 0.005278 \]

Discuss: Give the two-tailed \( P \)-value based on your previous answer.

(pval <- 2*sum(extreme_probs))

[1] 0.0146333

Discuss: State your conclusion.

Answer: Using a significance level of \( \alpha = 0.05 \), we reject \( H_{0} \) since \( P < 0.05 \). There is evidence that more people use the right ear than the left ear when listening to a stranger in the noisy nightclub.

Use binom.test to do a binomial test! It's more accurate. If our observed test statistic is \( X = 19 \) successes out of \( n = 25 \) trials, and our null hypothesized proportion is \( p_{0} = 0.5 \), then we have:

binom.test(19,
           n = 25,
           p = 0.5)

From proportions and binomial distributions…

Note: The binomial test is an example of a goodness-of-fit test.

Definition: A goodness-of-fit test is a method for comparing an observed frequency distribution with the frequency distribution that would be expected under a simple probability model governing the occurrence of different outcomes.

Definition: A model in this case is a simplified, mathematical representation that mimics how we think a natural process works.

Assignment Problem #21

A more recent study of Feline High-Rise Syndrom (FHRS) included data on the month in which each of 119 cats fell (Vnuk et al. 2004). The data are in the accompanying table. Can we infer that the rate of cat falling varies between months of the year?

A more recent study of Feline High-Rise Syndrom (FHRS) included data on the month in which each of 119 cats fell (Vnuk et al. 2004). The data are in the accompanying table. Can we infer that the rate of cat falling varies between months of the year?

Question: What are the null and alternative hypotheses?

Answer:
     \( H_{0} \): The frequency of cats falling is the same in each month.
     \( H_{A} \): The frequency of cats falling is not the same in each month.

Observed and Expected Frequencies

We want to use dplyr for practice, so…

if (!require(dplyr)) {
  install.packages("dplyr")
  library(dplyr)
}

Now load data and store as a tibble.

mydata <- read.csv("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter08/chap08q21FallingCatsByMonth.csv") %>%
  tbl_df()

Observed and Expected Frequencies

Let's peek at the data using glimpse.

glimpse(mydata)

Rows: 119
Columns: 1
$ month <chr> "January", "January", "January", "January", "February", "Februa…

We need frequencies for months of the year.

The data in this case is in tidy form, i.e. each row is an observation (a falling cat), and each column is a measurement (month).

Question: How do we get frequencies?

Observed and Expected Frequencies

You can use the table command…

table(mydata)

mydata
    April    August  December  February   January      July      June     March 
       10        13         5         6         4        19        14         8 
      May  November   October September 
        9         7        12        12 

…but the output is not a data frame!

Question: How can we get frequencies in data frame format?

Observed and Expected Frequencies

(mytable <- mydata %>% 
  group_by(month) %>%
  summarize(obs = n()))

# A tibble: 12 x 2
   month       obs
   <chr>     <int>
 1 April        10
 2 August       13
 3 December      5
 4 February      6
 5 January       4
 6 July         19
 7 June         14
 8 March         8
 9 May           9
10 November      7
11 October      12
12 September    12

Observed and Expected Frequencies

How do we get the months in the correct order?!?!

mytable$month <- factor(mytable$month, 
                        levels = c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"))

mytable %>% arrange(month)

• First, we reorder the factor levels.
• Then, we used arrange on month.

Observed and Expected Frequencies

How do we get the months in the correct order?!?!

# A tibble: 12 x 2
   month       obs
   <fct>     <int>
 1 January       4
 2 February      6
 3 March         8
 4 April        10
 5 May           9
 6 June         14
 7 July         19
 8 August       13
 9 September    12
10 October      12
11 November      7
12 December      5

Observed and Expected Frequencies

Now let's compute the expected frequencies using mutate.

(mytable <- mytable %>%
   mutate(exp = sum(obs)/12))

Observed and Expected Frequencies

Now let's compute the expected frequencies using mutate.

# A tibble: 12 x 3
   month       obs   exp
   <fct>     <int> <dbl>
 1 April        10  9.92
 2 August       13  9.92
 3 December      5  9.92
 4 February      6  9.92
 5 January       4  9.92
 6 July         19  9.92
 7 June         14  9.92
 8 March         8  9.92
 9 May           9  9.92
10 November      7  9.92
11 October      12  9.92
12 September    12  9.92

Observed and Expected Frequencies

Data validation: Do the observed and expected frequencies have the same sum? Let's use summarize.

mytable %>% summarize(sum(obs),
                      sum(exp))

# A tibble: 1 x 2
  `sum(obs)` `sum(exp)`
       <int>      <dbl>
1        119        119

YES!

Let's make a barplot. With ggplot2, you can plot directly from this tibble. For now with base R graphics, we need a matrix.

mymatrix <- mytable %>% 
  arrange(month) %>%
  select(-month) %>%
  as.matrix()

head(mymatrix)

Let's make a barplot. With ggplot2, you can plot directly from this tibble. For now with base R graphics, we need a matrix.

     obs      exp
[1,]   4 9.916667
[2,]   6 9.916667
[3,]   8 9.916667
[4,]  10 9.916667
[5,]   9 9.916667
[6,]  14 9.916667

Rownames!!

rownames(mymatrix) <- levels(mytable$month)
head(mymatrix)

rownames(mymatrix) <- levels(mytable$month)
head(mymatrix)

         obs      exp
January    4 9.916667
February   6 9.916667
March      8 9.916667
April     10 9.916667
May        9 9.916667
June      14 9.916667

That's better! Let's plot…

tl;dr

tl;dr

• You will have to manipulate the data to create the frequency distributions!
• Consult the supplemental information for next week's homework

Definition: The \( \chi^2 \) statistic measures the discrepancy between observed frequencies from the data and expected frequencies from the null hypothesis and is given by

\[ \chi^2 = \sum_{i}\frac{(Observed_{i} - Expected_{i})^2}{Expected_{i}} \]

Definition: The \( \chi^2 \) statistic measures the discrepancy between observed frequencies from the data and expected frequencies from the null hypothesis and is given by

\[ \chi^2 = \sum_{i}\frac{(Observed_{i} - Expected_{i})^2}{Expected_{i}} \]

Discuss: What would support the null hypothesis more: a small value or large value for \( \chi^{2} \)?

Answer: Small value for \( \chi^{2} \)

\( \chi^2 \) test statistic - computed in R

(mytable %>%
  mutate(tmp = ((obs-exp)^2)/exp) %>%
  summarize(chi2 = sum(tmp)))

# A tibble: 1 x 1
   chi2
  <dbl>
1  20.7

1x1 tibble?! Here's where “everything is a data frame” doesn't work.

Solution?

\( \chi^2 \) test statistic - computed in R

Use the tidyverse package and the pull() function.

mytable %>%
  mutate(tmp = ((obs-exp)^2)/exp) %>%
  summarize(chi2 = sum(tmp)) %>%
  pull(chi2)

[1] 20.66387

The observed \( \chi^2 \) test statistic is 20.6638655.

Is that good?

Definition: The number of degrees of freedom of a \( \chi^2 \) statistic specifies which \( \chi^2 \) distribution to use as the null distribution and is given by

df = Number of categories - 1 - Number of estimated parameters from data

This is analogous to the sample size in the null distribution for a mean and proportion. Remember, for the mean and proportion, the null distribution depends on the sample size.

A more recent study of Feline High-Rise Syndrom (FHRS) included data on the month in which each of 119 cats fell (Vnuk et al. 2004). The data are in the accompanying table. Can we infer that the rate of cat falling varies between months of the year?

\( P \)-value = 0.0370255

Discuss: Conclusion?

Conclusion: Reject \( H_{0} \), i.e. there is evidence that the frequency of cats falling is not the same in each month.

Statistical tables

\[ Pr[\chi^{2} \geq \chi_{0.05,11}^2] = Pr[\chi^{2} \geq 19.675] = 0.05 \] \[ \mathrm{Observed} \ \chi^2 = 20.6638655 \]

The sampling distribution of the \( \chi^2 \) statistic follows a \( \chi^2 \) distribution only approximately. Excellent approximation if the following is true:

• None of the categories have an expected frequency less than one.
• No more than 20% of the categories have expected frequencies less than five.

Note: Still good approximation if average expected value at least five.

Note: ALL STATISTICAL TESTS HAVE ASSUMPTIONS

Assumption common to all: random sample

• Proportional model (e.g. Example 8.1)
• Binomial distribution (e.g. Practice Problem #6)
• Poisson distribution (e.g. Practice Problem #2)