Rで一般化線形混合モデルの分析に挑戦！！
年末から少し勉強を始めていて、統計が得意な友だちにも何度も見てもらって、一般化線形混合モデルの分析をやってみるところまではマスターした！と思っていた私ですが、久しぶりに少し余裕が出てきたので、別のデータを分析してみようとしたのですが…全然うまくいきませんでした。
前回教えてもらったRのコードもすべて記録してあるのに…
半日悩んで、いろいろ検索しまくって、仕方ないので友だちにZoomで私のRstudioの画面を共有して、できているところまで説明している途中で気が付きました。
最初にパッケージをインストールするとき、パッケージの名前を” ”で囲むことを忘れていたんです。
そんな初歩的なミスを…とショックでしたが、自分で気が付けたし、良かったとしたところでしょう（笑）
ということで、私はまだ方法しか理解できておらず、分析の中身をふんわりとしか理解できていないので、せっかくZoomで友だちが見てくれているので、いろいろ教えてもらいました(^^)
前回の復習をして、ちゃんとわかっていないところを再確認し、さらに今回は分析する要因が増えた場合のやり方まで教わって、少しパワーアップです！
次回やろうとしたときに、パッケージをインストールできないという失敗をしないよう、しっかりここに書いておきます！笑
（前回挑戦した時より、パソコンが新しくなっているので、パッケージのインストールをするところからやってみる必要があったのです。）
lmer関数は、線形混合モデルや一般化線形混合モデルの推定計算をする関数だそうで、今回私は一般化線形混合モデルの分析に挑戦するのでちょっとだけ勉強しています(^^)
使い方はこんな感じだそうです。
lmer(Y ~ x1 + x2 (ID|factor), data = data)
まずは、パッケージのインストールと読み込みをします。
lme4パッケージでlmerを用いた分析ができるのですが、p値が計算されません。lmerTestパッケージも一緒に読み込んでおくと、p値を一緒に計算してもらえるようになります。
install.packages(”lme4”)#lme4パッケージのインストール
library(lme4)#lme4パッケージの読み込み
install.packages("lmerTest")#lmerTestパッケージのインストール
library(lmerTest)#lmerTestパッケージの読み込み
１．2水準の一般化線形混合モデルの分析
データを読み込みます。
エクセルで必要な表をコピーして、一番上の行をラベルとして読み込んでもらいます。
read.clip()でもクリップボードのデータを読み込むことができますが、これはanova君を読み込んでいないと使えません。今回は、read,tableという関数を使います。
data <- read.table("clipboard", header = T)#クリップボードのデータをdataに読み込んで一番上をラベルにする
読み込んだファイルはこんな感じでした。
IDを見るとわかりますが、参加者01のデータが2行あって、それをダミー変数で区別できるようになっています。
脳波のデータ（P300、N2、LPP ）が目的変数、性格特性の得点（Psychopathy、Machiavellian、Narcissism）が説明変数として、条件（0.5:Positive画像、-0.5:Negative画像）をダミー変数として区別するようにします（性格特性とダミー変数は、中心化してあります）。
そうして、性格特性が条件ごとの脳波データにどのように影響するのかについての分析を行います。
まずは、性格特性（Psychopathy）が条件ごとの脳波データ（P300）にどのように影響するのかについて分析してみます。
fit <- lmer(P300 ~ Psychopathy * (1|ID)+(-1+dummy|ID), data)#一般化線形混合モデルによる分析
summary(fit)#結果を表示させる
結果はこんな風にでてきました。
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: P300 ~ Psychopathy * (1 | ID) + (-1 + dummy | ID)
Data: data
REML criterion at convergence: 139.8
Scaled residuals:
Min 1Q Median 3Q Max
-1.49145 -0.37942 -0.09697 0.24262 1.81294
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 3.2540 1.8039
ID.1 dummy 0.5664 0.7526
Residual 0.7901 0.8889
Number of obs: 36, groups: ID, 18
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.2534 0.4502 16.0004 7.226 2.02e-06 ***
Psychopathy -0.0286 0.0962 16.0004 -0.297 0.77 ←有意差なし
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
Psychopathy 0.000
赤字のところを他の脳波成分や性格特性のデータに入れ替えて、分析をやってみます！！
２．3水準の一般化線形混合モデルの分析
ここまでのやり方は、条件が２つの場合のやり方です。
具体的には、呈示したPositive刺激とネガティブ刺激の差を見ていくものでした。
しかし、実はPositive刺激とNegative刺激だけでなく、Neutral刺激も呈示していたのです。
そうなると、ダミー変数をどう指定したらいいのかわからなかったのですが、やり方を教えてもらいました！
こういう場合は、ダミー変数を２つ用意すると良いそうです。
どんな仮説を設定しているかによって、ダミー変数の作り方が変わってくるようなのですが、今回のデータでは、統制条件であるNeutralがあるので、Neutral刺激とPositive刺激の比較、Neutral刺激とNegative刺激の比較ができるようにダミー変数を作りました。
読み込んだデータはこちら。
さっきのデータとほぼ一緒ですが、IDのところの参加者01が3行になっています。そして、ダミー変数が２つあります。上からPositive、Neutral、Negativeになているので、統制条件であるNeutralが常に０になっていて、比較対象のPositiveが1になるdummy1とNegativeが1になるdummy2があります。
先ほどとほぼ同じですが、脳波のデータ（P300、N2、LPP ）が目的変数、性格特性の得点（Psychopathy、Machiavellian、Narcissism）が説明変数として、条件（Positive、Neutral、Negative）をダミー変数として区別するようにします（性格特性とダミー変数は、中心化してあります）。
そうして、性格特性が3つの条件ごとの脳波データにどのように影響するのかについての分析を行います。
まずは、性格特性（Psychopathy）が3つの条件ごとの脳波データ（P300）にどのように影響するのかについて分析してみます。
fit <- lmer(P300 ~ dummy1 + dummy2 + Psychopathy + dummy1 : Psychopathy + dummy2 : Psychopathy +(1|ID)+(-1+dummy1|ID)+(-1+dummy2|ID), data) > #dummy変数が2つの場合の一般化線形混合モデルによる分析
summary(fit)#結果を表示させる
結果はこんな風にでてきました。
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: P300 ~ dummy1 + dummy2 + Psychopathy + dummy1:Psychopathy + dummy2:Psychopathy +
(1 | ID) + (-1 + dummy1 | ID) + (-1 + dummy2 | ID)
Data: data
REML criterion at convergence: 205.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.86581 -0.37546 -0.04477 0.37409 1.48376
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 2.6391 1.6245
ID.1 dummy1 1.3391 1.1572
ID.2 dummy2 0.0000 0.0000
Residual 0.8866 0.9416
Number of obs: 54, groups: ID, 18
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.43917 0.44257 21.44835 7.771 1.13e-07 ***
dummy1 -0.13250 0.41582 22.12539 -0.319 0.753
dummy2 -0.24078 0.31386 15.34604 -0.767 0.455
Psychopathy -0.06065 0.09456 21.44835 -0.641 0.528
dummy1:Psychopathy 0.08142 0.08885 22.12539 0.916 0.369
dummy2:Psychopathy -0.01752 0.06706 15.34604 -0.261 0.797
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) dummy1 dummy2 Psychp dmm1:P
dummy1 -0.268
dummy2 -0.355 0.377
Psychopathy 0.000 0.000 0.000
dmmy1:Psych 0.000 0.000 0.000 -0.268
dmmy2:Psych 0.000 0.000 0.000 -0.355 0.377
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
どこにも有意差はなかったので、ちょっと差がでていた条件の結果を紹介しますね。
こちらは、サイコパシー得点が3つの条件ごとのN2にどのように影響するのかについての分析を行った結果です。
> fit <- lmer(N2 ~ dummy1 + dummy2 + Psychopathy + dummy1 : Psychopathy + dummy2 : Psychopathy +(1|ID)+(-1+dummy1|ID)+(-1+dummy2|ID), data) > summary(fit)#結果を表示させる Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest'] Formula: N2 ~ dummy1 + dummy2 + Psychopathy + dummy1:Psychopathy + dummy2:Psychopathy + (1 | ID) + (-1 + dummy1 | ID) + (-1 + dummy2 | ID) Data: data REML criterion at convergence: 206.8 Scaled residuals: Min 1Q Median 3Q Max -1.81854 -0.37034 0.00532 0.52450 1.41051 Random effects: Groups Name Variance Std.Dev. ID (Intercept) 1.2283 1.1083 ID.1 dummy1 0.9072 0.9524 ID.2 dummy2 0.9536 0.9765 Residual 1.1235 1.0600 Number of obs: 54, groups: ID, 18 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) -2.27339 0.36146 18.30366 -6.289 5.8e-06 *** dummy1 0.06644 0.41861 10.00213 0.159 0.8770 dummy2 1.09028 0.42168 22.74326 2.586 0.0166 * Psychopathy 0.11356 0.07723 18.30366 1.470 0.1584 dummy1:Psychopathy -0.05828 0.08944 10.00213 -0.652 0.5293 dummy2:Psychopathy -0.09522 0.09010 22.74326 -1.057 0.3017 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) dummy1 dummy2 Psychp dmm1:P dummy1 -0.413 dummy2 -0.410 0.354 Psychopathy 0.000 0.000 0.000 dmmy1:Psych 0.000 0.000 0.000 -0.413 dmmy2:Psych 0.000 0.000 0.000 -0.410 0.354
dummy2のところに有意差が出ています。
dummy2には、NeutralとNegativeの比較をする変数が入っているので、dummy2が有意であるということは、「サイコパシー得点に関わらず、N2の値はNegative刺激に対するものよりもNeutral刺激に対するものの方が大きい」ということが言えます。
３．まとめ
今回使ったコード（スクリプト？）のまとめです。
これをコピーして、左上のペインに貼り付けて、「→Run」のボタンを押すと進んでいきます(^^)
ダミー変数が１つの場合
install.packages(”lme4”)#lme4パッケージのインストール
library(lme4)#lme4パッケージの読み込み
install.packages("lmerTest")#lmerTestパッケージのインストール
library(lmerTest)#lmerTestパッケージの読み込み
data <- read.table("clipboard", header = T)#クリップボードのデータをdataに読み込んで一番上をラベルにする
fit <- lmer(P300 ~ Psychopathy * (1|ID)+(-1+dummy|ID), data)#一般化線形混合モデルによる分析
summary(fit)#結果を表示させる
ダミー変数が２つの場合
install.packages(”lme4”)#lme4パッケージのインストール
library(lme4)#lme4パッケージの読み込み
install.packages("lmerTest")#lmerTestパッケージのインストール
library(lmerTest)#lmerTestパッケージの読み込み
data <- read.table("clipboard", header = T)#クリップボードのデータをdataに読み込んで一番上をラベルにする
fit <- lmer(P300 ~ dummy1 + dummy2 + Psychopathy + dummy1 : Psychopathy + dummy2 : Psychopathy +(1|ID)+(-1+dummy1|ID)+(-1+dummy2|ID), data) > #dummy変数が2つの場合の一般化線形混合モデルによる分析
summary(fit)#結果を表示させる
これで、次回から１人で分析ができる（はず）！！
