# ======== ライブラリ ======== #
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
library(lme4)
## Warning: package 'lme4' was built under R version 4.3.3
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(emmeans)
## Warning: package 'emmeans' was built under R version 4.3.3
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(effectsize)
## Warning: package 'effectsize' was built under R version 4.3.3
# ======== データ読み込み ======== #
df <- read_excel("~/Library/CloudStorage/Dropbox/土居_関口/2.Analysis_Doi/writing_results_scored_total.xlsx")
# 条件を統一・順序付け
df <- df %>%
mutate(
Condition = factor(case_when(
条件 == "ai-wcf" ~ "AI",
条件 == "model text" ~ "Model",
条件 == "control" ~ "Control",
TRUE ~ 条件
), levels = c("AI", "Model", "Control")),
ID = as.factor(学籍番号)
)
# ======== LMEを回す関数を定義 ======== #
run_lme <- function(data, criterion_name, pre_col, post_col) {
# データを縦長化
df_long <- data %>%
select(ID, Condition, !!sym(pre_col), !!sym(post_col)) %>%
pivot_longer(
cols = c(!!sym(pre_col), !!sym(post_col)),
names_to = "Time",
values_to = "Score"
) %>%
mutate(
Time = ifelse(Time == pre_col, "Pre", "Post"),
Time = factor(Time, levels = c("Pre", "Post"))
)
# 混合効果モデル
model <- lmer(Score ~ Condition * Time + (1 | ID), data = df_long)
cat("\n\n===== LME結果:", criterion_name, "=====\n")
print(anova(model, type = 3))
print(summary(model))
# EMM
emm <- emmeans(model, ~ Condition * Time)
# プロット
emm_df <- as.data.frame(emm)
p <- ggplot(emm_df, aes(x = Time, y = emmean, group = Condition, color = Condition)) +
geom_line(linewidth = 1.2) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.1, linewidth = 0.8) +
scale_color_manual(values = c("AI" = "#1f78b4", "Model" = "#33a02c", "Control" = "#e31a1c")) +
labs(
title = paste("Interaction Plot:", criterion_name, "(Condition × Time)"),
x = "Time (Pre / Post)",
y = "Estimated Mean Score (0–100)",
color = "Condition"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
axis.title = element_text(face = "bold"),
axis.text = element_text(size = 12),
legend.position = "top",
legend.title = element_text(face = "bold")
) +
coord_cartesian(ylim = c(0, 100))
print(p)
}
# ======== 総合スコアでLME実行 ======== #
run_lme(df, "Total Writing Score (0–100)", "Pre_Total", "Post_Total")
##
##
## ===== LME結果: Total Writing Score (0–100) =====
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Condition 2.93 1.47 2 91 0.0455 0.9556
## Time 1522.64 1522.64 1 91 47.2366 7.653e-10 ***
## Condition:Time 121.88 60.94 2 91 1.8905 0.1569
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Score ~ Condition * Time + (1 | ID)
## Data: df_long
##
## REML criterion at convergence: 1349.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09055 -0.46505 -0.03632 0.47142 2.39686
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 100.42 10.021
## Residual 32.23 5.678
## Number of obs: 188, groups: ID, 94
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 34.0000 2.1028 115.6991 16.169 < 2e-16 ***
## ConditionModel -0.5000 2.9269 115.6991 -0.171 0.865
## ConditionControl 0.6875 2.9269 115.6991 0.235 0.815
## TimePost 6.3333 1.4659 91.0000 4.320 3.97e-05 ***
## ConditionModel:TimePost 0.9167 2.0405 91.0000 0.449 0.654
## ConditionControl:TimePost -2.8333 2.0405 91.0000 -1.389 0.168
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CndtnM CndtnC TimPst CnM:TP
## ConditinMdl -0.718
## CondtnCntrl -0.718 0.516
## TimePost -0.349 0.250 0.250
## CndtnMdl:TP 0.250 -0.349 -0.180 -0.718
## CndtnCnt:TP 0.250 -0.180 -0.349 -0.718 0.516
# ==== 単純効果比較関数 ==== #
compare_prepost_by_condition <- function(data, criterion_name, pre_col, post_col) {
df_long <- data %>%
select(ID, Condition, !!sym(pre_col), !!sym(post_col)) %>%
pivot_longer(
cols = c(!!sym(pre_col), !!sym(post_col)),
names_to = "Time",
values_to = "Score"
) %>%
mutate(
Time = ifelse(Time == pre_col, "Pre", "Post"),
Time = factor(Time, levels = c("Pre", "Post"))
)
# LMEモデル
model <- lmer(Score ~ Condition * Time + (1 | ID), data = df_long)
emm <- emmeans(model, ~ Condition * Time)
cat("\n\n===== 単純効果比較: ", criterion_name, "=====\n")
# 各群内のPre vs Post
result <- pairs(emm, by = "Condition", adjust = "tukey")
print(result)
cat("\n推定平均（参考）:\n")
print(summary(emm))
}
# 総合スコアで実行
compare_prepost_by_condition(df, "Total Writing Score (0–100)", "Pre_Total", "Post_Total")
##
##
## ===== 単純効果比較: Total Writing Score (0–100) =====
## Condition = AI:
## contrast estimate SE df t.ratio p.value
## Pre - Post -6.33 1.47 91 -4.320 <.0001
##
## Condition = Model:
## contrast estimate SE df t.ratio p.value
## Pre - Post -7.25 1.42 91 -5.108 <.0001
##
## Condition = Control:
## contrast estimate SE df t.ratio p.value
## Pre - Post -3.50 1.42 91 -2.466 0.0155
##
## Degrees-of-freedom method: kenward-roger
##
## 推定平均（参考）:
## Condition Time emmean SE df lower.CL upper.CL
## AI Pre 34.0 2.10 116 29.8 38.2
## Model Pre 33.5 2.04 116 29.5 37.5
## Control Pre 34.7 2.04 116 30.7 38.7
## AI Post 40.3 2.10 116 36.2 44.5
## Model Post 40.8 2.04 116 36.7 44.8
## Control Post 38.2 2.04 116 34.2 42.2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
# --- データを縦長に変換（Total Score用） ---
df_long <- df %>%
select(ID, Condition, Pre_Total, Post_Total) %>%
pivot_longer(
cols = c(Pre_Total, Post_Total),
names_to = "Time",
values_to = "Score"
) %>%
mutate(
Time = factor(ifelse(Time == "Pre_Total", "Pre", "Post"), levels = c("Pre", "Post")),
Condition = factor(Condition, levels = c("AI", "Model", "Control"))
)
# --- LMEモデル構築 ---
model_total <- lmer(Score ~ Condition * Time + (1 | ID), data = df_long)
# --- 効果量算出 ---
eta_sq <- eta_squared(model_total, partial = TRUE)
cat("\n===== 効果量 (partial η²): Total Score (0–100) =====\n")
##
## ===== 効果量 (partial η²): Total Score (0–100) =====
print(eta_sq)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Eta2 (partial) | 95% CI
## ----------------------------------------------
## Condition | 9.98e-04 | [0.00, 1.00]
## Time | 0.34 | [0.22, 1.00]
## Condition:Time | 0.04 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].