---
title: "AI支援による評価報告書平易化RCTの再現分析"
author: "Masatsugu Mikami"
date: "2026-07-08"
output:
  html_document:
    toc: true
    toc_depth: 3
    number_sections: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 5
)

seed_value <- 123
target_r_version <- "4.5.1"
runtime_r_version <- paste(R.version$major, R.version$minor, sep = ".")

required_packages <- c(
  "readr", "dplyr", "forcats", "tibble", "tidyr", "purrr",
  "car", "emmeans", "effectsize", "broom", "mediation", "glmnet",
  "MASS", "knitr"
)

missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]

if (length(missing_packages) > 0) {
  stop(
    "Install the following packages before knitting this document: ",
    paste(missing_packages, collapse = ", ")
  )
}

invisible(lapply(required_packages, library, character.only = TRUE))
set.seed(seed_value)

if (!identical(runtime_r_version, target_r_version)) {
  warning(
    "This public analysis is specified for R version ",
    target_r_version,
    "; the current runtime is R version ",
    runtime_r_version,
    "."
  )
}
```

# 概要

この文書は、論文で報告したランダム化比較試験（RCT）の分析を再現するための
R Markdownです。査読対応として追加した以下の頑健性確認も含めています。

- stepwise-AICモデルに対するLASSO正則化による感度分析。
- `mediation::medsens`を用いた媒介感度分析。
- 5件法アウトカムを順序尺度として扱う順序ロジット感度分析。
- `Simplification x Document Type`交互作用による文書タイプ別の異質性確認。

このR Markdownファイルと同じディレクトリに`rawdata.csv`が置かれていることを
前提としています。

すべての分析はR v.4.5.1での実行を想定しています。乱数を用いる手続きでは、
下記の固定シードを用いています。

# 計算環境

```{r computational-environment}
knitr::kable(
  tibble::tibble(
    item = c("対象Rバージョン", "実行時Rバージョン", "乱数シード"),
    value = c(
      paste0("R v.", target_r_version),
      paste0("R v.", runtime_r_version),
      seed_value
    )
  )
)
```

# データ準備

```{r data-import}
rawdata <- readr::read_csv("rawdata.csv", show_col_types = FALSE)

education_levels <- c(
  "Junior High School",
  "High School (In progress)",
  "High School Graduate",
  "Vocational School (In progress)",
  "Vocational School Graduate",
  "Junior College (In progress)",
  "Junior College Graduate",
  "University (In progress)",
  "University Graduate",
  "Graduate School (Master's - In progress)",
  "Graduate School (Master's - Completed)",
  "Graduate School (Ph.D. - In progress)",
  "Graduate School (Ph.D. - Completed)",
  "Other",
  "Prefer not to say"
)

occupation_levels <- c(
  "Student",
  "Homemaker",
  "Self-employed",
  "Public Official",
  "Teaching Staff",
  "Medical/Healthcare",
  "Business Owner",
  "Corporate Executive",
  "Employee",
  "Part-time/Temporary",
  "Unemployed",
  "Agriculture/Forestry/Fisheries",
  "Professional (Lawyer/Tax Accountant, etc.)",
  "Other"
)

prefecture_levels <- c(
  "Hokkaido", "Aomori", "Iwate", "Miyagi", "Akita", "Yamagata",
  "Fukushima", "Ibaraki", "Tochigi", "Gunma", "Saitama", "Chiba",
  "Tokyo", "Kanagawa", "Niigata", "Toyama", "Ishikawa", "Fukui",
  "Yamanashi", "Nagano", "Gifu", "Shizuoka", "Aichi", "Mie",
  "Shiga", "Kyoto", "Osaka", "Hyogo", "Nara", "Wakayama",
  "Tottori", "Shimane", "Okayama", "Hiroshima", "Yamaguchi",
  "Tokushima", "Kagawa", "Ehime", "Kochi", "Fukuoka", "Saga",
  "Nagasaki", "Kumamoto", "Oita", "Miyazaki", "Kagoshima",
  "Okinawa"
)

