特ヲタ認定OT

#########################################
#ABABデザイン：誤差構造 × 分散構造最適判定コード#
#########################################

library(nlme)

library(openxlsx)

# --- 1. データ読み込み ---

dat=0

dat <- read.xlsx("ABABdata.xlsx", sheet = 1)

# --- 2. ACF / PACF（目視チェック） ---

cases <- unique(dat$case)

par(mfrow = c(2, length(cases)))

for (i in seq_along(cases)) {

  y <- dat$outcome[dat$case == cases[i]]

  acf(y, main = paste("ACF:", cases[i]))

  pacf(y, main = paste("PACF:", cases[i]))

}

# --- 3. 6モデルをフィットしてAIC比較 ---

select_best_of_6 <- function(dat) {

  fm_ind_const <- lme(outcome ~ 1 + trt, random = ~1|case,

                      data = dat, control = lmeControl(returnObject=TRUE))

  fm_ind_var   <- lme(outcome ~ 1 + trt, random = ~1|case,

                      weights = varIdent(form=~1|phase),

                      data = dat, control = lmeControl(returnObject=TRUE))

  fm_ar1_const <- lme(outcome ~ 1 + trt, random = ~1|case,

                      correlation = corAR1(0.01, ~session|case),

                      data = dat, control = lmeControl(returnObject=TRUE))

  fm_ar1_var   <- lme(outcome ~ 1 + trt, random = ~1|case,

                      correlation = corAR1(0.01, ~session|case),

                      weights = varIdent(form=~1|phase),

                      data = dat, control = lmeControl(returnObject=TRUE))

  fm_ma1_const <- lme(outcome ~ 1 + trt, random = ~1|case,

                      correlation = corARMA(0, ~session|case, p=0, q=1),

                      data = dat, control = lmeControl(returnObject=TRUE))

  fm_ma1_var   <- lme(outcome ~ 1 + trt, random = ~1|case,

                      correlation = corARMA(0, ~session|case, p=0, q=1),

                      weights = varIdent(form=~1|phase),

                      data = dat, control = lmeControl(returnObject=TRUE))

  # AIC表

  aic_tbl <- AIC(

    fm_ind_const, fm_ind_var,

    fm_ar1_const, fm_ar1_var,

    fm_ma1_const, fm_ma1_var

  )

  # 最適モデル

  best_model <- rownames(aic_tbl)[which.min(aic_tbl$AIC)]

  # 読みやすい名前

  readable <- c(

    fm_ind_const = "independent + constant variance",

    fm_ind_var   = "independent + phase variance",

    fm_ar1_const = "AR1 + constant variance",

    fm_ar1_var   = "AR1 + phase variance",

    fm_ma1_const = "MA1 + constant variance",

    fm_ma1_var   = "MA1 + phase variance"

  )

  # --- 〇×判定表 ---

  model_names <- names(readable)

  flag <- ifelse(model_names == best_model, "〇", "×")

  flag_table <- data.frame(

    Model = readable,

    Selected = flag

  )

  list(

    AIC_table = aic_tbl,

    best_model_code = best_model,

    best_model_label = readable[best_model],

    flag_table = flag_table

  )

}

# --- 4. 実行 ---

result6 <- select_best_of_6(dat)

# --- 5. 出力 ---

cat("\n--- AIC table ---\n")

print(result6$AIC_table)

cat("\n--- 〇× 判定表 ---\n")

print(result6$flag_table)

cat("\n最適モデル（コード名）：", result6$best_model_code, "\n")

cat("最適モデル（読みやすい名前）：", result6$best_model_label, "\n") ​ ​

以上です。

どうやら自己相関あたりをみているそうです。

目視的に判定もできるようですのでACFとPACFについては
グラフを確認できるようにもしました。

論文などで用いる場合は、グラフもしっかり目視で確認してください。

ちなみに判定基準はこちらです。

​

◯モデルの判定基準
# ACFがlag1から徐々に減衰 & PACFがlag1だけ強い → AR(1)

# ACFがlag1だけ突出 & PACFが徐々に減衰 → MA(1)

# ACFもPACFもほぼゼロ → independent（誤差は独立）

◯分散差の基準
# p < 0.05 → フェーズごとの分散差が有意 → Variance differs by phase を採用

# p >= 0.05 → 分散差は有意でない → Constant variance を採用