Effect Sizes and Statistical Inferences

Getting ready: load packages and some data.

library(tidyverse)
library(easystats)
library(pwr)

chem <- read_csv('winter_2016_senses_valence.csv') %>% 
  filter(Modality %in% c('Taste', 'Smell'))

Dance of the p-values… in R

We’re going to simulate a large (n = 100000) set of data as a population and then take some small (n = 30) samples from it.

set.seed(100)
x <- round(rnorm(100000, mean = 50, sd = 10), 0)
mean(x)

[1] 50.01609

sd(x)

[1] 9.995032

range(x)

[1]  5 94

hist(x)

Now to sample:

d <- tibble(df01 = sample(x, 30),
            df02 = sample(x, 30),
            df03 = sample(x, 30),
            df04 = sample(x, 30),
            df05 = sample(x, 30),
            df06 = sample(x, 30),
            df07 = sample(x, 30),
            df08 = sample(x, 30),
            df09 = sample(x, 30),
            df10 = sample(x, 30))

and the plot:

d %>% pivot_longer(df01:df10, names_to = "sample", values_to = "score") %>%
  group_by(sample) %>%
  summarise(mean = mean(score),
            n = n(),
            sd = sd(score),
            se = sd/sqrt(n),
            CI = 1.96*se) %>%
  ggplot(aes(y = mean, x = sample))+
  geom_hline(yintercept = 50, linetype = 2)+
  geom_pointrange(aes(x = sample, y = mean, ymin = mean-CI, ymax = mean+CI),
                  color = "blue")+
  scale_y_continuous(limits = c(35, 65), breaks = c(40, 50, 60))+
  theme_bw()+
  coord_flip()

and now a ‘dance of the t-values’

set.seed(42)
t.test(rnorm(10), rnorm(10))

    Welch Two Sample t-test

data:  rnorm(10) and rnorm(10)
t = 1.2268, df = 13.421, p-value = 0.241
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.5369389  1.9584459
sample estimates:
 mean of x  mean of y 
 0.5472968 -0.1634567 

t.test(rnorm(10), rnorm(10))

    Welch Two Sample t-test

data:  rnorm(10) and rnorm(10)
t = 0.36582, df = 17.977, p-value = 0.7188
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.8814711  1.2531201
sample estimates:
 mean of x  mean of y 
-0.1780795 -0.3639041 

t.test(rnorm(10), rnorm(10))

    Welch Two Sample t-test

data:  rnorm(10) and rnorm(10)
t = -0.082117, df = 16.257, p-value = 0.9356
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.0340503  0.9568318
sample estimates:
  mean of x   mean of y 
-0.02021535  0.01839391 

t.test(rnorm(10), rnorm(10))

    Welch Two Sample t-test

data:  rnorm(10) and rnorm(10)
t = 2.3062, df = 17.808, p-value = 0.03335
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.06683172 1.44707262
sample estimates:
 mean of x  mean of y 
 0.5390768 -0.2178754 

set.seed(42)
t.test(rnorm(10, mean = 1), rnorm(10, mean = 0))

    Welch Two Sample t-test

data:  rnorm(10, mean = 1) and rnorm(10, mean = 0)
t = 2.9527, df = 13.421, p-value = 0.01089
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.4630611 2.9584459
sample estimates:
 mean of x  mean of y 
 1.5472968 -0.1634567 

t.test(rnorm(10, mean = 1), rnorm(10, mean = 0))

    Welch Two Sample t-test

data:  rnorm(10, mean = 1) and rnorm(10, mean = 0)
t = 2.3345, df = 17.977, p-value = 0.03137
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.1185289 2.2531201
sample estimates:
 mean of x  mean of y 
 0.8219205 -0.3639041 

t.test(rnorm(10, mean = 1), rnorm(10, mean = 0))

    Welch Two Sample t-test

data:  rnorm(10, mean = 1) and rnorm(10, mean = 0)
t = 2.0448, df = 16.257, p-value = 0.05742
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.03405031  1.95683178
sample estimates:
 mean of x  mean of y 
0.97978465 0.01839391 

t.test(rnorm(10, mean = 1), rnorm(10, mean = 0))

    Welch Two Sample t-test

data:  rnorm(10, mean = 1) and rnorm(10, mean = 0)
t = 5.3528, df = 17.808, p-value = 4.513e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.066832 2.447073
sample estimates:
 mean of x  mean of y 
 1.5390768 -0.2178754 

Planning for Power in R

pwr.r.test(r = .30, sig.level = .05, power = .80, alternative = "two.sided")

     approximate correlation power calculation (arctangh transformation) 

              n = 84.07364
              r = 0.3
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

