History Tests

Statistical Significance

One-sided Cochran-Armitage trend test:

library(DescTools)

data_table <- matrix(
  c(
    63, 57, 57, 51, # Row of incorrect answers
    13, 19, 19, 25 # Row of correct answers
  ),
  nrow = 2,
  byrow = TRUE
)

colnames(data_table) <- c("Difficult_1", "Difficult_2", "Difficult_3", "Difficult_4")
rownames(data_table) <- c("False", "True")

CochranArmitageTest(data_table, alternative = "one.sided")

    Cochran-Armitage test for trend

data:  data_table
Z = -2.1325, dim = 4, p-value = 0.01648
alternative hypothesis: one.sided

Effect Size

Graph of probabilities predicted by logistic regression (a visual way to estimate the effect size):

library(dplyr)
library(ggplot2)
library(marginaleffects)

# Data
data <- data.frame(
  detail_level = c(rep(1, 76), rep(2, 76), rep(3, 76), rep(4, 76)),
  correct = c(
    rep(1, 13), rep(0, 63), # Level 1
    rep(1, 19), rep(0, 57), # Level 2
    rep(1, 19), rep(0, 57), # Level 3
    rep(1, 25), rep(0, 51) # Level 4
  )
)

# Build logistic regression
model <- glm(correct ~ detail_level, data = data, family = binomial)

# Predict probabilities at each detail level
all_levels <- predictions(
  model,
  newdata = datagrid(detail_level = 1:4)
)

# Build graph
ggplot(all_levels, aes(x = factor(detail_level), y = estimate)) +
  # Dashed line for random guessing
  geom_hline(yintercept = 0.25, linetype = "dashed", color = "gray50") +
  # Confidence intervals
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1, color = "gray30", linewidth = 0.8) +
  # Points of predicted probabilities
  geom_point(size = 4, color = "gray30") +
  # Add a thin line connecting points to emphasize the trend
  geom_line(aes(group = 1), linetype = "dotted", color = "gray30") +
  # Axis settings
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    breaks = seq(0.1, 0.5, by = 0.05),
    limits = c(0.1, 0.45)
  ) +
  theme_minimal() +
  labs(
    x = "Answer option detail level",
    y = "Probability of correct answer",
  ) +
  theme(text = element_text(size = 16, family = "serif"))

print(all_levels)

 detail_level Estimate Pr(>|z|)    S 2.5 % 97.5 %
            1    0.182   <0.001 31.1 0.122  0.263
            2    0.223   <0.001 49.8 0.175  0.280
            3    0.271   <0.001 39.5 0.220  0.328
            4    0.324   <0.001 11.1 0.241  0.420

Type: invlink(link)

Practically interesting effect size (the strategy of choosing the most detailed answer option allows guessing with a probability 0.0741 higher):

library(marginaleffects)

# Data
data <- data.frame(
  detail_level = c(rep(1, 76), rep(2, 76), rep(3, 76), rep(4, 76)),
  correct = c(
    rep(1, 13), rep(0, 63), # Level 1
    rep(1, 19), rep(0, 57), # Level 2
    rep(1, 19), rep(0, 57), # Level 3
    rep(1, 25), rep(0, 51) # Level 4
  )
)

# Build logistic regression
model <- glm(correct ~ detail_level, data = data, family = binomial)

# Calculate the difference in probabilities
results <- predictions(
  model,
  newdata = datagrid(detail_level = 4),
  hypothesis = "b1 - 0.25 = 0",
  type = "response"
)

print(results)

 Hypothesis Estimate Std. Error    z Pr(>|z|)   S  2.5 % 97.5 %
  b1-0.25=0   0.0741      0.046 1.61    0.107 3.2 -0.016  0.164

Type: response

Power Analysis

Post-hoc Power Analysis

Post-hoc power analysis for the Cochran-Armitage test:

library(multiCA)

power.CA.test(
  N = 304, # Sample size (number of answer options)
  n.prop = c(0.25, 0.25, 0.25, 0.25), # Proportion of group sizes at each level
  pvec = c(0.1711, 0.25, 0.25, 0.3289), # Actual proportions of correct answers at each level
  sig.level = 0.05, # Alpha significance level
  alternative = "greater" # Direction of testing
)

     Cochran-Armitage trend test 

              n = 304
         n.prop = 0.25, 0.25, 0.25, 0.25
              p = 0.1711, 0.2500, 0.2500, 0.3289
    alternative = greater
      sig.level = 0.05
          power = 0.6892607

A-priori Power Analysis

A-priori power analysis for the Cochran-Armitage test:

library(multiCA)

power.CA.test(
  n.prop = c(0.25, 0.25, 0.25, 0.25), # Proportion of group sizes at each level
  pvec = c(0.1711, 0.25, 0.25, 0.3289), # Actual proportions of correct answers at each level
  sig.level = 0.05, # Alpha significance level
  alternative = "greater", # Direction of testing
  power = 0.8 # Required power
)

     Cochran-Armitage trend test 

              n = 409.6055
         n.prop = 0.25, 0.25, 0.25, 0.25
              p = 0.1711, 0.2500, 0.2500, 0.3289
    alternative = greater
      sig.level = 0.05
          power = 0.8

A-priori power analysis for a practically interesting effect size using the Monte Carlo method. Statistical power calculation for a sample size of 172 questions (688 answer options):

library(marginaleffects)

# Function to calculate statistical power by sample size n_per_level
simulate_power_exact <- function(n_per_level, b0, b1, n_sim, alpha) {
  success_count <- 0

  for (i in 1:n_sim) {
    # 1. Generate data
    sim_detail_level <- rep(1:4, each = n_per_level)
    sim_log_odds <- b0 + b1 * sim_detail_level
    sim_probs <- 1 / (1 + exp(-sim_log_odds))
    sim_correct <- rbinom(length(sim_detail_level), size = 1, prob = sim_probs)

    sim_data <- data.frame(detail_level = sim_detail_level, correct = sim_correct)

    # 2. Train model
    fit <- glm(correct ~ detail_level, data = sim_data, family = binomial)

    # 3. Use marginaleffects to check effect significance
    results <- predictions(
      fit,
      newdata = datagrid(detail_level = 4),
      hypothesis = "b1 - 0.25 = 0",
      type = "response"
    )

    # 4. Extract one-sided p-value
    # Divide by 2 if the effect is in the expected direction (Estimate > 0)
    if (results$estimate > 0) {
      p_val_one_sided <- results$p.value / 2
    } else {
      p_val_one_sided <- 1 - (results$p.value / 2)
    }

    if (p_val_one_sided < alpha) {
      success_count <- success_count + 1
    }
  }

  return(success_count / n_sim)
}

# Calculate statistical power
set.seed(42)
estimated_power <- simulate_power_exact(
  n_per_level = 172, # Sample size per level (number of questions)
  b0 = coef(model)[1], # Constant (intercept) in regression
  b1 = coef(model)[2], # Coefficient for the detail_level variable
  n_sim = 1000, # Number of simulations
  alpha = 0.05 # Significance level α
)

print(estimated_power)

[1] 0.806