はじめに

一般化線型モデル (Generalized Linear Model; GLM) は線型モデル (Linear Model; LM) を正規分布以外の分布に拡張した統計モデルです。具体的には、GLMで扱う確率分布には二項分布、ポアソン分布、負の二項分布、ベータ分布、ガンマ分布などがあります。さらに、GLMに変量効果を含むモデルとして、一般化線形混合効果モデル (Generalized Linear Mixed Effect Model; GLMM) があります。

近年の英語教育学研究や心理言語学研究では、混合効果モデルを用いた分析が多く報告されています ( e.g., Baayen et al., 2008; Brysbaert & Stevens, 2018; Cunnings, 2012; Gries, 2021; Linck & Cunnings, 2015; Meteyard & Davies, 2020)。具体的には、二項分布に基づくロジスティック回帰分析、ポアソン分布もしくは負の二項分布に基づくポアソン回帰分析を扱うことが多いです。近年では目的変数が多値となるデータに対応するGLMMを扱った研究も登場しています。ここでは目的変数が2値となるロジスティック回帰分析について見ていきましょう。

2値ロジスティック回帰分析が用いられている理由

GLMMは分析手法の総称であり、目的変数の特徴に応じて用いられる具体的な名前が異なります。例えば、目的変数が語彙性判断課題における反応時間のような連続値の場合、反応時間として得られたデータを対数変換して重回帰分析 (ランダム効果を含むモデルでは線形混合効果モデルと呼ばれます) を実施します。詳細は「5. 一般化線形モデル・一般化加法モデル」のページも参照。

線形混合効果モデルでは右のようなデータセットが得られます。このようなデータをANOVAで分析する場合、参加者ごとに平均の反応時間を用いる方法 (F1分析) と項目ごとに平均反応時間を計算する方法 (F2分析) があります。混合効果モデルではこれらのF1分析とF2分析を1つの分析で両方実施可能であるというメリットがあります。詳細は「13. 混合効果モデルとサンプルサイズの設計」のページも参照。

第二言語習得研究で利用される混合効果モデルで有名なものの1つに2値ロジスティック回帰分析があります。2値という名前が示す通り、正解・不正解のような言語テストのデータの分析に主に用いられます。特に、語彙習得研究では、ItemとParticipantの2つのランダム効果を含めた混合効果モデルでは多くの場合、目的変数が2値となります。そのため、テストの合計点を分析する場合には線形モデルであるANOVAが用いられましたが、Itemの変量効果を含めて分析する場合には一般化線形モデルである2値ロジスティック回帰分析を用いることとなります。

データの作成

以下では、Rを用いてランダム効果を取り入れたロジスティック回帰分析を実施する方法を紹介します。データセットは、50人の参加者が2種類の語彙学習方法を用いて比較された架空のデータです。この研究では、目標語彙に集中する学習方法がより効果的であると仮定します。具体的には、注釈付きの読解よりも、目標語彙を使用して文章を作る方が学習効果が高いと考えられます。参加者は15個ずつの目標語を2つの方法で学び、L2からL1への翻訳テスト (recall test) と選択問題 (recognition test) でその知識を評価します。このテストは学習直後と1週間後の遅延で実施しました。この実験設定を前提に、以下の分析手法を提案します：

1. t検定：各学習方法で15語ずつ学習し、15点満点のテストで評価します。再生テストと再認テストの両方でt検定を行います。
2. 多変量分散分析 (Multivariate Analysis of Variance; MANOVA)：再生テストと再認テストの両方を目的変数とした多変量分散分析を実施します。
3. ランダム効果を含むロジスティック回帰分析：目標語と参加者をランダム効果とし、各目標語を観測値として、テストの正解・不正解を2値データで分析します。再生テストと再認テストそれぞれに対してロジスティック回帰分析を行い、結果を報告します。
4. 多変量ロジスティック回帰分析：多変量分散分析のアプローチを用いて、ランダム効果を含むロジスティック回帰分析を実施します。これはRのbrmsパッケージを使用して実装できます。

以下のコードを実行してデータセットを作成します。実際に作成されたデータセットはこちらです。

sample.csv

# パッケージの読み込み
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)

# csvとしてデータを出力
write.csv(dat, "sample.csv")

作成されたデータは50名の参加者に対して、2つの介入 (GlossingとWriting) を行っています。Glossingでは15語、Writingでは15語を学び、合計30語の目標語を学習します。その成果を直後と遅延で測定します。測定方法にはRecallとRecognitionの2種類のテストがあります。

このデータは以下のようにScoreの列が1か0の2値データとなっています。

> head(dat)
Participant Item Treatment Test Time Score
1 1 1 Glossing Recall Immediate 1
2 2 1 Glossing Recall Immediate 1
3 3 1 Glossing Recall Immediate 1
4 4 1 Glossing Recall Immediate 1
5 5 1 Glossing Recall Immediate 1
6 6 1 Glossing Recall Immediate 1

共変量の追加

重回帰分析のRコードと同様に、共変量としてProficiencyとAptitudeを加えます。

# Participantのユニークな値を抽出
participants <- unique(dat$Participant)

# 各Participantに対してProficiencyとAptitudeを生成
proficiency_values <- rnorm(length(participants), mean = 2500, sd = sqrt(500))
aptitude_values <- rnorm(length(participants), mean = 10, sd = sqrt(10))

# 新しいデータフレームの作成
participant_data <- data.frame(Participant = participants, Proficiency = proficiency_values, Aptitude = aptitude_values)

# 元のデータフレームに結合
dat <- merge(dat, participant_data, by = "Participant")

# データの変換
dat <- mutate(dat,
Participant = factor(Participant),
Item = factor(Item),
Proficiency = scale(Proficiency),
Aptitude = scale(Aptitude))

# RecallとRecognitionテストのデータの作成
dat_recall <- dat %>%
dplyr::filter(Test == "Recall")
dat_recognition <- dat %>%
dplyr::filter(Test == "Recall")

作成されたデータは以下のとおりです。今回はMANOVAと異なり、recallとrecognitionそれぞれにモデルを作成します。そのためデータは別々のオブジェクトとして格納しました。

> head(dat_recall)
Participant Item Treatment Test Time Score Proficiency Aptitude
1 1 1 Glossing Recall Immediate 1 1.822503 -0.2371038
2 1 18 Writing Recall Immediate 1 1.822503 -0.2371038
3 1 26 Writing Recall Delayed 1 1.822503 -0.2371038
4 1 14 Glossing Recall Immediate 0 1.822503 -0.2371038
5 1 5 Glossing Recall Delayed 0 1.822503 -0.2371038
6 1 10 Glossing Recall Immediate 1 1.822503 -0.2371038

> head(dat_recognition)
Participant Item Treatment Test Time Score Proficiency Aptitude
1 1 1 Glossing Recall Immediate 1 1.822503 -0.2371038
2 1 18 Writing Recall Immediate 1 1.822503 -0.2371038
3 1 26 Writing Recall Delayed 1 1.822503 -0.2371038
4 1 14 Glossing Recall Immediate 0 1.822503 -0.2371038
5 1 5 Glossing Recall Delayed 0 1.822503 -0.2371038
6 1 10 Glossing Recall Immediate 1 1.822503 -0.2371038