Тесты по истории

Статистическая значимость

Односторонний тест Кокрана-Армитажа на тренд:

library(DescTools)

data_table <- matrix(
  c(
    63, 57, 57, 51, # Строка неправильных ответов
    13, 19, 19, 25 # Строка правильных ответов
  ),
  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

Размер эффекта

График вероятностей, предсказанных логистической регрессией (наглядный способ оценить размер эффекта):

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

# Данные
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), # Уровень 1
    rep(1, 19), rep(0, 57), # Уровень 2
    rep(1, 19), rep(0, 57), # Уровень 3
    rep(1, 25), rep(0, 51) # Уровень 4
  )
)

# Строим логистическую регрессия
model <- glm(correct ~ detail_level, data = data, family = binomial)

# Предсказываем вероятности на каждом уровне детализации
all_levels <- predictions(
  model,
  newdata = datagrid(detail_level = 1:4)
)

# Строим график
ggplot(all_levels, aes(x = factor(detail_level), y = estimate)) +
  # Пунктирная линия случайного угадывания
  geom_hline(yintercept = 0.25, linetype = "dashed", color = "gray50") +
  # Доверительные интервалы
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1, color = "gray30", linewidth = 0.8) +
  # Точки предсказанных вероятностей
  geom_point(size = 4, color = "gray30") +
  # Добавим тонкую линию, соединяющую точки, чтобы подчеркнуть тренд
  geom_line(aes(group = 1), linetype = "dotted", color = "gray30") +
  # Настройка осей
  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 = "Уровень детализации варианта ответа",
    y = "Вероятность ответить правильно",
  ) +
  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)

Практически интересный размера эффекта (стратегия выбора наиболее детализированного варианта ответа позволяет угадывать с вероятностью на 0.0741 выше):

library(marginaleffects)

# Данные
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), # Уровень 1
    rep(1, 19), rep(0, 57), # Уровень 2
    rep(1, 19), rep(0, 57), # Уровень 3
    rep(1, 25), rep(0, 51) # Уровень 4
  )
)

# Строим логистическую регрессию
model <- glm(correct ~ detail_level, data = data, family = binomial)

# Считаем разницу вероятностей
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

Анализ мощности

Апостериорные анализы мощности

Апостериорный анализ мощности для теста Кокрана-Армитажа:

library(multiCA)

power.CA.test(
  N = 304, # Размер выборки (количество вариантов ответа)
  n.prop = c(0.25, 0.25, 0.25, 0.25), # Пропорция размеров групп на каждом уровне
  pvec = c(0.1711, 0.25, 0.25, 0.3289), # Фактические пропорции правильных ответов на каждом уровне
  sig.level = 0.05, # Уровень значимости альфа
  alternative = "greater" # Направление проверки
)

     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

Априорные анализы мощности

Априорный анализ мощности для теста Кокрана-Армитажа:

library(multiCA)

power.CA.test(
  n.prop = c(0.25, 0.25, 0.25, 0.25), # Пропорция размеров групп на каждом уровне
  pvec = c(0.1711, 0.25, 0.25, 0.3289), # Фактические пропорции правильных ответов на каждом уровне
  sig.level = 0.05, # Уровень значимости альфа
  alternative = "greater", # Направление проверки
  power = 0.8 # Необходимая мощность
)

     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

Априорный анализ мощности для практически интересного размера эффекта методом Монте-Карло. Расчёт статистической мощности для выборки размером 172 вопроса (688 вариантов ответа):

library(marginaleffects)

# Функция расчёта статистической мощность по размеру выборки 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. Генерируем данные
    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. Обучаем модель
    fit <- glm(correct ~ detail_level, data = sim_data, family = binomial)

    # 3. Используем marginaleffects для проверки значимости эффекта
    results <- predictions(
      fit,
      newdata = datagrid(detail_level = 4),
      hypothesis = "b1 - 0.25 = 0",
      type = "response"
    )

    # 4. Извлекаем одностороннее p-значение
    # Делим на 2, если эффект направлен в ожидаемую сторону (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)
}

# Расчитываем статистическую мощность
set.seed(42)
estimated_power <- simulate_power_exact(
  n_per_level = 172, # Размер выборки на каждый уровень (количество вопросов)
  b0 = coef(model)[1], # Константа (intercept) в регрессии
  b1 = coef(model)[2], # Коэффициент для переменной detail_level
  n_sim = 100, # Количество симуляций
  alpha = 0.05 # Уровень значимости α
)

print(estimated_power)

[1] 0.84