Two-way mixed repeated-measures ANOVA

afex::aov_ez() で mixed repeated-measures ANOVA を行い、その後に pairwise比較と各比較の効果量を出す基本例の記録。

# 初回のみ
# install.packages(c("readr", "dplyr", "afex", "emmeans"))

library(readr)
library(dplyr)
library(afex)
library(emmeans)

# データ読み込み
dat <- read_csv("data.csv")

# 列の型を指定
dat <- dat %>%
  mutate(
    subject   = factor(subject),
    Dose      = factor(Dose),
    Timepoint = factor(Timepoint, levels = c("T1", "T2", "T3"))
  )

# Mixed repeated measures ANOVA
# Dose: 群間要因
# Timepoint: 反復測定要因
fit_aov <- aov_ez(
  id = "subject",
  dv = "Value",
  between = "Dose",
  within = "Timepoint",
  data = dat,
  type = 3
)

# ANOVA表
fit_aov

(1). interactionが有意な場合
A. 各timepointで treatment群間比較
「どの時点でcontrolとdrugが違うか」を見る場合です。

emm_treat_by_time <- emmeans(fit, ~ treatment | timepoint)

pair_treat_by_time <- pairs(
  emm_treat_by_time,
  adjust = "bonferroni"
)

pair_treat_by_time

B. 各treatment群内で timepoint間比較
「各群でどのtimepoint間が違うか」を見る場合です。

emm_time_by_treat <- emmeans(fit, ~ timepoint | treatment)

pair_time_by_treat <- pairs(
  emm_time_by_treat,
  adjust = "bonferroni"
)

pair_time_by_treat

効果量の計算

library(dplyr)
library(emmeans)

effect_size_result <- bind_rows(
  lapply(levels(dat$Time), function(tt) {

    dat_t <- dat %>% filter(Time == tt)

    # このモデルは、各時点の効果量の算出に使用
    fit_t <- lm(Intensity ~ Treatment, data = dat_t)
    emm_t <- emmeans(fit_t, ~ Treatment)

    es_t <- eff_size(
      emm_t,
      sigma = sigma(fit_t),
      edf = df.residual(fit_t),
      method = "pairwise"
    )

    as.data.frame(es_t) %>%
      mutate(Time = tt, .before = 1)
  })
)

effect_size_result

# 初回のみ
# install.packages("effectsize")

library(dplyr)
library(tidyr)

wide <- dat %>%
  select(Subject, Treatment, Time, Intensity) %>%
  pivot_wider(
    id_cols = c(Subject, Treatment),
    names_from = Time,
    values_from = Intensity,
    names_prefix = "Time_"
  )

time_pairs <- combn(
  levels(dat$Time),
  2,
  simplify = FALSE
)

effect_size_time <- bind_rows(
  lapply(levels(dat$Treatment), function(tr) {

    group_data <- wide %>% filter(Treatment == tr)

    bind_rows(lapply(time_pairs, function(tp) {

      x <- group_data[[paste0("Time_", tp[1])]]
      y <- group_data[[paste0("Time_", tp[2])]]

      # 同じ個体の両時点の値がそろっているペアを使用
      keep <- complete.cases(x, y)
      x <- x[keep]
      y <- y[keep]

      es <- effectsize::cohens_d(
        x,
        y,
        paired = TRUE,
        ci = 0.95,
        alternative = "two.sided"
      )

      tibble(
        Treatment = tr,
        contrast = paste(tp[1], "-", tp[2]),
        n = length(x),
        mean_difference = mean(x - y),
        dz = es$Cohens_d,
        CI_low = es$CI_low,
        CI_high = es$CI_high
      )
    }))
  })
)

effect_size_time

p値は反復測定ANOVAモデルに基づく比較、効果量は実測値の対応差に基づく(d_z)という組み合わせです。論文では単に「Cohen’s d」ではなく、“Cohen’s (d_z), calculated using the standard deviation of paired differences”**と定義を明示する

(2). interactionなし、timepoint主効果のみ有意な場合
この場合は、群をまとめてtimepoint間比較します。

emm_time <- emmeans(fit, ~ timepoint)

pair_time <- pairs(
  emm_time,
  adjust = "bonferroni"
)

pair_time

効果量

# emm_time <- emmeans(fit_aov, ~ Time)
# pairs(emm_time, adjust = "holm")

library(dplyr)
library(tidyr)
library(emmeans)

# datにはSubject, Treatment, Time, Intensityがある想定
wide <- dat %>%
  select(Subject, Treatment, Time, Intensity) %>%
  pivot_wider(
    id_cols = c(Subject, Treatment),
    names_from = Time,
    values_from = Intensity,
    names_prefix = "Time_"
  )

time_pairs <- combn(
  levels(dat$Time),
  2,
  simplify = FALSE
)

effect_size_time <- bind_rows(
  lapply(time_pairs, function(tp) {

    pair_dat <- wide %>%
      mutate(
        difference =
          .data[[paste0("Time_", tp[1])]] -
          .data[[paste0("Time_", tp[2])]]
      ) %>%
      filter(!is.na(difference))

    # 各個体の時点間差をTreatmentで調整
    fit_diff <- lm(difference ~ Treatment, data = pair_dat)

    # Treatmentを等しい重みで平均した時点間差
    emm_diff <- emmeans(fit_diff, ~ 1, weights = "equal")

    # 推定された差を、その残差SDで標準化
    es <- eff_size(
      emm_diff,
      sigma = sigma(fit_diff),
      edf = df.residual(fit_diff),
      method = "identity"
    )

    es_table <- as.data.frame(
      summary(es, infer = c(TRUE, FALSE), level = 0.95)
    )

    tibble(
      contrast = paste(tp[1], "-", tp[2]),
      n = nrow(pair_dat),
      dz_adjusted = es_table$effect.size,
      CI_low = es_table$lower.CL,
      CI_high = es_table$upper.CL
    )
  })
)

effect_size_time

効果量は、各時点ペアの個体内差を、Treatmentで調整した差のSDで標準化した(d_z)型の指標として、以下で95%信頼区間とともに算出。各個体の差にlm(difference ~ Treatment)を当て、Treatmentを等しい重みで平均するため、上記のemmeans()と比較の対象がそろいます。