# ======== ライブラリ ======== #
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 31.53 15.76 2 52 0.4264 0.6551072
## Time 634.42 634.42 1 52 17.1619 0.0001266 ***
## Condition:Time 115.18 57.59 2 52 1.5579 0.2202507
## ---
## 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: 764.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.69424 -0.48653 0.00105 0.51205 2.96362
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 61.64 7.851
## Residual 36.97 6.080
## Number of obs: 110, groups: ID, 55
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 38.6500 2.2204 74.7792 17.406 < 2e-16 ***
## ConditionModel -0.5250 3.3307 74.7792 -0.158 0.87518
## ConditionControl -0.4921 3.1812 74.7792 -0.155 0.87748
## TimePost 5.6000 1.9227 52.0000 2.913 0.00527 **
## ConditionModel:TimePost 1.2750 2.8840 52.0000 0.442 0.66026
## ConditionControl:TimePost -3.6000 2.7546 52.0000 -1.307 0.19700
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CndtnM CndtnC TimPst CnM:TP
## ConditinMdl -0.667
## CondtnCntrl -0.698 0.465
## TimePost -0.433 0.289 0.302
## CndtnMdl:TP 0.289 -0.433 -0.201 -0.667
## CndtnCnt:TP 0.302 -0.201 -0.433 -0.698 0.465
# ==== 単純効果比較関数 ==== #
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 -5.60 1.92 52 -2.913 0.0053
##
## Condition = Model:
## contrast estimate SE df t.ratio p.value
## Pre - Post -6.88 2.15 52 -3.198 0.0024
##
## Condition = Control:
## contrast estimate SE df t.ratio p.value
## Pre - Post -2.00 1.97 52 -1.014 0.3153
##
## Degrees-of-freedom method: kenward-roger
##
## 推定平均（参考）:
## Condition Time emmean SE df lower.CL upper.CL
## AI Pre 38.6 2.22 74.8 34.2 43.1
## Model Pre 38.1 2.48 74.8 33.2 43.1
## Control Pre 38.2 2.28 74.8 33.6 42.7
## AI Post 44.2 2.22 74.8 39.8 48.7
## Model Post 45.0 2.48 74.8 40.1 49.9
## Control Post 40.2 2.28 74.8 35.6 44.7
##
## 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 | 0.02 | [0.00, 1.00]
## Time | 0.25 | [0.10, 1.00]
## Condition:Time | 0.06 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].