Pearson correlation

Purpose: Use Pearson correlation to quantify the strength and direction of a linear association between two continuous variables.

data <- read_example("pearson.csv")
cor.test(data$x, data$y, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  data$x and data$y
## t = 58.974, df = 10, p-value = 4.773e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9947113 0.9996114
## sample estimates:
##       cor 
## 0.9985654

plot(data$x, data$y, xlab = "x", ylab = "y")

Spearman rank correlation

Purpose: Use Spearman correlation to assess whether two variables have a monotonic association based on ranks, especially when normality or linearity is not plausible.

data <- read_example("spearman.csv")
cor.test(data$stress_rank, data$wellbeing_rank, method = "spearman", exact = FALSE)

## 
##  Spearman's rank correlation rho
## 
## data:  data$stress_rank and data$wellbeing_rank
## S = 571, p-value = 4.128e-12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.9964974

Kendall rank correlation

Purpose: Use Kendall correlation to evaluate rank-order association, often with smaller samples or many tied ranks.

data <- read_example("kendall.csv")
cor.test(data$adherence_rank, data$wellbeing_rank, method = "kendall", exact = FALSE)

## 
##  Kendall's rank correlation tau
## 
## data:  data$adherence_rank and data$wellbeing_rank
## z = 3.7029, p-value = 0.0002131
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.8181818

Chi-square test of association

Purpose: Use this test to examine whether two categorical variables are statistically associated in a contingency table.

data <- read_example("chiSquareAssociation.csv")
tab <- table(data$group, data$outcome)
tab

##    
##     improved not_improved
##   A        8            4
##   B        3            9

chisq.test(tab)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab
## X-squared = 2.6853, df = 1, p-value = 0.1013

Fisher exact test

Purpose: Use Fisher’s exact test to examine association between two categorical variables when sample sizes or expected cell counts are small.

data <- read_example("fisher.csv")
tab <- table(data$treatment, data$response)
tab

##    
##     no yes
##   A  2   3
##   B  4   1

fisher.test(tab)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  tab
## p-value = 0.5238
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.002560407 4.586187647
## sample estimates:
## odds ratio 
##  0.2033186

McNemar test

Purpose: Use McNemar’s test to compare paired binary responses, such as before-after classifications from the same participants.

data <- read_example("mcnemar.csv")
tab <- table(data$before, data$after)
tab

##      
##       no yes
##   no   2   6
##   yes  2   4

mcnemar.test(tab)

## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  tab
## McNemar's chi-squared = 1.125, df = 1, p-value = 0.2888

Log-linear model

Purpose: Use a log-linear model to analyse associations among several categorical variables through their frequency counts.

data <- read_example("logLinearModel.csv")
tab <- xtabs(count ~ therapy + outcome + gender, data = data)

if (has_pkg("MASS")) {
  fit <- MASS::loglm(~ therapy * outcome + gender, data = tab)
  summary(fit)
} else {
  message("Install MASS to run MASS::loglm().")
}

## Formula:
## ~therapy * outcome + gender
## attr(,"variables")
## list(therapy, outcome, gender)
## attr(,"factors")
##         therapy outcome gender therapy:outcome
## therapy       1       0      0               1
## outcome       0       1      0               1
## gender        0       0      1               0
## attr(,"term.labels")
## [1] "therapy"         "outcome"         "gender"          "therapy:outcome"
## attr(,"order")
## [1] 1 1 1 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                       X^2 df  P(> X^2)
## Likelihood Ratio 2.967460  5 0.7050026
## Pearson          2.956001  5 0.7067692

Path analysis / mediation

Purpose: Use path analysis or mediation to evaluate a theoretically specified indirect pathway between measured variables. This can estimate indirect effects, but causal interpretation still depends on design and assumptions.

data <- read_example("pathAnalysis.csv")

if (has_pkg("lavaan")) {
  model <- '
    mediator ~ a * predictor
    outcome ~ b * mediator + c_prime * predictor
    indirect := a * b
    total := c_prime + indirect
  '

  fit <- lavaan::sem(model, data = data, se = "bootstrap", bootstrap = 5000)
  summary(fit, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE,ci = TRUE)
} else {
  message("Install lavaan to run the mediation model.")
}

## Length  Class   Mode 
##      1 lavaan     S4

#Extract path coefficients

lavaan::parameterEstimates(
  fit,
  standardized = TRUE,
  ci = TRUE
)

##         lhs op              rhs    label     est    se      z pvalue ci.lower
## 1  mediator  ~        predictor        a   0.968 0.016 61.454  0.000    0.938
## 2   outcome  ~         mediator        b   0.867 0.220  3.935  0.000    0.332
## 3   outcome  ~        predictor  c_prime   0.147 0.210  0.701  0.483   -0.196
## 4  mediator ~~         mediator            0.603 0.180  3.351  0.001    0.232
## 5   outcome ~~          outcome            0.174 0.054  3.215  0.001    0.025
## 6 predictor ~~        predictor          139.500 0.000     NA     NA  139.500
## 7  indirect :=              a*b indirect   0.840 0.218  3.855  0.000    0.316
## 8     total := c_prime+indirect    total   0.987 0.015 63.732  0.000    0.962
##   ci.upper  std.lv std.all std.nox
## 1    1.001   0.968   0.998   0.084
## 2    1.221   0.867   0.851   0.851
## 3    0.647   0.147   0.149   0.013
## 4    0.933   0.603   0.005   0.005
## 5    0.228   0.174   0.001   0.001
## 6  139.500 139.500   1.000 139.500
## 7    1.196   0.840   0.849   0.072
## 8    1.024   0.987   0.998   0.084

#Extract regression path
pe <- lavaan::parameterEstimates(
  fit,
  standardized = TRUE,
  ci = TRUE
)

pe[pe$op == "~", ]

##        lhs op       rhs   label   est    se      z pvalue ci.lower ci.upper
## 1 mediator  ~ predictor       a 0.968 0.016 61.454  0.000    0.938    1.001
## 2  outcome  ~  mediator       b 0.867 0.220  3.935  0.000    0.332    1.221
## 3  outcome  ~ predictor c_prime 0.147 0.210  0.701  0.483   -0.196    0.647
##   std.lv std.all std.nox
## 1  0.968   0.998   0.084
## 2  0.867   0.851   0.851
## 3  0.147   0.149   0.013

#R2

lavaan::inspect(fit, "r2")

## mediator  outcome 
##    0.995    0.999

Structural equation modeling

Purpose: Use SEM to test models that combine latent variables, measurement models, and structural relations among variables.

data <- read_example("structuralEquationModeling.csv")

if (has_pkg("lavaan")) {
  model <- '
    stress =~ stress1 + stress2 + stress3
    recovery =~ sleep1 + sleep2 + sleep3
    recovery ~ stress
    performance ~ recovery + stress
  '

  fit <- lavaan::sem(model, data = data)
  if (lavaan::lavInspect(fit, "converged")) {
    print(summary(fit, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE))
    print(lavaan::parameterEstimates(fit, standardized = TRUE, ci = TRUE))
  } else {
    cat("This example SEM did not converge; its estimates should not be interpreted.\n")
  }
} else {
  message("Install lavaan to run SEM.")
}

## This example SEM did not converge; its estimates should not be interpreted.

Confirmatory factor analysis

Purpose: Use confirmatory factor analysis to test whether observed indicators fit a measurement model defined before looking at the data. CFA is appropriate when the expected factors and item-factor assignments are theory-driven.

data <- read_example("confirmatoryFactorAnalysis.csv")

if (has_pkg("lavaan")) {
  model <- '
    stress =~ stress1 + stress2 + stress3
    recovery =~ recovery1 + recovery2 + recovery3
    confidence =~ conf1 + conf2 + conf3
  '

  fit <- lavaan::cfa(model, data = data)
  summary(fit, standardized = TRUE, fit.measures = TRUE)
  lavaan::fitMeasures(fit, c("cfi", "tli", "rmsea", "srmr"))
} else {
  message("Install lavaan to run CFA.")
}

##   cfi   tli rmsea  srmr 
## 0.968 0.952 0.158 0.024

Discriminant analysis

Purpose: Use discriminant analysis to classify cases into known groups using continuous predictors and to describe which predictors separate the groups.

data <- read_example("discriminantAnalysis.csv")
data$group <- factor(data$group)

if (has_pkg("MASS")) {
  fit <- MASS::lda(group ~ anxiety + motivation + reaction_time, data = data)
  fit
  pred <- predict(fit)$class
  table(Predicted = pred, Observed = data$group)
  mean(pred == data$group)
} else {
  message("Install MASS to run linear discriminant analysis.")
}

## [1] 1

One-sample t-test

Purpose: Use a one-sample t-test to compare the mean of one continuous variable against a theoretical or reference value.

data <- read_example("oneSampleT.csv")
t.test(data$score, mu = 0)

## 
##  One Sample t-test
## 
## data:  data$score
## t = 84.126, df = 13, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  51.98692 54.72736
## sample estimates:
## mean of x 
##  53.35714

One-sample z-test

Purpose: Use a one-sample z-test to compare a sample mean against a reference mean when the population standard deviation is known.

data <- read_example("oneSampleZ.csv")

mu0 <- 50
sigma <- 15
x <- data$score
n <- length(x)
z <- (mean(x) - mu0) / (sigma / sqrt(n))
p <- 2 * pnorm(-abs(z))

c(
  sample_mean = mean(x),
  null_mean = mu0,
  population_sd = sigma,
  n = n,
  z = z,
  p = p
)

##   sample_mean     null_mean population_sd             n             z 
##  1.080667e+02  5.000000e+01  1.500000e+01  3.000000e+01  2.120295e+01 
##             p 
## 8.969653e-100

One-sample Wilcoxon signed-rank test

Purpose: Use the one-sample Wilcoxon signed-rank test to compare a median or typical value against a reference value when a parametric t-test is not appropriate.

data <- read_example("oneSampleWilcoxon.csv")
wilcox.test(data$rating, mu = 0)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  data$rating
## V = 120, p-value = 0.000592
## alternative hypothesis: true location is not equal to 0

Independent-samples t-test

Purpose: Use an independent-samples t-test to compare the means of two independent groups on a continuous outcome.

data <- read_example("independentT.csv")
data$group <- factor(data$group)
t.test(score ~ group, data = data, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  score by group
## t = -8.8299, df = 17.783, p-value = 6.452e-08
## alternative hypothesis: true difference in means between group control and group training is not equal to 0
## 95 percent confidence interval:
##  -6.933593 -4.266407
## sample estimates:
##  mean in group control mean in group training 
##                   19.7                   25.3

Mann-Whitney U test

Purpose: Use the Mann-Whitney U test to compare two independent groups when the outcome is ordinal or when parametric assumptions are not suitable.

data <- read_example("mannWhitney.csv")
data$group <- factor(data$group)
wilcox.test(pain_score ~ group, data = data)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  pain_score by group
## W = 1, p-value = 0.0002051
## alternative hypothesis: true location shift is not equal to 0

Paired-samples t-test

Purpose: Use a paired-samples t-test to compare two related measurements, such as pre-post scores from the same participants.

data <- read_example("pairedT.csv")
t.test(data$baseline_score, data$followup_score, paired = TRUE)

## 
##  Paired t-test
## 
## data:  data$baseline_score and data$followup_score
## t = -19.334, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -5.660577 -4.593590
## sample estimates:
## mean difference 
##       -5.127083

shapiro.test(data$difference)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$difference
## W = 0.98296, p-value = 0.7051

Wilcoxon signed-rank test

Purpose: Use the Wilcoxon signed-rank test to compare two paired measurements when the differences are ordinal or non-normal.

data <- read_example("wilcoxon.csv")
wilcox.test(data$before, data$after, paired = TRUE)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  data$before and data$after
## V = 0, p-value = 0.002964
## alternative hypothesis: true location shift is not equal to 0

Chi-square test for one variance

Purpose: Use this test to compare a sample variance with a hypothesized population variance under normality assumptions.

data <- read_example("chiSquareVariance.csv")

n <- length(data$score)
s2 <- var(data$score)
sigma2_0 <- 100
chi2 <- (n - 1) * s2 / sigma2_0
p <- 2 * min(pchisq(chi2, df = n - 1), 1 - pchisq(chi2, df = n - 1))

c(
  chi_square = chi2,
  df = n - 1,
  p = p
)

##  chi_square          df           p 
## 21.10666667 11.00000000  0.06456529

F-test for two variances

Purpose: Use the F-test to compare whether two independent groups have equal variances.

data <- read_example("varianceFTest.csv")
data$group <- factor(data$group)
var.test(score ~ group, data = data)

## 
##  F test to compare two variances
## 
## data:  score by group
## F = 0.059761, num df = 7, denom df = 7, p-value = 0.001418
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.01196431 0.29849903
## sample estimates:
## ratio of variances 
##         0.05976065

One-way ANOVA

Purpose: Use one-way ANOVA to compare the means of three or more independent groups defined by one factor.

data <- read_example("anova.csv")
data$group <- factor(data$group)

fit <- aov(score ~ group, data = data)
summary(fit)

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## group        2 265.33  132.67    85.9 6.1e-09 ***
## Residuals   15  23.17    1.54                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(fit)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = score ~ group, data = data)
## 
## $group
##                  diff        lwr       upr     p adj
## low-high    -9.333333 -11.197034 -7.469633 0.0000000
## medium-high -5.666667  -7.530367 -3.802966 0.0000029
## medium-low   3.666667   1.802966  5.530367 0.0003548

ANCOVA

Purpose: Use ANCOVA to compare group means on one continuous outcome while adjusting for one or more continuous covariates. It is useful when groups may differ on a baseline score or another variable that should be statistically controlled.

data <- read_example("ancova.csv")
data$group <- factor(data$group)

fit <- lm(post_anxiety ~ group + baseline_anxiety, data = data)
anova(fit)

## Analysis of Variance Table
## 
## Response: post_anxiety
##                  Df  Sum Sq Mean Sq F value    Pr(>F)    
## group             2 1113.56  556.78  1262.4 < 2.2e-16 ***
## baseline_anxiety  1  365.22  365.22   828.1 < 2.2e-16 ***
## Residuals        32   14.11    0.44                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

drop1(fit, test = "F")

## Single term deletions
## 
## Model:
## post_anxiety ~ group + baseline_anxiety
##                  Df Sum of Sq     RSS     AIC F value    Pr(>F)    
## <none>                          14.11 -25.711                      
## group             2   1018.82 1032.94 124.839  1155.0 < 2.2e-16 ***
## baseline_anxiety  1    365.22  379.33  90.776   828.1 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

MANOVA

Purpose: Use MANOVA when you want to compare two or more related continuous dependent variables across groups at the same time. It first tests whether the groups differ on the combined outcome profile before considering any follow-up tests for individual outcomes.

data <- read_example("manova.csv")
data$group <- factor(data$group)

fit <- manova(
  cbind(stress_score, wellbeing_score, performance_score) ~ group,
  data = data
)
summary(fit, test = "Pillai")

##            Df  Pillai approx F num Df den Df    Pr(>F)    
## group       2 0.47607   13.641      6    262 1.783e-13 ***
## Residuals 132                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit, test = "Wilks")

##            Df   Wilks approx F num Df den Df    Pr(>F)    
## group       2 0.56426   14.354      6    260 3.935e-14 ***
## Residuals 132                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(fit)

##  Response stress_score :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## group         2  721.9  360.97  13.091 6.507e-06 ***
## Residuals   132 3639.6   27.57                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response wellbeing_score :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## group         2  975.7  487.85  17.853 1.373e-07 ***
## Residuals   132 3607.0   27.33                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response performance_score :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## group         2 1716.7  858.35  34.958 6.555e-13 ***
## Residuals   132 3241.1   24.55                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

MANCOVA

Purpose: Use MANCOVA when groups are compared on two or more related continuous dependent variables while statistically controlling one or more continuous covariates. It is useful when baseline differences or participant characteristics should be adjusted before interpreting group differences.

data <- read_example("mancova.csv")
data$group <- factor(data$group)

fit <- manova(
  cbind(stress_score, wellbeing_score, performance_score) ~ group + age_years + baseline_score,
  data = data
)
summary(fit, test = "Pillai")

##                Df  Pillai approx F num Df den Df    Pr(>F)    
## group           2 1.97801   899.69      6     60 < 2.2e-16 ***
## age_years       1 0.70799    23.44      3     29 6.654e-08 ***
## baseline_score  1 0.94963   182.25      3     29 < 2.2e-16 ***
## Residuals      31                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit, test = "Wilks")

##                Df    Wilks approx F num Df den Df    Pr(>F)    
## group           2 0.000115   891.80      6     58 < 2.2e-16 ***
## age_years       1 0.292012    23.44      3     29 6.654e-08 ***
## baseline_score  1 0.050368   182.25      3     29 < 2.2e-16 ***
## Residuals      31                                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(fit)

##  Response stress_score :
##                Df Sum Sq Mean Sq  F value    Pr(>F)    
## group           2 990.50  495.25 1034.028 < 2.2e-16 ***
## age_years       1  32.19   32.19   67.199 2.935e-09 ***
## baseline_score  1 239.47  239.47  499.981 < 2.2e-16 ***
## Residuals      31  14.85    0.48                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response wellbeing_score :
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## group           2 1628.72  814.36 1565.61 < 2.2e-16 ***
## age_years       1   26.06   26.06   50.10 5.967e-08 ***
## baseline_score  1  222.98  222.98  428.68 < 2.2e-16 ***
## Residuals      31   16.12    0.52                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response performance_score :
##                Df  Sum Sq Mean Sq  F value    Pr(>F)    
## group           2 1602.89  801.44 1203.966 < 2.2e-16 ***
## age_years       1   27.58   27.58   41.435 3.561e-07 ***
## baseline_score  1  213.20  213.20  320.277 < 2.2e-16 ***
## Residuals      31   20.64    0.67                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Repeated-measures ANOVA

Purpose: Use repeated-measures ANOVA to compare means across three or more related conditions or time points from the same participants.

data <- read_example("repeatedAnova.csv")

long_data <- data.frame(
  id = rep(data$id, each = 3),
  condition = factor(rep(c("baseline", "week4", "week8"), times = nrow(data))),
  score = c(t(data[, c("baseline", "week4", "week8")]))
)

fit <- aov(score ~ condition + Error(id / condition), data = long_data)
summary(fit)

## 
## Error: id
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  1  7.127   7.127               
## 
## Error: id:condition
##           Df Sum Sq Mean Sq
## condition  2  149.1   74.57
## 
## Error: Within
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## condition  2  32.02  16.010    7.26 0.00342 **
## Residuals 24  52.92   2.205                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Two-way ANOVA

Purpose: Use two-way ANOVA to test the main and interaction effects of two or more independent factors on a continuous outcome.

data <- read_example("twoWayAnova.csv")
data$therapy <- factor(data$therapy)
data$age_group <- factor(data$age_group)

fit <- aov(stress_score ~ therapy * age_group, data = data)
summary(fit)

##                   Df Sum Sq Mean Sq F value   Pr(>F)    
## therapy            1  469.2   469.2  79.120 1.29e-10 ***
## age_group          1  255.0   255.0  43.002 1.26e-07 ***
## therapy:age_group  1    4.2     4.2   0.712    0.404    
## Residuals         36  213.5     5.9                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Two-way repeated-measures ANOVA

Purpose: Use two-way repeated-measures ANOVA to test main and interaction effects when the same participants are measured across combinations of repeated factors.

data <- read_example("twoWayRepeatedAnova.csv")

long_data <- data.frame(
  participant = rep(data$participant, each = 6),
  load = factor(rep(rep(c("low", "high"), each = 3), times = nrow(data))),
  time = factor(rep(c("pre", "post", "followup"), times = 2 * nrow(data))),
  score = c(t(data[, c("low_pre", "low_post", "low_followup", "high_pre", "high_post", "high_followup")]))
)

fit <- aov(score ~ load * time + Error(participant / (load * time)), data = long_data)
summary(fit)

## 
## Error: participant
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 23   1008   43.82               
## 
## Error: participant:load
##           Df Sum Sq Mean Sq F value Pr(>F)    
## load       1 158.34  158.34   754.6 <2e-16 ***
## Residuals 23   4.83    0.21                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: participant:time
##           Df Sum Sq Mean Sq F value Pr(>F)    
## time       2   1265   632.5    3656 <2e-16 ***
## Residuals 46      8     0.2                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: participant:load:time
##           Df Sum Sq Mean Sq F value Pr(>F)    
## load:time  2 21.097   10.55   149.9 <2e-16 ***
## Residuals 46  3.236    0.07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kruskal-Wallis test

Purpose: Use the Kruskal-Wallis test to compare three or more independent groups when the outcome is ordinal or parametric assumptions are not suitable.

data <- read_example("kruskalWallis.csv")
data$group <- factor(data$group)
kruskal.test(rank_score ~ group, data = data)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  rank_score by group
## Kruskal-Wallis chi-squared = 15.028, df = 2, p-value = 0.0005455

Friedman test

Purpose: Use the Friedman test to compare three or more related conditions when repeated-measures ANOVA assumptions are not suitable.

data <- read_example("friedman.csv")
friedman.test(as.matrix(data[, c("condition_a", "condition_b", "condition_c")]))

## 
##  Friedman rank sum test
## 
## data:  as.matrix(data[, c("condition_a", "condition_b", "condition_c")])
## Friedman chi-squared = 18.541, df = 2, p-value = 9.418e-05

Nonparametric two-way ANOVA via aligned rank transform

Purpose: Use aligned rank transform ANOVA to test factorial effects with nonparametric data while retaining an ANOVA-like interaction framework.

data <- read_example("nonparametricTwoWayAnova.csv")
data$therapy <- factor(data$therapy)
data$age_group <- factor(data$age_group)

if (has_pkg("ARTool")) {
  fit <- ARTool::art(stress_rating ~ therapy * age_group, data = data)
  anova(fit)
} else {
  message("Install ARTool to run aligned rank transform ANOVA.")
}

## Analysis of Variance of Aligned Rank Transformed Data
## 
## Table Type: Anova Table (Type III tests) 
## Model: No Repeated Measures (lm)
## Response: art(stress_rating)
## 
##                     Df Df.res F value     Pr(>F)    
## 1 therapy            1     36  72.966 3.5045e-10 ***
## 2 age_group          1     36  15.429 0.00037187 ***
## 3 therapy:age_group  1     36   0.000 1.00000000    
## ---
## Signif. codes:   0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nonparametric two-way repeated-measures ANOVA via aligned rank transform

Purpose: Use aligned rank transform methods for repeated factorial designs when the outcome does not meet parametric assumptions.

data <- read_example("nonparametricTwoWayRepeatedAnova.csv")

long_data <- data.frame(
  participant = factor(rep(data$participant, each = 6)),
  load = factor(rep(rep(c("low", "high"), each = 3), times = nrow(data))),
  time = factor(rep(c("pre", "post", "followup"), times = 2 * nrow(data))),
  stress_rating = c(t(data[, c("low_pre", "low_post", "low_followup", "high_pre", "high_post", "high_followup")]))
)

if (has_pkg("ARTool")) {
  fit <- ARTool::art(stress_rating ~ load * time + Error(participant), data = long_data)
  anova(fit)
} else {
  message("Install ARTool to run aligned rank transform ANOVA.")
}

## Analysis of Variance of Aligned Rank Transformed Data
## 
## Table Type: Repeated Measures Analysis of Variance Table (Type I) 
## Model: Repeated Measures (aov)
## Response: art(stress_rating)
## 
##             Error Df Df.res    F value Pr(>F)    
## 1 load      Withn  1    115 3.9364e+02 <2e-16 ***
## 2 time      Withn  2    115 3.0906e+02 <2e-16 ***
## 3 load:time Withn  2    115 8.7071e-30      1    
## ---
## Signif. codes:   0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Binomial test

Purpose: Use the binomial test to compare an observed number of successes in a two-level outcome with an expected probability.

data <- read_example("binomialTest.csv")
successes <- sum(data$response == "yes")
n <- length(data$response)
binom.test(successes, n, p = 0.50)

## 
##  Exact binomial test
## 
## data:  successes and n
## number of successes = 14, number of trials = 20, p-value = 0.1153
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.4572108 0.8810684
## sample estimates:
## probability of success 
##                    0.7

Chi-square goodness-of-fit test

Purpose: Use the chi-square goodness-of-fit test to compare observed category frequencies with expected frequencies from a theoretical distribution.

data <- read_example("chiSquareGoodness.csv")
observed <- table(data$preference)
observed

## 
##  A  B  C 
## 12  8  5

chisq.test(observed, p = rep(1 / length(observed), length(observed)))

## 
##  Chi-squared test for given probabilities
## 
## data:  observed
## X-squared = 2.96, df = 2, p-value = 0.2276

Linear regression

Purpose: Use linear regression to predict a continuous outcome from one or more predictors and estimate their associations. Regression coefficients should not be interpreted as causal effects without a suitable design.

data <- read_example("linearRegression.csv")

fit <- lm(
  race_time_min ~ training_volume_km_week + vo2max_ml_kg_min + age_years,
  data = data
)

summary(fit)

## 
## Call:
## lm(formula = race_time_min ~ training_volume_km_week + vo2max_ml_kg_min + 
##     age_years, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5961 -0.3601 -0.1200  0.3805  1.8638 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             81.77108    8.00704  10.212 1.37e-10 ***
## training_volume_km_week -0.04576    0.06357  -0.720   0.4780    
## vo2max_ml_kg_min        -0.65184    0.19190  -3.397   0.0022 ** 
## age_years                0.11934    0.02371   5.033 3.08e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7385 on 26 degrees of freedom
## Multiple R-squared:  0.9887, Adjusted R-squared:  0.9874 
## F-statistic: 759.4 on 3 and 26 DF,  p-value: < 2.2e-16

par(mfrow = c(2, 2))
plot(fit)

par(mfrow = c(1, 1))

Binary logistic regression

Purpose: Use binary logistic regression to predict a two-category outcome and estimate effects as log-odds or odds ratios.

data <- read_example("logisticRegression.csv")
data$injury <- factor(data$injury, levels = c("no", "yes"))

fit <- glm(
  injury ~ training_load_au + previous_injury + recovery_score + age_years,
  data = data,
  family = binomial
)

summary(fit)

## 
## Call:
## glm(formula = injury ~ training_load_au + previous_injury + recovery_score + 
##     age_years, family = binomial, data = data)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)         9.293e+02  4.908e+06       0        1
## training_load_au   -5.776e-01  4.081e+03       0        1
## previous_injuryyes -6.158e+00  1.336e+05       0        1
## recovery_score     -8.640e+00  3.829e+04       0        1
## age_years          -1.526e+00  4.149e+04       0        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4.3860e+01  on 31  degrees of freedom
## Residual deviance: 4.5513e-10  on 27  degrees of freedom
## AIC: 10
## 
## Number of Fisher Scoring iterations: 25

exp(coef(fit))

##        (Intercept)   training_load_au previous_injuryyes     recovery_score 
##                Inf       0.5612239961       0.0021169728       0.0001769126 
##          age_years 
##       0.2173359948

Multinomial logistic regression

Purpose: Use multinomial logistic regression to predict an outcome with more than two unordered categories.

data <- read_example("multinomialRegression.csv")

if (has_pkg("nnet")) {
  data$treatment_outcome <- relevel(factor(data$treatment_outcome), ref = "completion")
  fit <- nnet::multinom(
    treatment_outcome ~ baseline_symptom_score + motivation_score + previous_therapy + age_years,
    data = data
  )
  summary(fit)
} else {
  message("Install nnet to run multinomial logistic regression.")
}

## # weights:  18 (10 variable)
## initial  value 49.437553 
## iter  10 value 6.671487
## iter  20 value 0.022366
## iter  30 value 0.004528
## iter  40 value 0.002878
## iter  50 value 0.002641
## iter  60 value 0.001961
## iter  70 value 0.001823
## iter  80 value 0.001242
## iter  90 value 0.001146
## iter 100 value 0.000850
## final  value 0.000850 
## stopped after 100 iterations

## Call:
## nnet::multinom(formula = treatment_outcome ~ baseline_symptom_score + 
##     motivation_score + previous_therapy + age_years, data = data)
## 
## Coefficients:
##         (Intercept) baseline_symptom_score motivation_score previous_therapyyes
## dropout   -33.98122               5.509793        -27.93061           -8.623852
## switch     19.85701               2.638611        -12.04692           -5.902525
##         age_years
## dropout -4.558764
## switch  -2.382015
## 
## Std. Errors:
##         (Intercept) baseline_symptom_score motivation_score previous_therapyyes
## dropout    84.82273               47.49466         255.3361            279.6206
## switch     67.97662               44.96074         222.8516            276.9262
##         age_years
## dropout  91.49469
## switch   92.43814
## 
## Residual Deviance: 0.001700612 
## AIC: 20.0017

Linear mixed model

Purpose: Use a linear mixed model for continuous outcomes when observations are clustered or repeated, allowing random effects such as participant-level intercepts.

data <- read_example("linearMixedModel.csv")
data$participant <- factor(data$participant)
data$condition <- factor(data$condition)
data$team <- factor(data$team)

if (has_pkg("lme4") && has_pkg("lmerTest")) {
  fit <- lmerTest::lmer(score ~ condition + week + (1 | participant), data = data)
  summary(fit)

  if (has_pkg("performance")) {
    performance::r2_nakagawa(fit)
  }
} else {
  message("Install lme4 and lmerTest to run the linear mixed model.")
}

## # R2 for Mixed Models
## 
##   Conditional R2: 0.821
##      Marginal R2: 0.528

Generalized linear mixed model

Purpose: Use a generalized linear mixed model for non-continuous outcomes, such as binary responses, when observations are clustered or repeated.

data <- read_example("generalizedMixedModel.csv")
data$participant <- factor(data$participant)
data$condition <- factor(data$condition)
data$success_num <- ifelse(data$success %in% c("yes", "1", 1, "success"), 1, 0)

if (has_pkg("lme4")) {
  fit <- lme4::glmer(
    success_num ~ condition + week + (1 | participant),
    data = data,
    family = binomial
  )
  summary(fit)
  exp(lme4::fixef(fit))
} else {
  message("Install lme4 to run the generalized linear mixed model.")
}

##           (Intercept) conditionintervention                  week 
##          1.034478e-05          2.013264e+00          2.478716e+02

Ordinal mixed model

Purpose: Use an ordinal mixed model for ordered categorical outcomes collected in clustered or repeated-measures designs.

data <- read_example("ordinalMixedModel.csv")
data$participant <- factor(data$participant)
data$condition <- factor(data$condition)
data$rating <- ordered(data$rating)

if (has_pkg("ordinal")) {
  fit <- ordinal::clmm(rating ~ condition + day + (1 | participant), data = data)
  summary(fit)
  exp(coef(fit))
} else {
  message("Install ordinal to run the ordinal mixed model.")
}

##                   1|2                   2|3                   3|4 
##          3.095034e+00          6.828107e+06          1.059881e+15 
##                   4|5 conditionintervention                   day 
##          2.264493e+21          4.017512e+04          9.238886e+04

Factor analysis

Purpose: Use factor analysis to identify latent dimensions that explain patterns of correlations among several observed variables.

data <- read_example("factorAnalysis.csv")
items <- c(
  "stress_item_1", "stress_item_2", "stress_item_3",
  "mood_item_1", "mood_item_2", "mood_item_3",
  "sleep_item_1", "sleep_item_2"
)

if (has_pkg("psych")) {
  psych::fa.parallel(data[, items], fa = "fa")
  fit <- psych::fa(data[, items], nfactors = 3, rotate = "oblimin")
  print(fit$loadings, cutoff = .30)
} else {
  fit <- factanal(data[, items], factors = 3, rotation = "promax")
  print(fit)
}

## Parallel analysis suggests that the number of factors =  1  and the number of components =  NA

## 
## Loadings:
##               MR1    MR2    MR3   
## stress_item_1  0.969              
## stress_item_2  0.808              
## stress_item_3         0.789       
## mood_item_1   -0.969              
## mood_item_2   -0.735              
## mood_item_3   -0.492 -0.454       
## sleep_item_1   1.069              
## sleep_item_2          1.003       
## 
##                  MR1   MR2   MR3
## SS loadings    4.478 1.894 0.217
## Proportion Var 0.560 0.237 0.027
## Cumulative Var 0.560 0.797 0.824

Network analysis

Purpose: Use network analysis to explore direct relations among several variables while keeping each variable as a node. The resulting edges should be interpreted as conditional associations, not as evidence of causality.

data <- read_example("networkAnalysis.csv")
items <- c("stress", "sleep_quality", "rumination", "mood", "motivation", "performance")

cor_mat <- cor(data[, items])
round(cor_mat, 2)

##               stress sleep_quality rumination  mood motivation performance
## stress          1.00         -0.98       0.95 -0.96      -0.96       -0.98
## sleep_quality  -0.98          1.00      -0.96  0.98       0.96        0.98
## rumination      0.95         -0.96       1.00 -0.97      -0.95       -0.98
## mood           -0.96          0.98      -0.97  1.00       0.95        0.98
## motivation     -0.96          0.96      -0.95  0.95       1.00        0.98
## performance    -0.98          0.98      -0.98  0.98       0.98        1.00

if (has_pkg("qgraph")) {
  qgraph::qgraph(
    cor_mat,
    layout = "spring",
    labels = items,
    theme = "colorblind"
  )
} else {
  message("Install qgraph to draw the network plot.")
}

Cluster analysis

Purpose: Use cluster analysis to group people or objects into relatively homogeneous clusters based on similarity across variables.

data <- read_example("clusterAnalysis.csv")
variables <- c("training_load", "recovery_score", "sprint_time", "jump_height", "sleep_time")

scaled <- scale(data[, variables])
d <- dist(scaled, method = "euclidean")
fit <- hclust(d, method = "ward.D2")
plot(fit)

clusters <- cutree(fit, k = 3)
table(clusters)

## clusters
## 1 2 3 
## 4 4 4

set.seed(1)
kmeans_fit <- kmeans(scaled, centers = 3, nstart = 25)
table(kmeans_fit$cluster)

## 
## 1 2 3 
## 4 4 4

Multidimensional scaling

Purpose: Use multidimensional scaling to represent similarities or distances among objects in a low-dimensional visual space.

data <- read_example("multidimensionalScaling.csv")
variables <- c("calm", "joy", "anger", "sadness", "fear")

d <- dist(scale(data[, variables]))
fit <- cmdscale(d, k = 2, eig = TRUE)

plot(fit$points, xlab = "Dimension 1", ylab = "Dimension 2")
text(fit$points, labels = data$object, pos = 3)

if (has_pkg("MASS")) {
  nonmetric_fit <- MASS::isoMDS(d, k = 2)
  nonmetric_fit$stress
}

## initial  value 3.393709 
## final  value 0.000000 
## converged

## [1] 1.823673e-14

Odds ratio meta-analysis

Purpose: Use this meta-analysis when primary studies report binary events in two groups and the effect is summarized as an odds ratio.

data <- read_example("metaOddsRatio.csv")

if (has_pkg("metafor")) {
  fit <- metafor::rma(
    yi = log_odds_ratio,
    vi = variance,
    mods = ~ moderator_code,
    data = data,
    method = "REML"
  )
  summary(fit)
  predict(fit, transf = exp)
  metafor::forest(fit, slab = data$study_label)
} else {
  message("Install metafor to run this meta-analysis.")
}

Correlation coefficient meta-analysis

Purpose: Use this meta-analysis to synthesize correlations across studies, typically after transforming correlations to Fisher’s z scale.

data <- read_example("metaCorrelation.csv")

if (has_pkg("metafor")) {
  es <- metafor::escalc(measure = "ZCOR", ri = correlation, ni = sample_size, data = data)
  fit <- metafor::rma(yi, vi, mods = ~ moderator, data = es, method = "REML")
  summary(fit)
  predict(fit, transf = metafor::transf.ztor)
  metafor::forest(fit, slab = data$study_label)
} else {
  message("Install metafor to run this meta-analysis.")
}

Mean difference meta-analysis

Purpose: Use this meta-analysis when studies report group means, standard deviations, and sample sizes for a continuous outcome.

data <- read_example("metaMeanDifferences.csv")

if (has_pkg("metafor")) {
  es <- metafor::escalc(
    measure = "SMD",
    m1i = mean_treatment,
    sd1i = sd_treatment,
    n1i = n_treatment,
    m2i = mean_control,
    sd2i = sd_control,
    n2i = n_control,
    data = data
  )

  fit <- metafor::rma(yi, vi, mods = ~ moderator, data = es, method = "REML")
  summary(fit)
  metafor::forest(fit, slab = data$study_label)
} else {
  message("Install metafor to run this meta-analysis.")
}

Effect size meta-analysis

Purpose: Use this workflow when studies already provide an effect size with its variance or standard error.

data <- read_example("metaEffectSizes.csv")

if (has_pkg("metafor")) {
  fit <- metafor::rma(
    yi = effect_size,
    vi = variance,
    mods = ~ moderator,
    data = data,
    method = "REML"
  )

  summary(fit)
  metafor::forest(fit, slab = data$study_label)
  metafor::funnel(fit)
  metafor::regtest(fit, model = "rma")
} else {
  message("Install metafor to run this meta-analysis.")
}

## 
## Regression Test for Funnel Plot Asymmetry
## 
## Model:     mixed-effects meta-regression model
## Predictor: standard error
## 
## Test for Funnel Plot Asymmetry: z = -0.2219, p = 0.8244

Fixed-effect meta-analysis

Purpose: Use a fixed-effect meta-analysis when the goal is to estimate one common effect assumed to be shared by all included studies.

data <- read_example("metaEffectSizes.csv")

if (has_pkg("metafor")) {
  fit <- metafor::rma(
    yi = effect_size,
    sei = se,
    data = data,
    method = "FE"
  )

  summary(fit)
  metafor::forest(fit, slab = data$study_label)
} else {
  message("Install metafor to run this meta-analysis.")
}

Random-effects meta-analysis

Purpose: Use a random-effects meta-analysis when true effects are expected to vary across studies and the goal is to estimate the mean of that distribution.

data <- read_example("metaEffectSizes.csv")

if (has_pkg("metafor")) {
  fit <- metafor::rma(
    yi = effect_size,
    sei = se,
    data = data,
    method = "REML"
  )

  summary(fit)
  predict(fit)
  metafor::forest(fit, slab = data$study_label)
} else {
  message("Install metafor to run this meta-analysis.")
}

Subgroup meta-analysis

Purpose: Use subgroup meta-analysis to examine whether pooled effects differ across categorical study characteristics.

data <- read_example("metaEffectSizes.csv")

if (has_pkg("metafor")) {
  fit <- metafor::rma(
    yi = effect_size,
    sei = se,
    mods = ~ moderator,
    data = data,
    method = "REML"
  )

  summary(fit)

  by(
    data,
    data$moderator,
    function(d) summary(metafor::rma(yi = effect_size, sei = se, data = d, method = "REML"))
  )
} else {
  message("Install metafor to run this meta-analysis.")
}

## data$moderator: psychology
## 
## Random-Effects Model (k = 5; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
##   1.5074   -3.0148    0.9852   -0.2422   12.9852   
## 
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0242)
## tau (square root of estimated tau^2 value):      0
## I^2 (total heterogeneity / total variability):   0.00%
## H^2 (total variability / sampling variability):  1.00
## 
## Test for Heterogeneity:
## Q(df = 4) = 2.9246, p-val = 0.5705
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub      
##  -0.5400  0.0830  -6.5050  <.0001  -0.7028  -0.3773  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------ 
## data$moderator: sport
## 
## Random-Effects Model (k = 5; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
##  -3.6788    7.3575   11.3575   10.1301   23.3575   
## 
## tau^2 (estimated amount of total heterogeneity): 0.3193 (SE = 0.2598)
## tau (square root of estimated tau^2 value):      0.5651
## I^2 (total heterogeneity / total variability):   87.01%
## H^2 (total variability / sampling variability):  7.70
## 
## Test for Heterogeneity:
## Q(df = 4) = 30.4576, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub    
##   0.1261  0.2711  0.4651  0.6418  -0.4053  0.6574    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Meta-regression

Purpose: Use meta-regression to examine whether study-level moderators explain variation in effect sizes across studies.

data <- read_example("metaEffectSizes.csv")

if (has_pkg("metafor")) {
  fit <- metafor::rma(
    yi = effect_size,
    sei = se,
    mods = ~ moderator,
    data = data,
    method = "REML"
  )

  summary(fit)
} else {
  message("Install metafor to run meta-regression.")
}

## 
## Mixed-Effects Model (k = 10; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
##  -4.7349    9.4698   15.4698   15.7082   21.4698   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.1431 (SE = 0.0925)
## tau (square root of estimated tau^2 value):             0.3783
## I^2 (residual heterogeneity / unaccounted variability): 77.97%
## H^2 (unaccounted variability / sampling variability):   4.54
## R^2 (amount of heterogeneity accounted for):            44.47%
## 
## Test for Residual Heterogeneity:
## QE(df = 8) = 33.3822, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 6.4230, p-val = 0.0113
## 
## Model Results:
## 
##                 estimate      se     zval    pval    ci.lb    ci.ub     
## intrcpt          -0.5659  0.1895  -2.9859  0.0028  -0.9373  -0.1944  ** 
## moderatorsport    0.6901  0.2723   2.5344  0.0113   0.1564   1.2238   * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Publication-bias diagnostics

Purpose: Use publication-bias diagnostics to inspect funnel plot asymmetry and possible small-study effects after fitting a meta-analysis model.

data <- read_example("metaEffectSizes.csv")

if (has_pkg("metafor")) {
  fit <- metafor::rma(
    yi = effect_size,
    sei = se,
    data = data,
    method = "REML"
  )

  metafor::funnel(fit)
  metafor::regtest(fit, model = "rma")
  metafor::trimfill(fit)
} else {
  message("Install metafor to run publication-bias diagnostics.")
}

## 
## Estimated number of missing studies on the left side: 0 (SE = 2.0440)
## 
## Random-Effects Model (k = 10; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.2577 (SE = 0.1414)
## tau (square root of estimated tau^2 value):      0.5076
## I^2 (total heterogeneity / total variability):   86.48%
## H^2 (total variability / sampling variability):  7.40
## 
## Test for Heterogeneity:
## Q(df = 9) = 59.1993, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb   ci.ub    
##  -0.2279  0.1732  -1.3159  0.1882  -0.5674  0.1115    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Proportion meta-analysis

Purpose: Use proportion meta-analysis to synthesize event proportions or prevalence estimates across studies.

data <- read_example("metaProportion.csv")

if (has_pkg("metafor")) {
  es <- metafor::escalc(
    measure = "PLO",
    xi = event_frequency,
    ni = total_sample_size,
    data = data
  )

  fit <- metafor::rma(yi, vi, mods = ~ moderator_code, data = es, method = "REML")
  summary(fit)
  predict(fit, transf = metafor::transf.ilogit)
  metafor::forest(fit, slab = data$study_label)
} else {
  message("Install metafor to run this meta-analysis.")
}

Decision tree regression

Purpose: Use decision tree regression to predict a continuous outcome with transparent if-then decision rules.

data <- read_example("decisionTreeRegression.csv")
model_data <- subset(data, select = -client)

if (has_pkg("rpart")) {
  set.seed(1)
  train_id <- sample(seq_len(nrow(model_data)), size = floor(0.70 * nrow(model_data)))

  fit <- rpart::rpart(
    outcome_symptom_score ~ motivation_score + baseline_distress_score + sessions_attended,
    data = model_data[train_id, ],
    method = "anova",
    control = rpart::rpart.control(cp = .01)
  )

  pred <- predict(fit, newdata = model_data[-train_id, ])
  actual <- model_data$outcome_symptom_score[-train_id]
  c(
    RMSE = sqrt(mean((pred - actual)^2)),
    MAE = mean(abs(pred - actual)),
    R2 = cor(pred, actual)^2
  )
  plot(fit)
  text(fit)
} else {
  message("Install rpart to run decision tree regression.")
}

Random forest regression

Purpose: Use random forest regression to improve continuous-outcome prediction by averaging many decision trees and estimating variable importance.

data <- read_example("randomForestRegression.csv")
model_data <- subset(data, select = -athlete)

if (has_pkg("randomForest")) {
  set.seed(1)
  train_id <- sample(seq_len(nrow(model_data)), size = floor(0.70 * nrow(model_data)))

  fit <- randomForest::randomForest(
    race_time_min ~ .,
    data = model_data[train_id, ],
    ntree = 500,
    importance = TRUE
  )

  pred <- predict(fit, newdata = model_data[-train_id, ])
  actual <- model_data$race_time_min[-train_id]
  c(
    RMSE = sqrt(mean((pred - actual)^2)),
    MAE = mean(abs(pred - actual)),
    R2 = cor(pred, actual)^2
  )
  randomForest::importance(fit)
} else {
  message("Install randomForest to run random forest regression.")
}

##                           %IncMSE IncNodePurity
## training_volume_km_week 13.841384     293.01803
## vo2max_ml_kg_min        12.687973     247.91246
## tempo_pace_min_km       11.619106     230.47331
## long_run_km             11.077403     231.25752
## recovery_score          11.268020     250.28043
## sleep_hours              9.712273     170.88730
## age_years                7.212518      26.74854
## injury_days_last_month   6.256555      60.74022

k-nearest neighbors regression

Purpose: Use k-nearest neighbors regression to predict a continuous outcome from the outcomes of similar cases in the predictor space.

data <- read_example("knnRegression.csv")
model_data <- subset(data, select = -participant)

if (has_pkg("caret")) {
  set.seed(1)
  ctrl <- caret::trainControl(method = "cv", number = 10)

  fit <- caret::train(
    wellbeing_score ~ sleep_hours + stress_score + physical_activity_min_week + social_support_score,
    data = model_data,
    method = "knn",
    trControl = ctrl,
    preProcess = c("center", "scale"),
    tuneLength = 10
  )

  fit
} else {
  message("Install caret to run k-nearest neighbors regression.")
}

## k-Nearest Neighbors 
## 
## 60 samples
##  4 predictor
## 
## Pre-processing: centered (4), scaled (4) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 54, 52, 55, 54, 54, 56, ... 
## Resampling results across tuning parameters:
## 
##   k   RMSE       Rsquared   MAE      
##    5  0.8780874  0.9953210  0.6747143
##    7  0.9991272  0.9947135  0.7866497
##    9  1.2284369  0.9930314  0.9449735
##   11  1.3550122  0.9917154  0.9799919
##   13  1.5656150  0.9897858  1.1271978
##   15  1.8735785  0.9853155  1.3204365
##   17  2.1412111  0.9812369  1.4230182
##   19  2.4874152  0.9787332  1.6676629
##   21  2.8636398  0.9726694  1.9309240
##   23  3.2583268  0.9672477  2.1707557
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 5.

Decision tree classifier

Purpose: Use a decision tree classifier to predict categorical group membership with interpretable decision rules.

data <- read_example("decisionTreeClassifier.csv")
data$class <- factor(data$class)
model_data <- subset(data, select = -participant)

if (has_pkg("rpart")) {
  set.seed(1)
  train_id <- sample(seq_len(nrow(model_data)), size = floor(0.70 * nrow(model_data)))

  fit <- rpart::rpart(
    class ~ .,
    data = model_data[train_id, ],
    method = "class"
  )

  pred <- predict(fit, newdata = model_data[-train_id, ], type = "class")
  table(Predicted = pred, Observed = model_data$class[-train_id])
  mean(pred == model_data$class[-train_id])

  # plot(fit)
  # text(fit)
} else {
  message("Install rpart to run decision tree classification.")
}

## [1] 0.4444444

Random forest classifier

Purpose: Use a random forest classifier to improve categorical prediction by combining many trees and evaluating variable importance.

data <- read_example("randomForestClassifier.csv")
data$class <- factor(data$class)
model_data <- subset(data, select = -athlete)

if (has_pkg("randomForest")) {
  set.seed(1)
  train_id <- sample(seq_len(nrow(model_data)), size = floor(0.70 * nrow(model_data)))

  fit <- randomForest::randomForest(
    class ~ .,
    data = model_data[train_id, ],
    ntree = 500,
    importance = TRUE
  )

  pred <- predict(fit, newdata = model_data[-train_id, ])
  table(Predicted = pred, Observed = model_data$class[-train_id])
  mean(pred == model_data$class[-train_id])
  randomForest::importance(fit)
} else {
  message("Install randomForest to run random forest classification.")
}

##                       high       low  moderate MeanDecreaseAccuracy
## training_load_au 12.321786 10.713224  9.102714            14.019030
## recovery_score   12.485727 11.785023  9.094083            15.010969
## sleep_hours      12.345697 10.550427  8.753778            13.743853
## previous_injury   3.567065  3.940552 -3.461392             3.238400
## soreness_score   12.862177 12.492907  9.138869            15.110068
## age_years         8.292441  7.039137  2.418129             9.188735
##                  MeanDecreaseGini
## training_load_au        2.5642061
## recovery_score          2.8597106
## sleep_hours             2.5171404
## previous_injury         0.3698107
## soreness_score          3.0796034
## age_years               1.2720049

k-nearest neighbors classifier

Purpose: Use k-nearest neighbors classification to assign cases to categories based on the classes of similar observations.

data <- read_example("knnClassifier.csv")
data$class <- factor(data$class)
model_data <- subset(data, select = -participant)

if (has_pkg("caret")) {
  set.seed(1)
  ctrl <- caret::trainControl(method = "cv", number = 10)

  fit <- caret::train(
    class ~ .,
    data = model_data,
    method = "knn",
    trControl = ctrl,
    preProcess = c("center", "scale"),
    tuneLength = 10
  )

  fit
} else {
  message("Install caret to run k-nearest neighbors classification.")
}

## k-Nearest Neighbors 
## 
## 30 samples
##  5 predictor
##  3 classes: 'high_wellbeing', 'low_wellbeing', 'moderate_wellbeing' 
## 
## Pre-processing: centered (5), scaled (5) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 27, 27, 27, 27, 27, 27, ... 
## Resampling results across tuning parameters:
## 
##   k   Accuracy   Kappa
##    5  1.0000000  1.00 
##    7  1.0000000  1.00 
##    9  1.0000000  1.00 
##   11  1.0000000  1.00 
##   13  1.0000000  1.00 
##   15  1.0000000  1.00 
##   17  1.0000000  1.00 
##   19  0.6333333  0.45 
##   21  0.6333333  0.45 
##   23  0.4666667  0.20 
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 17.

Naive Bayes classifier

Purpose: Use a Naive Bayes classifier as a fast probabilistic model for categorical prediction, especially as a baseline model.

data <- read_example("naiveBayes.csv")
data$class <- factor(data$class)
model_data <- subset(data, select = -student)

if (has_pkg("e1071")) {
  set.seed(1)
  train_id <- sample(seq_len(nrow(model_data)), size = floor(0.70 * nrow(model_data)))

  fit <- e1071::naiveBayes(class ~ ., data = model_data[train_id, ])
  pred <- predict(fit, newdata = model_data[-train_id, ])

  table(Predicted = pred, Observed = model_data$class[-train_id])
  mean(pred == model_data$class[-train_id])
} else {
  message("Install e1071 to run naive Bayes classification.")
}

## [1] 1

Principal component analysis

Purpose: Use PCA to reduce several correlated continuous variables into fewer components that summarize the main sources of variance.

data <- read_example("principalComponentAnalysis.csv")
variables <- c(
  "attention_accuracy",
  "working_memory",
  "processing_speed",
  "inhibition_score",
  "cognitive_flexibility",
  "verbal_memory",
  "spatial_reasoning"
)

vars <- data[, variables]
fit <- prcomp(vars, center = TRUE, scale. = TRUE)

summary(fit)

## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5     PC6       PC7
## Standard deviation     2.6441 0.06736 0.05037 0.03283 0.01997 0.01451 3.081e-17
## Proportion of Variance 0.9988 0.00065 0.00036 0.00015 0.00006 0.00003 0.000e+00
## Cumulative Proportion  0.9988 0.99940 0.99976 0.99991 0.99997 1.00000 1.000e+00

fit$rotation

##                              PC1         PC2          PC3          PC4
## attention_accuracy    -0.3780990  0.16169560  0.319999090 -0.333217455
## working_memory        -0.3777388 -0.70127074 -0.008562141 -0.319983152
## processing_speed      -0.3780990  0.16169560  0.319999090 -0.333217455
## inhibition_score      -0.3779550 -0.09117910  0.493009067  0.772979685
## cognitive_flexibility -0.3781527  0.01351191 -0.173274585 -0.066199878
## verbal_memory         -0.3777806  0.64224748 -0.289534638  0.008492108
## spatial_reasoning     -0.3779260 -0.18695021 -0.661978493  0.271276436
##                               PC5          PC6           PC7
## attention_accuracy     0.30522047  0.155893669 -7.071068e-01
## working_memory        -0.49847680  0.120793591  1.830799e-16
## processing_speed       0.30522047  0.155893669  7.071068e-01
## inhibition_score      -0.09091677 -0.003872043  2.082403e-16
## cognitive_flexibility  0.09445931 -0.901936194  2.016684e-16
## verbal_memory         -0.57836888  0.162440800  1.785805e-16
## spatial_reasoning      0.46206366  0.311307457 -4.847109e-16

biplot(fit, scale = 0)