Hypothesis tests

Steen Flammild Harsted & Søren O’Neill

The Workflow

The Workflow

Rrrrr straight to the Danger Zone

The promise(s)




Base R - Stability is a two-edged sword

Most base R stat functions are not pipe friendly

starwars |> filter(species == "Human") |> t.test(mass ~ gender)


##
## Error: Can’t combine name and height .
## Run rlang::last_error() to see where the error occurred.
## In addition: Warning message:
## In is.na(y) : is.na() applied to non-(list or vector) of type ‘language’

Most base R stat functions are not pipe friendly

Place your data in the data argument ::: {.cell output-location=‘fragment’}

t.test(mass ~ gender, 
       data = starwars |> filter(species == "Human")
       )

    Welch Two Sample t-test

data:  mass by gender
t = -2.8744, df = 2.7808, p-value = 0.06986
alternative hypothesis: true difference in means between group feminine and group masculine is not equal to 0
95 percent confidence interval:
 -63.417398   4.648771
sample estimates:
 mean in group feminine mean in group masculine 
               56.33333                85.71765

:::


Use the _placeholder ::: {.cell output-location=‘fragment’}

starwars |>  
  filter(species == "Human") |> 
  t.test(mass ~ gender, data = _ )

    Welch Two Sample t-test

data:  mass by gender
t = -2.8744, df = 2.7808, p-value = 0.06986
alternative hypothesis: true difference in means between group feminine and group masculine is not equal to 0
95 percent confidence interval:
 -63.417398   4.648771
sample estimates:
 mean in group feminine mean in group masculine 
               56.33333                85.71765

:::

Some tests require specific input formats

…And do not have a data = argument

Does not work

chisq.test(sex, gender, data = starwars) # ERROR - no data = argument

chisq.test() expects vectors or tables, not data frames.


Works

chisq.test(starwars$sex, starwars$gender) # Accepts vectors

    Pearson's Chi-squared test

data:  starwars$sex and starwars$gender
X-squared = 77.883, df = 3, p-value < 2.2e-16
table(starwars$sex, starwars$gender) |> chisq.test() # Accepts table data

    Pearson's Chi-squared test

data:  table(starwars$sex, starwars$gender)
X-squared = 77.883, df = 3, p-value < 2.2e-16



rstatix

ggstatsplot

An extension of the ggplot2 package for creating graphics with details from statistical tests included in the information-rich plots themselves.

ggstatsplot
Summary of of available plots and types of statistical analyses
Function Plot Description
type =
parametric nonparametric robust bayesian
ggbetweenstats violin/box/dot plots Between group/condition comparisons Yes Yes Yes Yes
ggwithinstats violin/box/dot plots Within group/condition comparisons Yes Yes Yes Yes
gghistostats histograms Distribution of a numeric variable Yes Yes Yes Yes
ggdotplotstats dot plots/charts Distribution about labeled numeric variable Yes Yes Yes Yes
ggscatterstats scatterplots Correlation between two variables Yes Yes Yes Yes
ggcorrmat correlation matrices Correlations between multiple variables Yes Yes No Yes
ggpiestats pie charts Categorical data Yes Yes No Yes
ggbarstats bar charts Categorical data Yes Yes Yes Yes
ggcoefstats dot-and-whisker plots Random-effects meta-analysis Yes No Yes Yes
Read more here: https://indrajeetpatil.github.io/ggstatsplot/index.html

ggbetweenstats() - Between group/condition comparisons

sw_human <- 
  starwars |> 
  filter(species == "Human")

plot <- 
  sw_human |> 
  ggbetweenstats(x = gender, y = height,
                 type = "parametric")

plot

ggstatsplot

The output is a ggplot object that you can modify further ::: {.cell output-location=‘fragment’}

plot +
  labs(
    title = "Between gender height differences of Human Starwars characters", 
    y = "Height", 
    x = "Gender") +
  scale_y_continuous(labels = function(.) paste(.,"cm"))

:::

ggstatsplot

Changing standardized effect size measure to Cohens D. ::: {.cell output-location=‘fragment’}

sw_human |> 
  ggbetweenstats(
    x = gender, y = height,
    type = "parametric", 
    effsize.type = "biased"  # <- This lines gives Cohen D
    ) +
  labs(
    title = "Between gender height differences of Human Starwars characters", 
    y = "Height", 
    x = "Gender") +
  scale_y_continuous(labels = function(.) paste(.,"cm"))

:::

ggbetweenstats() - Between group/condition comparisons

starwars |> 
  filter(species == "Human") |> 
  
  # Plot
  ggbetweenstats(x = gender, y = height,
                 type = "nonparametric")+
  
  labs(
    title = "Between gender height differences of human Starwars characters")

ggbetweenstats() - Between group/condition comparisons

starwars |> 
  add_count(species) |>  
  filter(n > 2) |>  
  
  # Plot  
  ggbetweenstats(x = species, y = height,
                 type = "nonparametric")+
  
  labs(title = "Between species height differences of Starwars characters")

ggbarstats()

palmerpenguins::penguins|> 
  ggbarstats(
    x = species, y = sex, 
    bf.message = FALSE)

ggbarstats()

palmerpenguins::penguins|> 
  ggbarstats(
    x = species, y = island, 
    bf.message = FALSE)

ggcormat() type = “parametric”

diamonds |> 
  ggcorrmat(
    type = "parametric")

ggcormat() type = “nonparametric”

diamonds |> 
  ggcorrmat(
    type = "nonparametric")

gtsummary

Publication-ready analytical and summary tables


Remember, gtsummary defaults to nonparametric descriptive statistics and tests.


See section Tables

rstatix

