Understanding R-Output

Packages

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data

Attaching package: 'TH.data'

The following object is masked from 'package:MASS':

    geyser

Goal

The goal of this activity is to make the connections between what we have done this semester with pre-packaged functions in R that are commonly used to analyze data. It is critical that we be able to understand where this information is coming from, and how to talk about it. This information is not useful if we use R as a “magic box” that gives us values.

Data

For the first couple questions, we are going to use the Cholesterol data set. These data show cholesterol reduction in patients across five treatments. The data set is called cholesterol in the multcomp package. Run ?cholesterol in the Console to read more about these data.

Below, we are making our data frame a tibble. Tibbles are usually easier to work with and give us more information upfront about our data.

data("cholesterol")

cholesterol <- cholesterol |>
  tibble()

First, we are going to look at the cholesterol reduction between the trt 1time and 2times. Let’s create a new data set to only have these two treatments.

new_data <- cholesterol |>
  filter(trt %in% c("1time", "2times")) 

new_data |>
  group_by(trt) |>
  summarise(count = n())
# A tibble: 2 × 2
  trt    count
  <fct>  <int>
1 1time     10
2 2times    10

Why would a t-test be appropriate for these data?

a t-test would be appropriate because we are working with a quantitative response and a categorical explanatory variable with 2-levels

t-test

Below is the code used to conduct a t-test for this situation. Let’s assume that we want to test if the 1time treatment had lower cholesterol reduction than 2times. Before running this code, what assumptions would we need to check?

Run the following code below.

t.test(response ~ trt, data = new_data , var.equal = FALSE, alternative = "less")

    Welch Two Sample t-test

data:  response by trt
t = -2.4097, df = 17.382, p-value = 0.01365
alternative hypothesis: true difference in means between group 1time and group 2times is less than 0
95 percent confidence interval:
       -Inf -0.9605613
sample estimates:
 mean in group 1time mean in group 2times 
             5.78197              9.22497

Note: The degrees of freedom for this test is different than what we learned at the beginning of the semester. We learned to take the smaller of \(n_1\) or \(n_2\) and subtract 1. This is a more conservative approach than the complicated formula R uses (and as seen on our slides).

Use the code chunk below to confirm the sample means for each group. Then, calculate the corresponding statistics necessary to recreate the t-statistic for this problem.

new_data |>
  group_by(trt) |>
  summarise(mean = mean(response))
# A tibble: 2 × 2
  trt     mean
  <fct>  <dbl>
1 1time   5.78
2 2times  9.22

Connection in class

pt(-2.4097, df = 18, lower.tail = TRUE)
[1] 0.01344294

Anova

Now, instead of looking at two treatments, let’s assume that we want to look at all five at once in the cholesterol group. Before running the Anova test, please list the assumptions that you would check.

Now, run the appropriate code to conduct the Anova test.

m1 <- aov(response ~ trt, data = cholesterol)
summary(m1)
            Df Sum Sq Mean Sq F value   Pr(>F)    
trt          4 1351.4   337.8   32.43 9.82e-13 ***
Residuals   45  468.8    10.4                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Let’s discuss where each of these values came from / what they represent.

df: trt_df comes from the number of cholesterol groups-1; residual df is sample size - number of groups

sums of square is all about variability. The trt sums of squares is the variability between the group means vs the overall mean. The residual error is the within group error, quantifying how spread out the observations are away from their group center.

the mean squares are the sums of squares divided by the df.

The F stat is the ratio of the sums of squares!

“By hand” way

pf(32.43, 4, 45, lower.tail = FALSE)
[1] 9.832611e-13

Chi-square

For this example, we are going to use the survey data set. This data frame contains the responses of 237 Statistics I students at the University of Adelaide to a number of questions.

survey <- tibble(survey) |> #get rid of NA values for this example
  na.omit()

glimpse(survey)
Rows: 168
Columns: 12
$ Sex    <fct> Female, Male, Male, Female, Male, Female, Male, Male, Female, F…
$ Wr.Hnd <dbl> 18.5, 19.5, 20.0, 18.0, 17.7, 17.0, 20.0, 18.5, 17.0, 19.5, 18.…
$ NW.Hnd <dbl> 18.0, 20.5, 20.0, 17.7, 17.7, 17.3, 19.5, 18.5, 17.2, 20.2, 18.…
$ W.Hnd  <fct> Right, Left, Right, Right, Right, Right, Right, Right, Right, R…
$ Fold   <fct> R on L, R on L, Neither, L on R, L on R, R on L, R on L, R on L…
$ Pulse  <int> 92, 104, 35, 64, 83, 74, 72, 90, 80, 66, 89, 74, 78, 72, 72, 64…
$ Clap   <fct> Left, Left, Right, Right, Right, Right, Right, Right, Right, Ne…
$ Exer   <fct> Some, None, Some, Some, Freq, Freq, Some, Some, Freq, Some, Fre…
$ Smoke  <fct> Never, Regul, Never, Never, Never, Never, Never, Never, Never, …
$ Height <dbl> 173.00, 177.80, 165.00, 172.72, 182.88, 157.00, 175.00, 167.00,…
$ M.I    <fct> Metric, Imperial, Metric, Imperial, Imperial, Metric, Metric, M…
$ Age    <dbl> 18.250, 17.583, 23.667, 21.000, 18.833, 35.833, 19.000, 22.333,…

Specifically, we are going to look at if a student’s writing hand has an impact on how they position their hands when they clap. See the help file ?survey for more information.

Before conducting this test, what assumptions would you check?

Independence and expected counts

Write code to look at the number of levels each each variable has. Now, write the code below to conduct the chi-square test for these data.

survey |>
  group_by(W.Hnd, Clap) |>
  summarise(n = n())
`summarise()` has grouped output by 'W.Hnd'. You can override using the
`.groups` argument.
# A tibble: 6 × 3
# Groups:   W.Hnd [2]
  W.Hnd Clap        n
  <fct> <fct>   <int>
1 Left  Left        5
2 Left  Neither     4
3 Left  Right       3
4 Right Left       23
5 Right Neither    29
6 Right Right     104

The code below can be used to calculate the expected counts.

chisq.test(survey$W.Hnd, survey$Clap)$expected
Warning in chisq.test(survey$W.Hnd, survey$Clap): Chi-squared approximation may
be incorrect
            survey$Clap
survey$W.Hnd Left   Neither     Right
       Left     2  2.357143  7.642857
       Right   26 30.642857 99.357143

What do we notice?

We notice that some of the expected counts are < 5! We do get a warning (which is nice). Not all R functions are going to tell you when assumptions are violated.