pwr.t.test(d = 1, power = .8, type = "two.sample", alternative = "two.sided")

     Two-sample t test power calculation 

              n = 16.71472
              d = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

pwr.anova.test(k = 3, f = .25, power = .8)           

     Balanced one-way analysis of variance power calculation 

              k = 3
              n = 52.3966
              f = 0.25
      sig.level = 0.05
          power = 0.8

NOTE: n is number in each group

pwr.f2.test(u = 6, f2 = 0.15, power = .8) #whole model, based on R2

     Multiple regression power calculation 

              u = 6
              v = 90.30998
             f2 = 0.15
      sig.level = 0.05
          power = 0.8

95% CI for groups

chem %>% group_by(Modality) %>%
  summarise(mean = mean(Val),
            n = n(),
            sd = sd(Val),
            se = sd/sqrt(n),
            CI = 1.96*se)

plotting means and CIs

chem %>% group_by(Modality) %>%
  summarise(mean = mean(Val),
            n = n(),
            sd = sd(Val),
            se = sd/sqrt(n),
            CI = 1.96*se) %>%
  ggplot(aes(x = Modality, y = mean))+
  geom_point()+
  geom_errorbar(aes(ymin = mean-CI, ymax = mean+CI))+
  theme_bw()

Or we can try a different way, using stat_summary()

library(Hmisc)

chem %>% ggplot(aes(x = Modality, y = Val))+
  stat_summary(fun.data = "mean_cl_normal",
               size= 0.3, geom = "pointrange")+
  theme_bw()

NA

Cohen’s d

cohens_d(Val ~ Modality, data = chem)

Cohen's d |         95% CI
--------------------------
-1.07     | [-1.58, -0.55]

- Estimated using pooled SD.

Equivalence Testing

library(TOSTER)

chem_summary <- chem %>% group_by(Modality) %>%
  summarise(mean = mean(Val),
            n = n(),
            sd = sd(Val),
            se = sd/sqrt(n),
            CI = 1.96*se)

t_TOST(Val ~ Modality, data = chem,
       eqb = .4)

Welch Two Sample t-test

The equivalence test was non-significant, t(44.94) = 0.78, p = 0.22
The null hypothesis test was significant, t(44.94) = -4.19, p < 0.01
NHST: reject null significance hypothesis that the effect is equal to zero 
TOST: don't reject null equivalence hypothesis

TOST Results 

Effect Sizes 

Note: SMD confidence intervals are an approximation. See vignette("SMD_calcs").

Back to regression

We’ll now turn to some of the examples from Winter Ch. 11

library(broom)

icon <- read_csv('perry_winter_2017_iconicity.csv')

icon <- icon %>% mutate(SER_z = scale(SER),
                        CorteseImag_z = scale(CorteseImag),
                        Syst_z = scale(Syst),
                        Freq_z = scale(Freq))

Now to run a model and clean up the output…

icon_mdl_z <- lm(Iconicity ~ SER_z + CorteseImag_z +
                   Syst_z + Freq_z, data = icon)

tidy(icon_mdl_z) %>% 
  mutate(p.value = base::format.pval(p.value, 4),
         estimate = round(estimate, 2),
         std.error = round(std.error, 2),
         statistic = round(statistic, 2))

Or we can use the model_parameters() funcation from easystats’s parameters package.

model_parameters(icon_mdl_z)

Parameter     | Coefficient |   SE |         95% CI | t(984) |      p
---------------------------------------------------------------------
(Intercept)   |        1.15 | 0.03 | [ 1.08,  1.21] |  33.34 | < .001
SER z         |        0.53 | 0.04 | [ 0.45,  0.61] |  12.52 | < .001
CorteseImag z |       -0.42 | 0.04 | [-0.50, -0.34] | -10.72 | < .001
Syst z        |        0.06 | 0.03 | [-0.01,  0.12] |   1.79 | 0.074 
Freq z        |       -0.36 | 0.11 | [-0.58, -0.14] |  -3.20 | 0.001 

plotting time!

mycoefs <- tidy(icon_mdl_z, conf.int = T) %>%
  filter(term != '(Intercept)')

mycoefs %>% ggplot(aes(x = reorder(term, estimate), y = estimate))+
  geom_point()+
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2)+
  geom_hline(yintercept = 0, linetype = 2)+
  labs(x = NULL, y = 'Beta coefficient estimate')+
  coord_flip()+
  theme_minimal()

or…

model_parameters(icon_mdl_z) |> plot()