analysis_data <- rawdata |>
  dplyr::transmute(
    respondent_id = respid,
    document_id = Quota1,
    group = dplyr::case_when(
      Quota1 == 1 ~ "A1 (統制群)",
      Quota1 == 2 ~ "A2 (介入群)",
      Quota1 == 3 ~ "B1 (統制群)",
      Quota1 == 4 ~ "B2 (介入群)",
      TRUE ~ NA_character_
    ),
    document_type = factor(
      dplyr::if_else(Quota1 %in% c(3, 4), "B", "A"),
      levels = c("A", "B")
    ),
    simplification = ifelse(Quota1 %in% c(2, 4), 1, 0),
    simplification_factor = factor(
      ifelse(Quota1 %in% c(2, 4), "Simplified", "Original"),
      levels = c("Original", "Simplified")
    ),
    comprehension = Q2,
    acceptance = Q3,
    acceptance_ord = ordered(Q3, levels = 1:5),
    age = SC2_1,
    gender = factor(SC1, levels = c(1, 2), labels = c("Male", "Female")),
    education = factor(Q1, levels = seq_along(education_levels), labels = education_levels),
    occupation = factor(SC5, levels = seq_along(occupation_levels), labels = occupation_levels),
    prefecture = factor(SC3, levels = seq_along(prefecture_levels), labels = prefecture_levels)
  ) |>
  dplyr::mutate(
    group = factor(
      group,
      levels = c("A1 (統制群)", "B1 (統制群)", "A2 (介入群)", "B2 (介入群)")
    )
  ) |>
  tidyr::drop_na() |>
  droplevels()

knitr::kable(
  tibble::tibble(
    n = nrow(analysis_data),
    mean_age = mean(analysis_data$age),
    sd_age = sd(analysis_data$age),
    male_n = sum(analysis_data$gender == "Male"),
    female_n = sum(analysis_data$gender == "Female")
  ),
  digits = 3
)
```

# 記述統計とランダム化チェック

```{r group-counts}
group_counts <- analysis_data |>
  dplyr::count(group, name = "n")

knitr::kable(group_counts)
```

```{r table-primary-variables}
format_count_pct <- function(x, denominator) {
  paste0(x, " (", round(100 * x / denominator, 1), "%)")
}

likert_table <- function(data, variable, label) {
  counts <- table(data$group, data[[variable]])
  out <- as.data.frame.matrix(t(counts))
  out$level <- rownames(out)
  out$variable <- label
  out <- tibble::as_tibble(out)
  out |>
    dplyr::relocate(variable, level) |>
    dplyr::mutate(
      dplyr::across(
        dplyr::all_of(levels(data$group)),
        ~ format_count_pct(.x, sum(.x))
      )
    )
}

acceptance_table <- likert_table(analysis_data, "acceptance", "納得度")
comprehension_table <- likert_table(analysis_data, "comprehension", "理解度")

knitr::kable(
  dplyr::bind_rows(acceptance_table, comprehension_table),
  caption = "実験条件別の度数と割合。"
)
```

```{r randomization-checks}
set.seed(seed_value)

randomization_checks <- tibble::tibble(
  variable = c("年齢", "性別", "学歴", "職業", "居住都道府県"),
  test = c(
    "Kruskal-Wallis検定",
    "Fisher正確検定（シミュレーションp値）",
    "Fisher正確検定（シミュレーションp値）",
    "Fisher正確検定（シミュレーションp値）",
    "Fisher正確検定（シミュレーションp値）"
  ),
  p_value = c(
    kruskal.test(age ~ group, data = analysis_data)$p.value,
    fisher.test(table(analysis_data$gender, analysis_data$group), simulate.p.value = TRUE, B = 100000)$p.value,
    fisher.test(table(analysis_data$education, analysis_data$group), simulate.p.value = TRUE, B = 100000)$p.value,
    fisher.test(table(analysis_data$occupation, analysis_data$group), simulate.p.value = TRUE, B = 100000)$p.value,
    fisher.test(table(analysis_data$prefecture, analysis_data$group), simulate.p.value = TRUE, B = 100000)$p.value
  )
)

knitr::kable(randomization_checks, digits = 3)
```

シミュレーションによるFisher正確検定のp値は乱数に依存します。正確に再現するには、
上記の乱数種とシミュレーション回数を用います。

```{r distributional-tests}
distributional_tests <- tibble::tibble(
  outcome = c("納得度", "理解度"),
  test = "Kruskal-Wallis検定",
  p_value = c(
    kruskal.test(acceptance ~ group, data = analysis_data)$p.value,
    kruskal.test(comprehension ~ group, data = analysis_data)$p.value
  )
)

