Rで完全実施要因計画

library(dplyr)
library(ggplot2)
seed <- 20250510

factor_A_levels <- factor(c("低温", "高温"), levels = c("低温", "高温"))
factor_B_levels <- factor(c("低圧", "高圧"), levels = c("低圧", "高圧"))
# factor_A_levels <- factor(c("低温", "高温"))
# factor_A_levels <- relevel(factor_A_levels,"低温")
# factor_B_levels <- factor(c("低圧", "高圧"))
# factor_B_levels <- relevel(factor_B_levels,"低圧")

# 計画行列（デザインマトリックス）の作成
doe_plan_unique <- expand.grid(
  A = factor_A_levels,
  B = factor_B_levels
)

print("基本的な実験計画 (各条件1回):")
print(doe_plan_unique)

# 各実験条件での繰り返し回数
n_replicates <- 3

# 繰り返しを考慮した実験計画を作成
# 各行をn_replicates回繰り返す
doe_plan_replicated <- doe_plan_unique[rep(seq_len(nrow(doe_plan_unique)), each = n_replicates), ]
rownames(doe_plan_replicated) <- NULL # 行名をリセット

print(paste("繰り返しを含む実験計画 (各条件", n_replicates, "回):"))
print(doe_plan_replicated)

# 応答Y（製品の収率）をシミュレートします。
# 真のモデルを以下のように仮定します:
# Y = 切片 + Aの効果 + Bの効果 + AB交互作用効果 + 誤差

# 真のパラメータを設定
intercept_true <- 60 # ベースラインの収率
effect_A_high_true <- 15 # Aが"高温"であることによる収率増
effect_B_high_true <- 10 # Bが"高圧"であることによる収率増
effect_AB_interaction_true <- 8 # Aが"高温"かつBが"高圧"であることによる追加の収率増 (交互作用)
error_sd_true <- 2.5 # 実験誤差の標準偏差

# 乱数のシードを設定
set.seed(seed)

# 応答Yを計算する関数
generate_response <- function(A_level, B_level) {
  y <- intercept_true
  if (A_level == "高温") {
    y <- y + effect_A_high_true
  }
  if (B_level == "高圧") {
    y <- y + effect_B_high_true
  }
  if (A_level == "高温" && B_level == "高圧") {
    y <- y + effect_AB_interaction_true
  }
  # 正規分布に従う誤差を加える
  y <- y + rnorm(1, mean = 0, sd = error_sd_true)
  return(y)
}

# 各実験条件に対して応答Yを生成
sim_data <- doe_plan_replicated %>%
  mutate(Y = mapply(generate_response, A, B))

print("シミュレーションデータ (応答Yを含む):")
print(sim_data)

# 因子A、因子B、およびそれらの交互作用A:Bをモデルに含めます。
# A*B は A + B + A:B と同じ。
model <- lm(Y ~ A * B, data = sim_data)

# 分析結果のサマリーを表示
print("モデルのサマリー (summary):")
summary(model)

# 推定係数の95%信頼区間を表示
print("推定係数の95%信頼区間:")
confint(object = model, level = 0.95)