Solutions: Exercises 3–4

Exercise 3 — t-tests

library(tidyverse)
library(effectsize)

employees <- read_csv("data/employee_survey.csv", show_col_types = FALSE)
experiment <- read_csv("data/training_experiment.csv", show_col_types = FALSE)

cor.test(employees$workload, employees$burnout, use = "complete.obs")

    Pearson's product-moment correlation

data:  employees$workload and employees$burnout
t = 19.429, df = 478, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6111010 0.7114717
sample estimates:
    cor 
0.66427 
experiment |>
  group_by(group) |>
  summarise(
    n = n(),
    mean_post = mean(post_score),
    sd_post = sd(post_score),
    .groups = "drop"
  )
# A tibble: 2 × 4
  group        n mean_post sd_post
  <chr>    <int>     <dbl>   <dbl>
1 Control    115      65.1    11.3
2 Training   105      68.5    10.1
two_sided <- t.test(
  post_score ~ group,
  data = experiment,
  alternative = "two.sided"
)
two_sided

    Welch Two Sample t-test

data:  post_score by group
t = -2.3536, df = 217.85, p-value = 0.01948
alternative hypothesis: true difference in means between group Control and group Training is not equal to 0
95 percent confidence interval:
 -6.2358887 -0.5517717
sample estimates:
 mean in group Control mean in group Training 
              65.10522               68.49905 
training_post <- experiment |>
  filter(group == "Training") |>
  pull(post_score)

control_post <- experiment |>
  filter(group == "Control") |>
  pull(post_score)

one_sided <- t.test(
  training_post,
  control_post,
  alternative = "greater"
)
one_sided

    Welch Two Sample t-test

data:  training_post and control_post
t = 2.3536, df = 217.85, p-value = 0.009742
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 1.01182     Inf
sample estimates:
mean of x mean of y 
 68.49905  65.10522 
cohens_d(post_score ~ group, data = experiment)
Cohen's d |         95% CI
--------------------------
-0.32     | [-0.58, -0.05]

- Estimated using pooled SD.
training_only <- experiment |>
  filter(group == "Training")

t.test(
  training_only$post_score,
  training_only$pre_score,
  paired = TRUE,
  alternative = "two.sided"
)

    Paired t-test

data:  training_only$post_score and training_only$pre_score
t = 19.965, df = 104, p-value < 2.2e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 6.501155 7.935035
sample estimates:
mean difference 
       7.218095

The directional alternative must be justified before inspecting the relevant outcome. The one-sided p-value is not a legitimate “upgrade” chosen because the two-sided result is inconvenient.

Exercise 4 — measurement

library(tidyverse)
library(psych)
library(lavaan)

employees <- read_csv("data/employee_survey.csv", show_col_types = FALSE) |>
  mutate(
    leadership_4 = 6 - leadership_4_reverse,
    engagement_4 = 6 - engagement_4_reverse
  )

leadership_items <- employees |>
  select(leadership_1, leadership_2, leadership_3, leadership_4)

alpha(leadership_items)

Reliability analysis   
Call: alpha(x = leadership_items)

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean  sd median_r
      0.87      0.87    0.84      0.63 6.8 0.0095    3 1.2     0.64

    95% confidence boundaries 
         lower alpha upper
Feldt     0.85  0.87  0.89
Duhachek  0.85  0.87  0.89

 Reliability if an item is dropped:
             raw_alpha std.alpha G6(smc) average_r S/N alpha se   var.r med.r
leadership_1      0.82      0.82    0.76      0.60 4.6    0.014 0.00152  0.63
leadership_2      0.84      0.84    0.78      0.65 5.5    0.012 0.00032  0.65
leadership_3      0.85      0.85    0.79      0.65 5.6    0.012 0.00034  0.66
leadership_4      0.83      0.83    0.77      0.62 4.9    0.013 0.00290  0.65

 Item statistics 
               n raw.r std.r r.cor r.drop mean  sd
leadership_1 480  0.87  0.87  0.82   0.76  3.0 1.3
leadership_2 472  0.84  0.84  0.76   0.71  3.0 1.4
leadership_3 480  0.83  0.83  0.75   0.70  3.0 1.4
leadership_4 472  0.86  0.86  0.79   0.74  2.9 1.4