knitr::kable(distributional_tests, digits = 6)
```

したがって、理解度の分布差は`p < .001`として報告します。

# ANCOVA

論文で報告したANCOVAは、元のR Markdown分析と同じく、`aov()`による逐次平方和
（Type I）に基づいています。

```{r ancova}
ancova_model <- aov(
  acceptance ~ document_type + simplification_factor + age + gender +
    occupation + education + prefecture,
  data = analysis_data
)

summary(ancova_model)
effectsize::eta_squared(ancova_model, partial = TRUE)
```

```{r ancova-assumptions}
shapiro_acceptance <- shapiro.test(analysis_data$acceptance)
levene_acceptance <- car::leveneTest(acceptance ~ group, data = analysis_data)

shapiro_acceptance
levene_acceptance
```

```{r ancova-emmeans}
emmeans_simplification <- emmeans::emmeans(
  ancova_model,
  specs = ~ simplification_factor,
  rg.limit = 75000
)

emmeans_simplification
```

# 回帰モデル

```{r full-regression-without-comprehension}
regression_without_comprehension <- lm(
  acceptance ~ document_type + simplification_factor + age + gender +
    occupation + education + prefecture,
  data = analysis_data
)

summary(regression_without_comprehension)
confint(regression_without_comprehension)["simplification_factorSimplified", ]
car::vif(regression_without_comprehension)
```

調整済みGVIFの最大値は以下のとおりです。

```{r gvif-maximum}
vif_without_comprehension <- car::vif(regression_without_comprehension)
max_adjusted_gvif <- max(vif_without_comprehension[, "GVIF^(1/(2*Df))"])
max_adjusted_gvif
```

この結果は、多重共線性が低いことを示しています。調整済みGVIFの最大値は約1.33です。

```{r full-regression-with-comprehension}
regression_with_comprehension <- lm(
  acceptance ~ document_type + comprehension + simplification_factor + age +
    gender + occupation + education + prefecture,
  data = analysis_data
)

summary(regression_with_comprehension)
confint(regression_with_comprehension)[
  c("comprehension", "simplification_factorSimplified", "age"),
]
effectsize::standardize_parameters(regression_with_comprehension)
```

```{r stepwise-aic}
stepwise_model <- step(
  regression_with_comprehension,
  direction = "both",
  trace = 0
)

formula(stepwise_model)
AIC(stepwise_model)
```

```{r parsimonious-model-table}
parsimonious_model <- lm(
  acceptance ~ comprehension + simplification_factor + age,
  data = analysis_data
)

parsimonious_summary <- summary(parsimonious_model)
parsimonious_confint <- confint(parsimonious_model)

standardized_beta <- {
  x <- model.matrix(parsimonious_model)[, -1]
  b <- coef(parsimonious_model)[-1]
  sy <- sd(parsimonious_model$model[[1]])
  sx <- apply(x, 2, sd)
  b * sx / sy
}

table2 <- broom::tidy(parsimonious_model, conf.int = TRUE) |>
  dplyr::filter(term != "(Intercept)") |>
  dplyr::mutate(
    standardized_beta = standardized_beta[term],
    p_value = p.value
  ) |>
  dplyr::select(term, estimate, standardized_beta, conf.low, conf.high, p_value)

knitr::kable(table2, digits = 3)

tibble::tibble(
  statistic = c("Adjusted R-squared", "F", "df1", "df2", "p", "N"),
  value = c(
    parsimonious_summary$adj.r.squared,
    unname(parsimonious_summary$fstatistic[1]),
    unname(parsimonious_summary$fstatistic[2]),
    unname(parsimonious_summary$fstatistic[3]),
    pf(
      parsimonious_summary$fstatistic[1],
      parsimonious_summary$fstatistic[2],
      parsimonious_summary$fstatistic[3],
      lower.tail = FALSE
    ),
    nobs(parsimonious_model)
  )
) |>
  knitr::kable(digits = 3)
