まず練習にあたり必要なパッケージをインストールしましょう。
ここはよく分からなくても大丈夫です。
```R
# 必要なパッケージをインストール・読み込み
if (!require("lme4")) install.packages("lme4")
if (!require("ggplot2")) install.packages("ggplot2")
library(lme4)
library(lmerTest)
library(ggplot2)
```
20人分のデータを作成します。
```R
# 20人の被験者データ、各被験者につき2回の測定値を直接入力
# ID: 被験者ID
# visit: 測定回 (1回目か2回目)
# systolic: 収縮期血圧
bp_data <- data.frame(
ID = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5,
6, 6, 7, 7, 8, 8, 9, 9, 10, 10,
11, 11, 12, 12, 13, 13, 14, 14, 15, 15,
16, 16, 17, 17, 18, 18, 19, 19, 20, 20),
visit = c(1, 2, 1, 2, 1, 2, 1, 2, 1, 2,
1, 2, 1, 2, 1, 2, 1, 2, 1, 2,
1, 2, 1, 2, 1, 2, 1, 2, 1, 2,
1, 2, 1, 2, 1, 2, 1, 2, 1, 2),
systolic = c(128, 135, 115, 118, 142, 138, 125, 129, 137, 142,
118, 122, 130, 133, 140, 145, 119, 122, 131, 134,
147, 152, 122, 125, 139, 136, 133, 138, 126, 129,
145, 150, 129, 133, 136, 133, 122, 128, 138, 143)
)
```
データを確認しましょう。
```R
print(bp_data)
```
このようなデータができます。
![[Pasted image 20250330054915.png|200]]
データ型を変更します
```R
bp_data$ID <- factor(bp_data$ID) # IDを因子型に変換
bp_data$visit <- factor(bp_data$visit) # visitを因子型に変換
```
プロットしましょう。各個人について、2回目の血圧のほうが高くなっていることが分かります。
```R
# 各被験者の2回の測定値を線で繋いだプロット
ggplot(bp_data, aes(x = factor(visit), y = systolic, group = ID, color = factor(ID))) +
geom_point() +
geom_line() +
labs(title = "Change in Systolic Blood Pressure per Subject",
x = "Visit",
y = "Systolic Blood Pressure (mmHg)",
color = "Subject ID") +
theme_minimal()
```
![[Pasted image 20250330055952.png|400]]
記述統計もしておきましょう
```R
# 記述統計
cat("全体の平均収縮期血圧:", mean(bp_data$systolic), "mmHg\n")
cat("1回目の平均収縮期血圧:", mean(bp_data$systolic[bp_data$visit == 1]), "mmHg\n")
cat("2回目の平均収縮期血圧:", mean(bp_data$systolic[bp_data$visit == 2]), "mmHg\n")
```
いよいよ混合効果モデルを作ってみましょう
```R
# IDをランダム効果として、visitを固定効果として扱う
# ランダム切片モデル
random_intercept_model <- lmer(systolic ~ visit + (1|ID), data = bp_data)
summary(random_intercept_model)
```
実行結果です。
```
> summary(random_intercept_model)
Linear mixed model fit by REML. t-tests use Satterthwaite's method
[lmerModLmerTest]
Formula: systolic ~ visit + (1 | ID)
Data: bp_data
REML criterion at convergence: 239.7
Scaled residuals:
Min 1Q Median 3Q Max
-1.5835 -0.2813 -0.0923 0.4833 1.7719
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 80.553 8.975
Residual 4.541 2.131
Number of obs: 40, groups: ID, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 131.1000 2.0627 20.0409 63.558 < 2e-16 ***
visit2 3.1500 0.6739 19.0000 4.675 0.000165 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
visit2 -0.163
```
結果が色々でてきています。論文に使われる結果として必要なのは以下の部分です。
```
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 131.1000 2.0627 20.0409 63.558 < 2e-16 ***
visit2 3.1500 0.6739 19.0000 4.675 0.000165 ***
```
`visit1` が `(Intercept)` です。血圧が `131.1` であることが分かります。
`visit2` で、血圧が `3.15` 高くなっていることが分かります。
固定効果の信頼区間やランダム効果を確認しましょう。結果は各自確認してください。
```R
# 固定効果の信頼区間
confint(random_intercept_model, method = "Wald")
# ランダム効果の確認
ranef(random_intercept_model)
```
シンプルなｔ検定をしてみましょう。こちらも結果は各自確認してください。
```R
# 単純な対応のあるt検定
t_test_result <- t.test(systolic ~ visit, data = bp_data, paired = TRUE)
print(t_test_result)
```