Provides a simple and intuitive pipe-friendly framework, coherent with the ‘tidyverse’ design philosophy, for performing basic statistical tests, including t-test, Wilcoxon test, ANOVA, Kruskal-Wallis and correlation analyses.


Good if you know the test you want to apply.


https://rpkgs.datanovia.com/rstatix/

rstatix

starwars |> 
  add_count(species) |>  
  filter(n > 2) |>  
  kruskal_test(height ~ species)
# A tibble: 1 × 6
  .y.        n statistic    df       p method        
* <chr>  <int>     <dbl> <int>   <dbl> <chr>         
1 height    48      10.5     2 0.00519 Kruskal-Wallis

see more tests here: https://rpkgs.datanovia.com/rstatix/

rstatix

soldiers |> 
  group_by(sex) |> 
  kruskal_test(heightcm ~ race)
# A tibble: 2 × 7
  sex    .y.          n statistic    df        p method        
* <chr>  <chr>    <int>     <dbl> <int>    <dbl> <chr>         
1 Female heightcm  2056      122.     5 1.33e-24 Kruskal-Wallis
2 Male   heightcm  4152      253.     6 7.62e-52 Kruskal-Wallis

see more tests here: https://rpkgs.datanovia.com/rstatix/

statsExpressions

A consistent syntax to do statistical analysis with tidy data (in pipe-friendly manner)

The Workhorse that produces the statistical estimates and expressions for ggstatsplot

https://indrajeetpatil.github.io/statsExpressions/

statsExpressions

Functions:

  • one_sample_test()
  • two_sample_test()
  • oneway_anova()
  • corr_test()
  • contingency_table()

type =

  • “Parametric”
  • “Non-parametric”
  • “Robust”
  • “Bayesian”

statsExpressions

The output is always a tibble

starwars |> 
  filter(species == "Human") |> 
  two_sample_test(x = gender, y = height,
                       type = "parametric",
                       paired = FALSE)
# A tibble: 1 × 18
  parameter1 parameter2 mean.parameter1 mean.parameter2 statistic df.error
  <chr>      <chr>                <dbl>           <dbl>     <dbl>    <dbl>
1 height     gender                164.            182.     -3.93     7.81
  p.value method                  alternative effectsize estimate conf.level
    <dbl> <chr>                   <chr>       <chr>         <dbl>      <dbl>
1 0.00460 Welch Two Sample t-test two.sided   Hedges' g     -1.66       0.95
  conf.low conf.high conf.method conf.distribution n.obs expression
     <dbl>     <dbl> <chr>       <chr>             <int> <list>    
1    -2.80    -0.478 ncp         t                    30 <language>

statsExpressions

The expression column contains the pre-formatted text with statistical details.

stat <- starwars |> 
  filter(species == "Human") |> 
  two_sample_test(x = gender, y = height,
                  type = "parametric",
                  paired = FALSE)

stat$expression[[1]]
list(italic("t")["Welch"] * "(" * 7.81 * ")" == "-3.93", italic(p) == 
    "4.60e-03", widehat(italic("g"))["Hedges"] == "-1.66", CI["95%"] ~ 
    "[" * "-2.80", "-0.48" * "]", italic("n")["obs"] == "30")


starwars |> 
  filter(species == "Human") |> 
  ggplot(aes(x = gender, y = height))+
  geom_boxplot(width = 0.25)+
  geom_jitter(width = 0.1)+
  labs(
    subtitle = stat$expression[[1]]
  )

statsExpressions

You can acces the statistical information using pull() or $

stat |> 
  pull(method)
[1] "Welch Two Sample t-test"


The stat$method gave us a p.value of stat |> pull(p.value)

Will translate to:

The Welch Two Sample t-test gave us a p.value of 0.004597

Normal distribution?

Normality

Plots

starwars %>% 
  ggplot(aes(sample = height)) + 
  geom_qq_line(color = "red", linewidth = 1)+
  geom_qq()


::: {.cell output-location=‘column-fragment’}

starwars %>% 
  ggplot(aes(x = height)) + 
  geom_histogram(binwidth = 10)

:::
::: {.cell output-location=‘column-fragment’}

starwars %>% 
  ggplot(aes(x = height)) + 
  geom_density(fill = "grey")

:::

Normality

Test

sharpiro_test from rstatix ::: {.cell output-location=‘column-fragment’}

starwars |> shapiro_test(height)
# A tibble: 1 × 3
  variable statistic          p
  <chr>        <dbl>      <dbl>
1 height       0.875 0.00000105

:::

Normality within groups

Plots

starwars %>% 
  ggplot(aes(sample = height)) + 
  geom_qq_line(color = "red", linewidth = 1) +
  geom_qq() + 
  facet_wrap(vars(gender))


::: {.cell output-location=‘column-fragment’}

starwars %>% 
  ggplot(aes(x = height, fill = gender)) + 
  geom_histogram(binwidth = 10) +
  facet_wrap(vars(gender))

:::
::: {.cell output-location=‘column-fragment’}

starwars %>% 
  ggplot(aes(x = height, fill = gender)) + 
  geom_density() +
  facet_wrap(vars(gender), scales = "free")

:::

Normality within groups

Test

starwars %>% 
  group_by(gender) %>% 
  shapiro_test(height)
# A tibble: 3 × 4
  gender    variable statistic          p
  <chr>     <chr>        <dbl>      <dbl>
1 feminine  height       0.846 0.0153    
2 masculine height       0.857 0.00000357
3 <NA>      height       0.837 0.186

Equal variance?

starwars %>% 
  levene_test(height ~ gender)
# A tibble: 1 × 4
    df1   df2 statistic     p
  <int> <int>     <dbl> <dbl>
1     1    75      1.22 0.273

Thanks!