【初心者向け】対応のある2群における3条件以上の統計解析の実践ガイド：Two-way repeated way ANOVA

# tidyverse パッケージを読み込む
library(tidyverse)

# データ作成
set.seed(123) # 再現性のためのシード

data <- tibble(
    ID = rep(1:10, each = 8), # 被験者ID
    Culture = rep(c("LowGlc", "HighGlc"), each = 4, times = 10), # 育成条件
    Condition = rep(c("Control", "Intervention1", "Intervention2", "Intervention3"), times = 2 * 10) # 治療条件
)

# 育成条件ごとにLDL値を調整
data <- data %>%
    mutate(
        LDL = case_when(
            Culture == "LowGlc" & Condition == "Control" ~ rnorm(n(), mean = 130, sd = 10),
            Culture == "LowGlc" & Condition == "Intervention1" ~ rnorm(n(), mean = 120, sd = 10),
            Culture == "LowGlc" & Condition == "Intervention2" ~ rnorm(n(), mean = 125, sd = 10),
            Culture == "LowGlc" & Condition == "Intervention3" ~ rnorm(n(), mean = 110, sd = 10),
            Culture == "HighGlc" & Condition == "Control" ~ rnorm(n(), mean = 140, sd = 10),
            Culture == "HighGlc" & Condition == "Intervention1" ~ rnorm(n(), mean = 130, sd = 10),
            Culture == "HighGlc" & Condition == "Intervention2" ~ rnorm(n(), mean = 135, sd = 10),
            Culture == "HighGlc" & Condition == "Intervention3" ~ rnorm(n(), mean = 125, sd = 10)
        )
    )

# 必要なパッケージを読み込む
library(tidyverse)
library(rstatix)

# Two-way repeated measures ANOVA
anova_results <- data %>%
  anova_test(
    dv = LDL, # 従属変数
    wid = ID, # 対応の単位
    within = c(Culture, Condition) # 独立変数
  )

# 結果を表示
anova_results

# 各条件（LowGlc, HighGlc）で治療条件のペアワイズt検定をまとめて実行
post_hoc_condition_by_culture <- data %>%
    group_by(Culture) %>% # Cultureごとにグループ化
    group_split() %>% # グループごとにデータを分割
    map(~ {
        test_result <- pairwise_t_test(
            .x, # グループごとのデータ
            formula = LDL ~ Condition, # 治療条件ごとの比較
            p.adjust.method = "none" # 補正は後で行う
        )
        
        # 結果にCultureを追加
        test_result$Culture <- unique(.x$Culture)
        
        return(test_result)
    })

# 結果を1つのデータフレームにまとめる
culture_results <- bind_rows(post_hoc_condition_by_culture)

# 全体でp値補正
culture_results <- culture_results %>%
    mutate(p.adj = p.adjust(p, method = "bonferroni")) # 全体でBonferroni補正

# 結果の表示
culture_results

library(tidyverse)
library(rstatix)

# 各条件（Control, Intervention1, Intervention2, Intervention3）で育成条件間のペアワイズt検定をまとめて実行
post_hoc_condition_culture_all <- data %>%
    group_by(Condition) %>% # Conditionごとにグループ化
    group_split() %>% # グループごとにデータを分割
    map(~ {
        test_result <- pairwise_t_test(
            .x, # グループごとのデータ
            formula = LDL ~ Culture, # 育成条件ごとの比較
            p.adjust.method = "none" # 補正は後で行う
        )
        
        # 結果にConditionを追加
        test_result$Condition <- unique(.x$Condition)
        
        return(test_result)
    })

# 結果を1つのデータフレームにまとめる
condition_results <- bind_rows(post_hoc_condition_culture_all)

# 全体でp値補正
condition_results <- all_results %>%
    mutate(p.adj = p.adjust(p, method = "bonferroni")) # 全体でBonferroni補正

# 結果を表示
condition_results

# p値のフォーマット関数
format_p_value <- function(p) {
    if (p < 0.001) {
        return("<0.001")
    } else {
        return(sprintf("= %.3f", p))
    }
}

# フォーマットを適用してp値を表示し、Culture毎に位置を調整
culture_results <- culture_results %>%
  mutate(
    adjusted_p_value = sapply(p.adj, format_p_value),  # p値のフォーマット
    x1 = as.numeric(factor(group1, levels = c("Control", "Intervention1", "Intervention2", "Intervention3"))),  # group1を数値化
    x2 = as.numeric(factor(group2, levels = c("Control", "Intervention1", "Intervention2", "Intervention3"))),  # group2を数値化
    xmin = ifelse(Culture == "HighGlc", x1 - 0.2, x1 + 0.2),  # HighGlcは左に、LowGlcは右に配置
    xmax = ifelse(Culture == "HighGlc", x2 - 0.2, x2 + 0.2)   # xminとxmaxで位置調整
  )
    
# フォーマットを適用してp値を表示し、位置を調整    
condition_results <- condition_results %>%
  mutate(
    adjusted_p_value = sapply(p.adj, format_p_value),
    x = as.numeric(factor(Condition, levels = unique(Condition))),  # x軸位置を数値に変換
    xmin = x - 0.2,  # xminをxから計算
    xmax = x + 0.2   # xmaxをxから計算
  )

# グラフの作成
library(ggpubr)
library(ggprism)

# ggboxplotを使用して箱ひげ図を作成
ggboxplot(
    data,
    x = "Condition", y = "LDL", color = "Culture", fill = "Culture", 
    add = "jitter"                          # ドットプロットを追加
) +
    # Cultureごとの色を指定
    scale_color_manual(values = c("black", "red")) +  # 枠線の色
    scale_fill_manual(values = c("gray30", "pink")) +  # 塗りつぶしの色

    
    # 補正したP値をグラフに書き込む(色はCulture毎で変更)
    add_pvalue(
        data = culture_results,
        xmin = "xmin",
        xmax = "xmax",
        label = "p {adjusted_p_value}",
        y.position = 170, # P値を記載するY軸の高さの設定
        textsize = 6, 
        colour = "Culture", # 育成条件でP値を色分けする設定
        tip_length = 0.01,
        step.increase = 0.1, # P値の高さが段階的に増加するように調整
        bracket.colour = "black"
    )  +
    
    # 補正したP値を治療条件毎に書き込む
    add_pvalue(
        data = condition_results,
        xmin = "xmin",
        xmax = "xmax",
        label = "p {adjusted_p_value}",
        y.position = 160, # P値を記載するY軸の高さの設定
        textsize = 6,
        bracket.color = "black"
    )+
    
    # ラベルの設定
    labs(x = "治療条件", y = "LDL値", color = "育成条件", fill = "育成条件") +
    
    # 見た目・レイアウトの調整
    theme_classic()+
    theme(
        strip.background = element_blank(),   # ファセットタイトルの背景を消す
        strip.text = element_blank(),          # ファセットタイトルのテキストを消す
        axis.line = element_line(color = "black", size = 0.3, linetype = "solid"),
        axis.text = element_text(color = "black", size = 12),
        axis.title = element_text(size = 14),
        axis.ticks.length = unit(-2, "mm"),
        axis.ticks = element_line(color = "black", size = 0.3),
    )