```

# 媒介分析

媒介モデルでは、年齢、性別、学歴、職業、居住都道府県を共変量として投入します。
`mediate()`のbootstrapでまれな因子水準がリサンプルから落ちることによるエラーを
避けるため、因子共変量は明示的にダミー変数化しています。これは、同じ因子を
`lm()`に直接投入することと代数的に等価です。

```{r mediation-models, cache=TRUE}
covariate_matrix <- model.matrix(
  ~ age + gender + education + occupation + prefecture,
  data = analysis_data
)[, -1]

covariate_matrix <- covariate_matrix[
  ,
  apply(covariate_matrix, 2, sd) > 0,
  drop = FALSE
]
colnames(covariate_matrix) <- make.names(colnames(covariate_matrix), unique = TRUE)

mediation_data <- data.frame(
  acceptance = analysis_data$acceptance,
  comprehension = analysis_data$comprehension,
  simplification = analysis_data$simplification,
  covariate_matrix,
  check.names = FALSE
)

rhs_covariates <- paste(colnames(covariate_matrix), collapse = " + ")

mediator_model <- lm(
  as.formula(paste("comprehension ~ simplification +", rhs_covariates)),
  data = mediation_data
)

outcome_model <- lm(
  as.formula(paste("acceptance ~ simplification + comprehension +", rhs_covariates)),
  data = mediation_data
)

coef(summary(mediator_model))["simplification", , drop = FALSE]
coef(summary(outcome_model))[c("comprehension", "simplification"), , drop = FALSE]

set.seed(seed_value)
mediation_fit <- mediation::mediate(
  model.m = mediator_model,
  model.y = outcome_model,
  treat = "simplification",
  mediator = "comprehension",
  control.value = 0,
  treat.value = 1,
  boot = TRUE,
  sims = 2000
)

summary(mediation_fit)

mediation_report <- tibble::tibble(
  effect = c("ACME", "ADE", "総合効果", "媒介割合"),
  estimate = c(
    mediation_fit$d0,
    mediation_fit$z0,
    mediation_fit$tau.coef,
    mediation_fit$n0
  ),
  ci_low = c(
    mediation_fit$d0.ci[1],
    mediation_fit$z0.ci[1],
    mediation_fit$tau.ci[1],
    mediation_fit$n0.ci[1]
  ),
  ci_high = c(
    mediation_fit$d0.ci[2],
    mediation_fit$z0.ci[2],
    mediation_fit$tau.ci[2],
    mediation_fit$n0.ci[2]
  ),
  p_value = c(
    mediation_fit$d0.p,
    mediation_fit$z0.p,
    mediation_fit$tau.p,
    mediation_fit$n0.p
  )
)

knitr::kable(mediation_report, digits = 4)
```

上記の乱数種に基づく再実行でも、ADEは統計的に有意ではありません。

# 媒介感度分析

```{r mediation-sensitivity, cache=TRUE}
sensitivity_fit <- mediation::medsens(
  mediation_fit,
  rho.by = 0.01,
  effect.type = "indirect"
)

summary(sensitivity_fit)

sensitivity_values <- tibble::tibble(
  rho_at_acme_zero = sensitivity_fit$err.cr.d,
  R2star_product = sensitivity_fit$R2star.d.thresh,
  R2tilde_product = sensitivity_fit$R2tilde.d.thresh,
  r2_mediator_model = sensitivity_fit$r.square.m,
  r2_outcome_model = sensitivity_fit$r.square.y
)

knitr::kable(sensitivity_values, digits = 3)
```

補遺用の任意プロット:

```{r mediation-sensitivity-plot, eval=FALSE}
plot(
  sensitivity_fit,
  sens.par = "rho",
  main = "媒介感度分析"
)
```

# LASSO感度分析

```{r lasso-sensitivity}
x_lasso <- model.matrix(
  acceptance ~ simplification + comprehension + age + gender +
    education + occupation + prefecture,
  data = analysis_data
)[, -1]

y_lasso <- analysis_data$acceptance

set.seed(seed_value)
cv_lasso <- glmnet::cv.glmnet(
  x = x_lasso,
  y = y_lasso,
  alpha = 1,
  nfolds = 10,
  standardize = TRUE
)

