Statistics 2

Jeff Stevens

2025-05-09

Set-up

Install and load packages

First, install these packages

install.packages(c("lme4", "lmerTest"))

And load tidyverse

library(tidyverse)
library(papaja)
library(lme4)
library(easystats)

Create full data set

set.seed(1234) # this just means we'll always get the same random numbers
num_subjects <- 18
rtnorm <- function(n, mean = 0, sd = 1, min = 1, max = 5) {
  bounds <- pnorm(c(min, max), mean, sd)
  u <- runif(n, bounds[1], bounds[2])
  qnorm(u, mean, sd)}
my_data <- data.frame(subject = rep(1:num_subjects, each = 3),  # create subject column
                      gender = factor(rep(c("Female", "Male", "Nonbinary"), each = num_subjects)), # create gender column
                      age = factor(rep(c("Old", "Young"), each = num_subjects / 2)),  # create age column
                      time = rep(c(1:3), times = num_subjects / 2),  # create time column
                      response = rtnorm(num_subjects * 3, mean = 3.5, sd = 0.75),  # create response column
                      rt = round(rpois(num_subjects * 3, lambda = 100) * 10 + 1000, 1),  # create response time column
                      binary = sample(x = c(0, 1), size = num_subjects * 3, replace = TRUE))  # create binary response column
head(my_data)

  subject gender age time response   rt binary
1       1 Female Old    1 2.586046 1890      0
2       1 Female Old    2 3.706176 2060      1
3       1 Female Old    3 3.681445 2100      1
4       2 Female Old    1 3.708237 2330      1
5       2 Female Old    2 4.250138 1830      0
6       2 Female Old    3 3.740756 2120      1

Additional data

Create data set for just time 1

my_data1 <- my_data |>
  filter(time == 1)
head(my_data1)

  subject gender   age time response   rt binary
1       1 Female   Old    1 2.586046 1890      0
2       2 Female   Old    1 3.708237 2330      1
3       3 Female   Old    1 1.746820 2070      1
4       4 Female Young    1 3.505189 2140      0
5       5 Female Young    1 3.055292 1950      1
6       6 Female Young    1 4.181727 1940      0

Create within-subject/paired data

my_data_paired <- my_data |>
  filter(time != 3) |>
  pivot_wider(id_cols = subject, names_from = time, names_prefix = "time_", values_from = response)
head(my_data_paired)

# A tibble: 6 × 3
  subject time_1 time_2
    <int>  <dbl>  <dbl>
1       1   2.59   3.71
2       2   3.71   4.25
3       3   1.75   2.94
4       4   3.51   3.85
5       5   3.06   4.47
6       6   4.18   3.06

Regression

Linear relationships

