これまでt検定とMANOVAをRで実装する方法を見てきました。これらはいずれも説明変数がカテゴリカルであるという特徴がありました。しかし、実際のデータでは説明変数に連続値を取る場合があります。例えば、2種類の語彙学習方法 (Glossing vs. Writing) で学んだ語を、RecallとRecognitionの2種類のテストで測定し、学習者の語彙サイズも測定したとします。語彙学習にはもともと持っていた学習者の語彙サイズが影響すると考えられます。これは学習者のメンタルレキシコンが豊富である方が語彙のネットワークが構築されているため、新しい語彙を学習しやすいという理由です。日本人英語学習者を対象に開発された語彙サイズテストとして、Hamada et al. (2021) のVST-NJ8があります。学習指導要領では高校卒業までに5000語を学習することが明示されています。そこで、平均2500語、分散500の語彙サイズを新しくデータに付け加えます。
また、学習者の言語適正である記憶力やワーキングメモリの容量が影響する可能性もあります。語彙学習の分野ではPaul MearaのLLAMA Bが記憶力を測定するテストとして用いられています。LLAMA Bは20個のイラストの名前を2分で記憶するテストです。
語彙サイズと言語適性を共変量として組み込んだモデルを作成し、重回帰分析をRで実行してみたいと思います。
# パッケージの読み込み
library("tidyverse")
library("psych")
library("rstan")
library("brms")
library("rstantools")
library("bayesplot")
library("tidybayes")
library("loo")
library("bayestestR")
library("performance")
library("emmeans")
library("lme4")
library("lmerTest")
library("DHARMa")
library("sjPlot")
library("effectsize")
theme_set(theme_bw())
# 処理の高速化
rstan_options(auto_write = T)
options(mc.cores=parallel::detectCores())
# サンプルデータの作成
# パラメータの設定
n_participants <- 50
n_items <- 30
treatments <- c("Glossing", "Writing")
tests <- c("Recall", "Recognition")
times <- c("Immediate", "Delayed")
# データフレームの作成
set.seed(123) # 再現性のための乱数種の設定
dat <- expand.grid(Participant = 1:n_participants,
Item = 1:n_items,
Treatment = treatments,
Test = tests,
Time = times)
# アイテムの割り当てをTreatmentに基づいて修正
dat <- dat[dat$Treatment == "Glossing" & dat$Item <= 15 | dat$Treatment == "Writing" & dat$Item > 15, ]
# スコアの割り当て
dat$Score <- ifelse(dat$Treatment == "Glossing",
rbinom(nrow(dat), 1, ifelse(dat$Test == "Recognition", 0.5, 0.2)), # Glossing: Recognition > Recall
rbinom(nrow(dat), 1, ifelse(dat$Test == "Recognition", 0.9, 0.5))) # Writing: Recognition > Recall
# Timeによるスコアの調整
dat$Score <- ifelse(dat$Time == "Immediate",
dat$Score + sample(0:1, nrow(dat), replace = T), # Immediateのスコアを上昇
dat$Score)
# スコアが1を超えないように調整
dat$Score <- pmin(dat$Score, 1)
# 参加者ごとにスコアを集計
dat_summarized <- dat %>%
mutate(Participant = factor(Participant),
Item = factor(Item)) %>%
group_by(Participant, Treatment, Test, Time) %>%
summarize(Score = sum(Score), .groups = 'drop')
# データの整形
wide_scores <- dat_summarized %>%
pivot_wider(
names_from = Test,
values_from = Score,
names_prefix = "Score_"
)
# 各参加者に対してランダムなProficiencyとAptitudeの値を生成
participants <- unique(wide_scores$Participant)
proficiency_values <- rnorm(length(participants), mean = 2500, sd = sqrt(500))
aptitude_values <- rnorm(length(participants), mean = 10, sd = sqrt(10))
# ProficiencyとAptitudeデータを参加者と結合
proficiency_data <- data.frame(Participant = participants, Proficiency = proficiency_values)
aptitude_data <- data.frame(Participant = participants, Aptitude = aptitude_values)
# 元のデータフレームにProficiencyデータを結合
wide_scores <- merge(wide_scores, proficiency_data, by = "Participant")
# 結合したデータフレームにAptitudeデータを結合
wide_scores <- merge(wide_scores, aptitude_data, by = "Participant")
作成されたデータは以下の通りです。
> head(wide_scores)
Participant Treatment Time Score_Recall Score_Recognition Proficiency Aptitude
1 1 Glossing Immediate 10 11 2543.191 6.446539
2 1 Glossing Delayed 3 7 2543.191 6.446539
3 1 Writing Immediate 14 15 2543.191 6.446539
4 1 Writing Delayed 8 12 2543.191 6.446539
5 10 Glossing Immediate 3 12 2507.429 10.789955
6 10 Glossing Delayed 3 7 2507.429 10.789955
以下では共変量を含む重回帰分析を実施するコードを示します。重回帰分析における説明変数は変換を行います。連続変数はz変換、カテゴリカル変数はsum-contrastを実施します。z変換を行う目的は連続変数の単位を揃えることです。語彙サイズは平均2500、Aptitudeは平均10となるため、単位が全く異なります。そのため、重回帰分析の前にz変換を行います。カテゴリカル変数はRではダミーコーディングされることが多いですが、交互作用を検討する際には条件Aを-1、条件Bを1とし、足し合わせて0になるような設定です。
mutate関数を使うとデータセットの右側に新しく列を加える関数で、イコールの左に列の名前、右側にその内容を表します。scale関数はz変換を行うものです。
# 説明変数の変換
# 連続変数はz変換
wide_scores <- mutate(wide_scores,
Proficiency = scale(Proficiency),
Aptitude = scale(Aptitude))
カテゴリカル変数のコーディングにはいくつかの選択肢があります。sum contrastには以下のようなコーディングを行います。今回のような名義変数の場合にはダミーコーディングやサムコーディングが使えます。それぞれ2水準であることから、引数には2を用います。
# カテゴリカル変数は中心化
contrasts(wide_scores$Treatment) <- contr.sum(2)
contrasts(wide_scores$Time) <- contr.sum(2)
同じカテゴリカル変数でも順序尺度の場合はOrthogonal Polynomial Coding (contr.poly) が適切です。
contr.helmert(n, contrasts = T, sparse = F)
contr.poly(n, scores = 1:n, contrasts = T, sparse = F)
contr.sum(n, contrasts = T, sparse = F)
contr.treatment(n, base = 1, contrasts = T, sparse = F)
contr.SAS(n, contrasts = T, sparse = F)
重回帰分析をはじめとした線型モデルはlm関数を使って実装することができます。t検定と同様に、チルダ (~) の左側に目的変数、右側に説明変数を持ちます。説明変数が複数ある場合、プラスで繋ぐものは主効果を、アスタリスクで繋ぐものは交互作用をそれぞれ表しています。今回は変数選択にステップワイズ法を用いて、情報量基準でモデルを比較していきます。具体的なコードは以下の通り。