nonzero_coefs <- function(cvfit, s) {
  b <- as.matrix(coef(cvfit, s = s))
  out <- tibble::tibble(
    term = rownames(b),
    coefficient = as.numeric(b[, 1])
  )
  dplyr::filter(out, term != "(Intercept)", coefficient != 0)
}

lasso_lambda_min <- nonzero_coefs(cv_lasso, "lambda.min")
lasso_lambda_1se <- nonzero_coefs(cv_lasso, "lambda.1se")

cv_lasso$lambda.min
knitr::kable(lasso_lambda_min, digits = 4)

cv_lasso$lambda.1se
knitr::kable(lasso_lambda_1se, digits = 4)

stepwise_terms <- c("comprehension", "simplification", "age")

tibble::tibble(
  term = stepwise_terms,
  retained_at_lambda_min = stepwise_terms %in% lasso_lambda_min$term,
  retained_at_lambda_1se = stepwise_terms %in% lasso_lambda_1se$term
) |>
  knitr::kable()
```

`lambda.min`では、stepwise-AICによる節約モデルに含まれる3つの予測変数が
すべて保持されます。一方、より保守的な`lambda.1se`では理解度のみが保持されます。
これは、平易化と年齢の係数の大きさが相対的に小さいことと整合的です。

# 順序ロジット感度分析

```{r ordered-logit}
ordered_logit_model <- MASS::polr(
  acceptance_ord ~ comprehension + simplification + age,
  data = analysis_data,
  method = "logistic",
  Hess = TRUE
)

ordered_logit_table <- coef(summary(ordered_logit_model))
ordered_logit_p <- 2 * pnorm(
  abs(ordered_logit_table[, "t value"]),
  lower.tail = FALSE
)

ordered_logit_results <- tibble::tibble(
  term = rownames(ordered_logit_table),
  estimate = ordered_logit_table[, "Value"],
  std_error = ordered_logit_table[, "Std. Error"],
  t_value = ordered_logit_table[, "t value"],
  p_value = ordered_logit_p,
  odds_ratio = c(exp(coef(ordered_logit_model)), rep(NA_real_, 4))
)

knitr::kable(ordered_logit_results, digits = 4)
```

# 文書タイプ別の異質性チェック

```{r document-type-interaction}
interaction_model <- lm(
  acceptance ~ simplification * document_type + age + gender +
    occupation + education + prefecture,
  data = analysis_data
)

main_effect_model <- lm(
  acceptance ~ simplification + document_type + age + gender +
    occupation + education + prefecture,
  data = analysis_data
)

interaction_model_with_comprehension <- lm(
  acceptance ~ comprehension + simplification * document_type + age +
    gender + occupation + education + prefecture,
  data = analysis_data
)

main_effect_model_with_comprehension <- lm(
  acceptance ~ comprehension + simplification + document_type + age +
    gender + occupation + education + prefecture,
  data = analysis_data
)

interaction_terms <- grep(
  "simplification.*:document_type|document_type.*:simplification",
  rownames(coef(summary(interaction_model)))
)

interaction_terms_with_comprehension <- grep(
  "simplification.*:document_type|document_type.*:simplification",
  rownames(coef(summary(interaction_model_with_comprehension)))
)

coef(summary(interaction_model))[interaction_terms, , drop = FALSE]
anova(main_effect_model, interaction_model)

coef(summary(interaction_model_with_comprehension))[
  interaction_terms_with_comprehension,
  ,
  drop = FALSE
]
anova(main_effect_model_with_comprehension, interaction_model_with_comprehension)
```

# 追加確認の要約

追加確認の結果は、主要な解釈を支持しています。

- 納得度と理解度の分布差はいずれも`p < .001`で統計的に有意です。
- 多重共線性は低く、調整済みGVIFの最大値は約1.33です。
- 媒介分析の再実行でも、ADEは統計的に有意ではありません。
- LASSOでは、`lambda.min`で節約モデルの3変数がすべて保持されます。一方、
  より保守的な`lambda.1se`では理解度のみが保持されます。
- 順序ロジットと文書タイプ別交互作用の確認は、実験結果の主要な解釈を変更しません。

# セッション情報

以下は、このR Markdownをknitした時点の実行環境です。公開用分析の対象環境は
R v.4.5.1です。

```{r session-info}
sessionInfo()
```