my_data |>
  ggplot(aes(x = rt, y = response)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

Scaling

my_data |>
  mutate(zscore.resp = c(scale(response)),
         centered.resp = c(scale(response, scale = FALSE)))

   subject    gender   age time response   rt binary zscore.resp centered.resp
1        1    Female   Old    1 2.586046 1890      0  -1.2901200   -0.79941951
2        1    Female   Old    2 3.706176 2060      1   0.5175702    0.32071104
3        1    Female   Old    3 3.681445 2100      1   0.4776579    0.29597946
4        2    Female   Old    1 3.708237 2330      1   0.5208956    0.32277157
5        2    Female   Old    2 4.250138 1830      0   1.3954266    0.86467242
6        2    Female   Old    3 3.740756 2120      1   0.5733757    0.35529076
7        3    Female   Old    1 1.746820 2070      1  -2.6444793   -1.63864466
8        3    Female   Old    2 2.939893 2000      0  -0.7190741   -0.44557238
9        3    Female   Old    3 3.791165 1850      1   0.6547275    0.40570017
10       4    Female Young    1 3.505189 2140      0   0.1932133    0.11972414
11       4    Female Young    2 3.846467 1890      0   0.7439742    0.46100166
12       4    Female Young    3 3.561680 2120      0   0.2843785    0.17621439
13       5    Female Young    1 3.055292 1950      1  -0.5328411   -0.33017360
14       5    Female Young    2 4.471764 1900      0   1.7530924    1.08629913
15       5    Female Young    3 3.076101 1880      1  -0.4992586   -0.30936428
16       6    Female Young    1 4.181727 1940      0   1.2850242    0.79626189
17       6    Female Young    2 3.062909 2060      0  -0.5205482   -0.32255634
18       6    Female Young    3 3.019951 2010      0  -0.5898741   -0.36551394
19       7      Male   Old    1 2.821508 1910      1  -0.9101262   -0.56395731
20       7      Male   Old    2 2.939105 1890      1  -0.7203460   -0.44636049
21       7      Male   Old    3 3.127480 1820      1  -0.4163417   -0.25798507
22       8      Male   Old    1 3.098276 2160      0  -0.4634726   -0.28718960
23       8      Male   Old    2 2.741061 1920      1  -1.0399528   -0.64440403
24       8      Male   Old    3 2.182593 1940      0  -1.9412203   -1.20287206
25       9      Male   Old    1 2.905940 2060      1  -0.7738678   -0.47952517
26       9      Male   Old    2 4.110660 2250      0   1.1703346    0.72519478
27       9      Male   Old    3 3.526215 1990      0   0.2271446    0.14074956
28      10      Male Young    1 4.435598 1930      0   1.6947258    1.05013243
29      10      Male Young    2 4.165372 1950      1   1.2586295    0.77990649
30      10      Male Young    3 2.229546 1880      0  -1.8654475   -1.15591964
31      11      Male Young    1 3.398072 1990      0   0.0203456    0.01260710
32      11      Male Young    2 3.016262 2020      0  -0.5958273   -0.36920281
33      11      Male Young    3 3.102464 2110      1  -0.4567136   -0.28300139
34      12      Male Young    1 3.492437 2210      1   0.1726329    0.10697155
35      12      Male Young    2 2.805802 1990      1  -0.9354727   -0.57966324
36      12      Male Young    3 3.988283 2030      0   0.9728409    0.60281829
37      13 Nonbinary   Old    1 2.860744 2010      0  -0.8468060   -0.52472115
38      13 Nonbinary   Old    2 3.001753 2010      0  -0.6192424   -0.38371190
39      13 Nonbinary   Old    3 4.906017 1920      0   2.4538981    1.52055154
40      14 Nonbinary   Old    1 4.102393 1920      1   1.1569929    0.71692759
41      14 Nonbinary   Old    2 3.577096 2300      0   0.3092573    0.19163047
42      14 Nonbinary   Old    3 3.752572 2010      1   0.5924444    0.36710665
43      15 Nonbinary   Old    1 3.117500 2010      1  -0.4324485   -0.26796556
44      15 Nonbinary   Old    2 3.705261 1920      0   0.5160925    0.31979538
45      15 Nonbinary   Old    3 3.154576 1980      0  -0.3726131   -0.23088875
46      16 Nonbinary Young    1 3.482685 2080      0   0.1568957    0.09721999
47      16 Nonbinary Young    2 3.813095 2060      0   0.6901183    0.42763001
48      16 Nonbinary Young    3 3.451423 2050      1   0.1064435    0.06595742
49      17 Nonbinary Young    1 2.967141 1950      0  -0.6750998   -0.41832381
50      17 Nonbinary Young    2 4.001500 2010      1   0.9941707    0.61603528
51      17 Nonbinary Young    3 2.406927 1940      1  -1.5791853   -0.97853806
52      18 Nonbinary Young    1 3.113022 2020      1  -0.4396747   -0.27244329
53      18 Nonbinary Young    2 3.895621 2170      0   0.8233006    0.51015604
54      18 Nonbinary Young    3 3.487366 2100      1   0.1644498    0.10190087

Linear models

(lm_linear1 <- lm(response ~ rt, data = my_data))


Call:
lm(formula = response ~ rt, data = my_data)

Coefficients:
(Intercept)           rt  
  2.5053515    0.0004374

summary(lm_linear1)  # like aov, use summary to see regression table


Call:
lm(formula = response ~ rt, data = my_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6640 -0.3803  0.0358  0.3970  1.5608 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.5053515  1.5242206   1.644    0.106
rt          0.0004374  0.0007564   0.578    0.566

Residual standard error: 0.6236 on 52 degrees of freedom
Multiple R-squared:  0.006391,  Adjusted R-squared:  -0.01272 
F-statistic: 0.3344 on 1 and 52 DF,  p-value: 0.5655

lm_linear2 <- lm(response ~ age + rt, data = my_data)  # uses same formula syntax as aov
summary(lm_linear2)


Call:
lm(formula = response ~ age + rt, data = my_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.60409 -0.39789 -0.02336  0.40970  1.61928 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.4653239  1.5329157   1.608    0.114
ageYoung    0.1187116  0.1705917   0.696    0.490
rt          0.0004278  0.0007603   0.563    0.576

Residual standard error: 0.6267 on 51 degrees of freedom
Multiple R-squared:  0.01574,   Adjusted R-squared:  -0.02286 
F-statistic: 0.4077 on 2 and 51 DF,  p-value: 0.6673

Main effects

lm_linear2 <- lm(response ~ age + rt, data = my_data)  # uses same formula syntax as aov
summary(lm_linear2)


Call:
lm(formula = response ~ age + rt, data = my_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.60409 -0.39789 -0.02336  0.40970  1.61928 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.4653239  1.5329157   1.608    0.114
ageYoung    0.1187116  0.1705917   0.696    0.490
rt          0.0004278  0.0007603   0.563    0.576

Residual standard error: 0.6267 on 51 degrees of freedom
Multiple R-squared:  0.01574,   Adjusted R-squared:  -0.02286 
F-statistic: 0.4077 on 2 and 51 DF,  p-value: 0.6673

Interactions

lm_linear3 <- lm(response ~ age * rt, data = my_data)
summary(lm_linear3)


Call:
lm(formula = response ~ age * rt, data = my_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.60466 -0.39666 -0.02092  0.41116  1.62014 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.446e+00  1.866e+00   1.311    0.196
ageYoung     1.802e-01  3.338e+00   0.054    0.957
rt           4.374e-04  9.262e-04   0.472    0.639
ageYoung:rt -3.055e-05  1.656e-03  -0.018    0.985

Residual standard error: 0.6329 on 50 degrees of freedom
Multiple R-squared:  0.01574,   Adjusted R-squared:  -0.04331 
F-statistic: 0.2666 on 3 and 50 DF,  p-value: 0.8492

Null models

lm_linear0 <- lm(response ~ 1, data = my_data)  # use 1 for null model
summary(lm_linear0)


Call:
lm(formula = response ~ 1, data = my_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.63864 -0.38008  0.08159  0.39605  1.52055 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.38547    0.08432   40.15   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6196 on 53 degrees of freedom

Coefficients

coefficients(lm_linear3)

  (Intercept)      ageYoung            rt   ageYoung:rt 
 2.446122e+00  1.802079e-01  4.373711e-04 -3.055264e-05

Fitted values

fitted.values(lm_linear3)

       1        2        3        4        5        6        7        8 
3.272753 3.347106 3.364601 3.465197 3.246511 3.373349 3.351480 3.320864 
       9       10       11       12       13       14       15       16 
3.255258 3.496921 3.395217 3.488785 3.419626 3.399285 3.391149 3.415558 
      17       18       19       20       21       22       23       24 
3.464376 3.444035 3.281501 3.272753 3.242137 3.390844 3.285874 3.294622 
      25       26       27       28       29       30       31       32 
3.347106 3.430207 3.316490 3.411490 3.419626 3.391149 3.435899 3.448103 
      33       34       35       36       37       38       39       40 
3.484717 3.525399 3.435899 3.452171 3.325238 3.325238 3.285874 3.285874 
      41       42       43       44       45       46       47       48 
3.452075 3.325238 3.325238 3.285874 3.312117 3.472512 3.464376 3.460308 
      49       50       51       52       53       54 
3.419626 3.444035 3.415558 3.448103 3.509126 3.480649

Residuals

residuals(lm_linear3)

           1            2            3            4            5            6 
-0.686707701  0.359069748  0.316843323  0.243040072  1.003626491  0.367407201 
           7            8            9           10           11           12 
-1.604659663 -0.380971398  0.535906822  0.008267876  0.451250020  0.072894499 
          13           14           15           16           17           18 
-0.364334353  1.072479300 -0.315047739  0.766169322 -0.401467128 -0.424083797 
          19           20           21           22           23           24 
-0.459992931 -0.333648684 -0.114657284 -0.292567997 -0.544813360 -1.112028810 
          25           26           27           28           29           30 
-0.441166453  0.680452972  0.209724251  1.024108052  0.745745736 -1.161603092 
          31           32           33           34           35           36 
-0.037826393 -0.431840854 -0.382253103 -0.032962014 -0.630096729  0.536112064 
          37           38           39           40           41           42 
-0.464493881 -0.323484632  1.620142209  0.816518262  0.125020112  0.427333918 
          43           44           45           46           47           48 
-0.207738295  0.419386054 -0.157540346  0.010172832  0.348719223 -0.008885181 
          49           50           51           52           53           54 
-0.452484557  0.557465415 -1.008630628 -0.335081337  0.386495222  0.006717343

Linear mixed models

Add random effects

lmer_linear1 <- lmer(response ~ age * time + (1 | subject), 
                     data = my_data)  
summary(lmer_linear1)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ age * time + (1 | subject)
   Data: my_data

REML criterion at convergence: 103.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3463 -0.6396 -0.1627  0.6310  2.1673 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 0.0000   0.0000  
 Residual             0.3641   0.6034  
Number of obs: 54, groups:  subject, 18

Fixed effects:
              Estimate Std. Error t value
(Intercept)     2.7791     0.3072   9.046
ageYoung        1.0341     0.4345   2.380
time            0.2731     0.1422   1.920
ageYoung:time  -0.4568     0.2011  -2.271

Correlation of Fixed Effects:
            (Intr) ageYng time  
ageYoung    -0.707              
time        -0.926  0.655       
ageYoung:tm  0.655 -0.926 -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Linear mixed models

Add random slopes

lmer_linear2 <- lmer(response ~ age * time + (time | subject), 
                     data = my_data)
summary(lmer_linear2)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ age * time + (time | subject)
   Data: my_data

REML criterion at convergence: 103.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1795 -0.6739 -0.0800  0.6680  1.9714 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 subject  (Intercept) 0.23732  0.4872        
          time        0.05669  0.2381   -1.00
 Residual             0.32768  0.5724        
Number of obs: 54, groups:  subject, 18

Fixed effects:
              Estimate Std. Error t value
(Intercept)     2.7791     0.3337   8.329
ageYoung        1.0341     0.4719   2.192
time            0.2731     0.1565   1.744
ageYoung:time  -0.4568     0.2214  -2.064

Correlation of Fixed Effects:
            (Intr) ageYng time  
ageYoung    -0.707              
time        -0.944  0.667       
ageYoung:tm  0.667 -0.944 -0.707

Mixed model tests

To get p-values, use `lmerTest` package

summary(lmerTest::lmer(response ~ age * time + (1 | subject), 
                       data = my_data))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: response ~ age * time + (1 | subject)
   Data: my_data

REML criterion at convergence: 103.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3463 -0.6396 -0.1627  0.6310  2.1673 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 0.0000   0.0000  
 Residual             0.3641   0.6034  
Number of obs: 54, groups:  subject, 18

Fixed effects:
              Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)     2.7791     0.3072 50.0000   9.046  4.2e-12 ***
ageYoung        1.0341     0.4345 50.0000   2.380   0.0212 *  
time            0.2731     0.1422 50.0000   1.920   0.0606 .  
ageYoung:time  -0.4568     0.2011 50.0000  -2.271   0.0275 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) ageYng time  
ageYoung    -0.707              
time        -0.926  0.655       
ageYoung:tm  0.655 -0.926 -0.707
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Generalized linear models

glm_poisson <- glm(rt ~ gender, data = my_data, family = poisson)  # change family for generalized linear model
summary(glm_poisson)


Call:
glm(formula = rt ~ gender, family = poisson, data = my_data)

Coefficients:
                 Estimate Std. Error  z value Pr(>|z|)    
(Intercept)      7.604784   0.005260 1445.711   <2e-16 ***
genderMale      -0.002493   0.007444   -0.335    0.738    
genderNonbinary  0.008815   0.007423    1.188    0.235    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 333.68  on 53  degrees of freedom
Residual deviance: 331.12  on 51  degrees of freedom
AIC: 847.06

Number of Fisher Scoring iterations: 3

Logistic regression

summary(glm(binary ~ age, data = my_data, family = binomial))


Call:
glm(formula = binary ~ age, family = binomial, data = my_data)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.07411    0.38516   0.192    0.847
ageYoung    -0.44880    0.54933  -0.817    0.414

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 74.563  on 53  degrees of freedom
Residual deviance: 73.892  on 52  degrees of freedom
AIC: 77.892

Number of Fisher Scoring iterations: 4

Generalized linear mixed models

glmer_binomial <- glmer(binary ~ age * time + (1 | subject), 
                        data = my_data, family = binomial)  # add random effects
summary(glmer_binomial)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: binary ~ age * time + (1 | subject)
   Data: my_data

      AIC       BIC    logLik -2*log(L)  df.resid 
     82.7      92.7     -36.4      72.7        49 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.1608 -0.9283 -0.6514  0.9626  1.5351 

Random effects:
 Groups  Name        Variance Std.Dev.
 subject (Intercept) 0        0       
Number of obs: 54, groups:  subject, 18

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)     0.5216     1.0275   0.508    0.612
ageYoung       -1.8477     1.4985  -1.233    0.218
time           -0.2235     0.4747  -0.471    0.638
ageYoung:time   0.6924     0.6847   1.011    0.312

Correlation of Fixed Effects:
            (Intr) ageYng time  
ageYoung    -0.686              
time        -0.926  0.635       
ageYoung:tm  0.642 -0.929 -0.693
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Extract model statistics

papaja::apa_print(glmer_binomial)$full_result$ageYoung

[1] "$\\hat{\\beta} = -1.85$, 95\\% CI $[-4.78, 1.09]$, $z = -1.23$, $p = .218$"

papaja::apa_table(papaja::apa_print(glmer_binomial)$table)

(#tab:papaja)
Term	\(\hat{\beta}\)	95% CI	\(z\)	\(p\)
Intercept	0.52	[-1.49, 2.54]	0.51	.612
AgeYoung	-1.85	[-4.78, 1.09]	-1.23	.218
Time	-0.22	[-1.15, 0.71]	-0.47	.638
AgeYoung \(\times\) Time	0.69	[-0.65, 2.03]	1.01	.312

Extract model statistics

cocoon::format_stats(lm_linear3)

[1] "_R_^2^ = -0.043, _F_(3, 50) = 0.267, _p_ = .849"

cocoon::format_stats(lm_linear3, term = "ageYoung")

[1] "_β_ = 0.180, SE = 3.338, _t_ = 0.054, _p_ = .957"

cocoon::format_stats(glmer_binomial, term = "ageYoung")

[1] "_β_ = -1.848, SE = 1.498, _z_ = -1.233, _p_ = .218"

Model selection

AIC

AIC(lm_linear0)

[1] 104.5467

AIC(lm_linear1)

[1] 106.2005

AIC(lm_linear2)

[1] 107.6902

AIC(lm_linear3)

[1] 109.6898

Comparing models

Use `{performance}` package from `easystats`

library(performance)
compare_performance(lm_linear0, lm_linear1, lm_linear2, lm_linear3)

# Comparison of Model Performance Indices

Name       | Model | AIC (weights) | AICc (weights) | BIC (weights) |    R2
---------------------------------------------------------------------------
lm_linear0 |    lm | 104.5 (0.581) |  104.8 (0.630) | 108.5 (0.837) | 0.000
lm_linear1 |    lm | 106.2 (0.254) |  106.7 (0.244) | 112.2 (0.136) | 0.006
lm_linear2 |    lm | 107.7 (0.121) |  108.5 (0.098) | 115.6 (0.024) | 0.016
lm_linear3 |    lm | 109.7 (0.044) |  110.9 (0.029) | 119.6 (0.003) | 0.016

Name       | R2 (adj.) |  RMSE | Sigma
--------------------------------------
lm_linear0 |     0.000 | 0.614 | 0.620
lm_linear1 |    -0.013 | 0.612 | 0.624
lm_linear2 |    -0.023 | 0.609 | 0.627
lm_linear3 |    -0.043 | 0.609 | 0.633

Testing models

Test model performance with Bayes factors and/or likelihood ratio tests

test_performance(lm_linear0, lm_linear1, lm_linear2, lm_linear3)

Name       | Model |    BF |   Omega2 | p (Omega2) |       LR | p (LR)
----------------------------------------------------------------------
lm_linear0 |    lm |       |          |            |          |       
lm_linear1 |    lm | 0.162 | 5.31e-03 |      0.488 |     0.35 |  0.488
lm_linear2 |    lm | 0.176 | 9.36e-03 |      0.467 |     0.51 |  0.490
lm_linear3 |    lm | 0.136 | 5.83e-06 |      0.997 | 3.67e-04 |  0.827
Models were detected as nested (in terms of fixed parameters) and are compared in sequential order.

Assumption checks

Show residuals

ggplot(my_data, aes(x = rt, y = response)) +
  geom_segment(aes(xend = rt, yend = lm_linear1$fitted.values), color = "grey") +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

Normality of residuals

plot(density(residuals(lm_linear1)))  # plot distribution of residuals

Q-Q plot

ggplot(lm_linear1, aes(sample = rstandard(lm_linear1))) +   # plot qqplot
  geom_qq() +
  stat_qq_line()

Homogeneity of variance

Levene’s test

library(car)
leveneTest(residuals(lm_linear1) ~ my_data$age)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.5944 0.4442
      52

Check all assumptions

check_model(lm_linear1)

Effect sizes

Use the {effectsize} package

library(effectsize)

Cohen’s d

(my_d <- cohens_d(response ~ age, data = my_data1))

Cohen's d |        95% CI
-------------------------
-0.89     | [-1.85, 0.10]

- Estimated using pooled SD.

Eta squared

aov_test <- aov(response ~ gender * age, data = my_data)  # use * for both main effects and interaction
summary(aov_test)

            Df Sum Sq Mean Sq F value Pr(>F)
gender       2  0.698  0.3492   0.891  0.417
age          1  0.196  0.1959   0.500  0.483
gender:age   2  0.649  0.3245   0.828  0.443
Residuals   48 18.807  0.3918

eta_squared(aov_test)

# Effect Size for ANOVA (Type I)

Parameter  | Eta2 (partial) |       95% CI
------------------------------------------
gender     |           0.04 | [0.00, 1.00]
age        |           0.01 | [0.00, 1.00]
gender:age |           0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

Standardized coefficients

coefficients(lm_linear3)

  (Intercept)      ageYoung            rt   ageYoung:rt 
 2.446122e+00  1.802079e-01  4.373711e-04 -3.055264e-05

standardize_parameters(lm_linear3)

# Standardization method: refit

Parameter        | Std. Coef. |        95% CI
---------------------------------------------
(Intercept)      |      -0.10 | [-0.49, 0.30]
age [Young]      |       0.19 | [-0.37, 0.75]
rt               |       0.08 | [-0.26, 0.42]
age [Young] × rt |  -5.58e-03 | [-0.61, 0.60]

Convert between statistics and effect sizes

d_to_r(d = my_d)

    Cohens_d        CI     CI_low    CI_high
1 -0.4049657 0.4290568 -0.6782923 0.04962899

r_to_d(0.5)

[1] 1.154701

F_to_d(f = 1.578, df = 1, df_error = 44)

d    |        95% CI
--------------------
0.38 | [-0.22, 0.97]

Bayes factors

Use the {BayesFactor} package

library(BayesFactor)

Correlation

correlationBF(my_data$response, my_data$rt)

Bayes factor analysis
--------------
[1] Alt., r=0.333 : 0.356984 ±0%

Against denominator:
  Null, rho = 0 
---
Bayes factor type: BFcorrelation, Jeffreys-beta*

One-sample t-test

ttestBF(my_data$response, mu = 3)

Bayes factor analysis
--------------
[1] Alt., r=0.707 : 691.9069 ±0%

Against denominator:
  Null, mu = 3 
---
Bayes factor type: BFoneSample, JZS

Two-sample, independent t-test

ttestBF(formula = response ~ age, data = my_data1)

Bayes factor analysis
--------------
[1] Alt., r=0.707 : 1.295968 ±0%

Against denominator:
  Null, mu1-mu2 = 0 
---
Bayes factor type: BFindepSample, JZS

Two-sample, paired t-test

ttestBF(my_data_paired$time_1, my_data_paired$time_2, paired = TRUE)

Bayes factor analysis
--------------
[1] Alt., r=0.707 : 0.845526 ±0.02%

Against denominator:
  Null, mu = 0 
---
Bayes factor type: BFoneSample, JZS

ANOVA

anovaBF(response ~ age, data = my_data)

Bayes factor analysis
--------------
[1] age : 0.3378884 ±0.01%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

Repeated measures ANOVA

my_data <- mutate(my_data, time = as.factor(time))
anovaBF(response ~ age * time, data = my_data, whichRandom = "subject")

Bayes factor analysis
--------------
[1] age                   : 0.3378884 ±0.01%
[2] time                  : 0.3325388 ±0.01%
[3] age + time            : 0.1115577 ±1.42%
[4] age + time + age:time : 0.1466879 ±1.37%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

Regression

regressionBF(response ~ rt, data = my_data)

Bayes factor analysis
--------------
[1] rt : 0.3147939 ±0%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

Linear models

lmBF(response ~ rt, data = my_data)

Bayes factor analysis
--------------
[1] rt : 0.3147939 ±0%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

lmBF(response ~ age * time, data = my_data, whichRandom = "subject")

Bayes factor analysis
--------------
[1] age * time : 0.1460568 ±1.54%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

Let’s code!

Statistics 2

Statistics 2

Set-up

Install and load packages

Create full data set

Additional data

Create data set for just time 1

Create within-subject/paired data

Regression

Linear relationships

Scaling

Linear models

Main effects

Interactions

Null models

Coefficients

Fitted values

Residuals

Linear mixed models

Add random effects

Linear mixed models

Add random slopes

Mixed model tests

To get p-values, use lmerTest package

Generalized linear models

Logistic regression

Generalized linear mixed models

Extract model statistics

Extract model statistics

Model selection

AIC

Comparing models

Use {performance} package from easystats

Testing models

Test model performance with Bayes factors and/or likelihood ratio tests

Assumption checks

Show residuals

Normality of residuals

Q-Q plot

Homogeneity of variance

Levene’s test

Check all assumptions

Effect sizes

Effect sizes

Cohen’s d

Eta squared

Standardized coefficients

Convert between statistics and effect sizes

Bayes factors

Bayes factors

Correlation

One-sample t-test

Two-sample, independent t-test

Two-sample, paired t-test

ANOVA

Repeated measures ANOVA

Regression

Linear models

Let’s code!

To get p-values, use `lmerTest` package

Use `{performance}` package from `easystats`