Non missing response frequency for each item
                1    2    3    4    5 miss
leadership_1 0.17 0.22 0.23 0.20 0.18 0.00
leadership_2 0.20 0.21 0.21 0.19 0.19 0.02
leadership_3 0.19 0.22 0.21 0.19 0.19 0.00
leadership_4 0.19 0.22 0.24 0.17 0.18 0.02
item_names <- c(
  "leadership_1", "leadership_2", "leadership_3", "leadership_4",
  "engagement_1", "engagement_2", "engagement_3", "engagement_4",
  "role_clarity_1", "role_clarity_2", "role_clarity_3", "role_clarity_4"
)

measurement_data <- employees |>
  select(all_of(item_names)) |>
  drop_na()

set.seed(2026)
explore_rows <- sample(
  seq_len(nrow(measurement_data)),
  floor(nrow(measurement_data) / 2)
)

explore_data <- measurement_data[explore_rows, ]
confirm_data <- measurement_data[-explore_rows, ]

pca_fit <- prcomp(explore_data, center = TRUE, scale. = TRUE)
summary(pca_fit)
Importance of components:
                         PC1    PC2    PC3     PC4     PC5     PC6     PC7
Standard deviation     2.388 1.2532 1.1163 0.74414 0.70257 0.65685 0.65367
Proportion of Variance 0.475 0.1309 0.1038 0.04615 0.04113 0.03595 0.03561
Cumulative Proportion  0.475 0.6059 0.7098 0.75590 0.79704 0.83299 0.86860
                           PC8     PC9    PC10    PC11    PC12
Standard deviation     0.60261 0.59243 0.57230 0.52889 0.50540
Proportion of Variance 0.03026 0.02925 0.02729 0.02331 0.02129
Cumulative Proportion  0.89886 0.92811 0.95540 0.97871 1.00000
round(pca_fit$rotation[, 1:3], 2)
                 PC1   PC2   PC3
leadership_1   -0.32  0.24 -0.30
leadership_2   -0.29  0.25 -0.34
leadership_3   -0.29  0.21 -0.32
leadership_4   -0.31  0.28 -0.28
engagement_1   -0.29  0.13  0.37
engagement_2   -0.30  0.09  0.36
engagement_3   -0.30  0.09  0.33
engagement_4   -0.28  0.19  0.47
role_clarity_1 -0.26 -0.43 -0.12
role_clarity_2 -0.27 -0.42 -0.05
role_clarity_3 -0.27 -0.42 -0.03
role_clarity_4 -0.27 -0.38 -0.05
fa.parallel(explore_data, fa = "fa", plot = FALSE)
Parallel analysis suggests that the number of factors =  3  and the number of components =  NA 
efa_fit <- fa(
  explore_data,
  nfactors = 3,
  rotate = "oblimin",
  fm = "minres"
)
print(efa_fit$loadings, cutoff = .30)

Loadings:
               MR1    MR2    MR3   
leadership_1    0.875              
leadership_2    0.813              
leadership_3    0.726              
leadership_4    0.820              
engagement_1                  0.700
engagement_2                  0.673
engagement_3                  0.637
engagement_4                  0.900
role_clarity_1         0.773       
role_clarity_2         0.775       
role_clarity_3         0.788       
role_clarity_4         0.680       

                 MR1   MR2   MR3
SS loadings    2.639 2.319 2.171
Proportion Var 0.220 0.193 0.181
Cumulative Var 0.220 0.413 0.594
cfa_model <- '
  Leadership =~ leadership_1 + leadership_2 + leadership_3 + leadership_4
  Engagement =~ engagement_1 + engagement_2 + engagement_3 + engagement_4
  RoleClarity =~ role_clarity_1 + role_clarity_2 + role_clarity_3 + role_clarity_4
'

cfa_fit <- cfa(
  cfa_model,
  data = confirm_data,
  ordered = item_names,
  estimator = "WLSMV",
  std.lv = TRUE
)

fitMeasures(cfa_fit, c("cfi", "tli", "rmsea", "srmr"))
  cfi   tli rmsea  srmr 
1.000 1.008 0.000 0.031

A concise distinction:

  • PCA summarizes observed variance with weighted components.
  • EFA explores a plausible latent factor structure.
  • CFA evaluates a measurement structure specified before inspecting the confirmatory fit.