ベイズマックスへようこそ。
統計に関するあれやこれやを自由気ままに投稿します。
全記事一覧はこちらから。

2020年09月14日

混合効果モデルにおいて個々の変量効果の点推定値から計算した分散が真値よりも過小推定される理由

 線形混合モデルもしくは一般化線形混合モデル (GLMM) を使って変量効果を推定するときに,REML(制限付き最尤推定法)で推定された変量効果の分散と,個々の変量効果の点推定値を使って計算した分散が一致しない理由についての備忘録です。Rでいえば,summary(model)やVarCorr(model)で確認できる変量効果の分散の大きさと,var(ranef(model))で計算した変量効果の分散の大きさが一致しない理由についての解説です。

 まずは具体例から。Rで架空のデータを作って線形混合モデルの当てはめをし,その推定結果を確認してみます。説明変数がひとつだけのシンプルなモデルで,切片と傾きのそれぞれについて固定効果と変量効果を考えます。真値として,切片の固定効果を0.5,傾きの固定効果を1.0とし,変量効果に関しては切片の分散を1^2 (=1),傾きの分散を1.5^2(=2.25)と指定しました。

線形混合モデルの当てはめ
library(tidyverse)
library(MuMIn)
library(lmerTest)
set.seed(8823)

## 真値の設定 ----

beta_fixed <- c(0.5,1.0) #固定効果の真値
sigma_beta_rand <- c(1.0,1.5) #変量効果の標準偏差
sigma_err <- 5.0 #残差の標準偏差

## データの作成 ----

N_id<-500 #グループの数(e.g., 実験参加者の人数)
N_trial <- runif(N_id,20,30) %>% round(0) #各グループの値の数(e.g., 参加者ごとの試行数)

N <- sum(N_trial) #データの総数

ID <- NULL
for(i in 1:N_id) ID <- c(ID,rep(i,N_trial[i]))

beta0_i <- rnorm(N_id,beta_fixed[1],sigma_beta_rand[1])
beta1_i <- rnorm(N_id,beta_fixed[2],sigma_beta_rand[2])

X <- rnorm(N,0,1) #説明変数
Err <- rnorm(N,0,sigma_err) #残差

Y <- beta0_i[ID] + beta1_i[ID] * X1 + Err  #観測変数を作る

## 線形混合モデルを当てはめてパラメータを推定する

model <- lmer(Y ~ 1 + X + (1+X|ID))

summary(model)
# > summary(model)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
# Formula: Y ~ 1 + X + (1 + X | ID)
# 
# REML criterion at convergence: 76647.6
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -3.8479 -0.6567  0.0041  0.6601  3.9094 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  ID       (Intercept)  1.079   1.039         
#           X            1.829   1.352    -0.13
#  Residual             25.098   5.010         
# Number of obs: 12504, groups:  ID, 500
# 
# Fixed effects:
#              Estimate Std. Error        df t value Pr(>|t|)    
# (Intercept)   0.47940    0.06504 495.92464   7.371 7.15e-13 ***
# X             0.95096    0.07657 488.35016  12.419  < 2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Correlation of Fixed Effects:
#   (Intr)
# X -0.067

VarCorr(model)
# > VarCorr(model)
# Groups   Name        Std.Dev. Corr  
# ID       (Intercept) 1.0386         
#          X           1.3524   -0.128
# Residual             5.0098 

 この通り,ちゃんと真値に近い推定値が得られました。特に変量効果の分散に関していえば,切片は1.04^2 (真値は1.0^2),傾きは1.35^2(真値は1.5^2)と,多少のずれはありますがそれなりにちゃんと推定できているようです。ところで,ranef()という関数を使うとグループごと(e.g., 実験参加者ごと)の変量効果の推定値を確認することができます。心理学であれば,例えば実験参加者ごとの変量効果の大きさが別の関心のある変数と相関するかどうかを調べたい,といった場合があり得ると思いますが,そういうときにはranef()を使うと便利です。ranef()で出力されるのは正規分布から得られた個々の変量効果の実現値なので,直観的には,それらの値から計算した分散は,先ほど確認した変量効果の分散と一致しそうです。確認してみましょう。

ranef(model)の分散を確認してみる
ranef(model)$ID %>% head() #ranef()の中身を確認する。
# > ranef(model)$ID %>% head()
#   (Intercept)          X
# 1  -0.1284018  1.5457369
# 2   0.1054594  1.8687398
# 3  -0.8152035 -1.2718659
# 4   0.4388866 -0.3965585
# 5   0.9940244  0.4130416
# 6   0.2910409  0.2689888

ranef(model)$ID %>% apply(2,var) #変量効果の分散を計算する。
# > ranef(model)$ID %>% apply(2,var)
# (Intercept)           X 
#  0.5532284   1.1435624 #真値は1.0と2.25
ranef(model)$ID %>% apply(2,sd) #変量効果の標準偏差を計算する。
#> ranef(model)$ID %>% apply(2,sd) #変量効果の標準偏差を計算する。
# (Intercept)           X 
#   0.7437933   1.0693748 #真値は1.0と1.5

 おやー?やけに分散が小さく推定されてしまいました。真値よりも明らかに小さいし,summary()やVarCorr()で確認した値とも一致しません。真の分布のグラフとranef()から得た値のヒストグラムを重ねてプロットしたものが下の図です。明らかに分散が小さいですね。

hist_ranef.png

 どうして分散が一致しない(過小推定される)のか気になって調べてみたら,StackExchangeに同じ質問を見つけることができました。

 かみ砕いて説明するとこんな感じです:ranef()を使うとグループごと(参加者ごと)の変量効果の点推定値を得ることができますが,当然ながら新たに別の標本からデータを得て同じ分析を行ったら異なる推定値が得られます。このように,個々の変量効果の値は何らかの標本分布に従って変動するので,その変動を加味せずに分散を計算したら真値よりも小さな値になってしまうのは当然ということです。裏を返せば,ちゃんと推定量の標準誤差を考慮して計算すれば真値に近い値が得られるはずです。armというパッケージのse.ranef()という関数を使えば個々の変量効果の標準誤差を計算することができるので,これを利用して標準誤差込みの分散を改めて計算してみましょう。

標準誤差を加味した分散を計算する
library(arm)

se <- se.ranef(model)$ID
head(se) #中身の確認
#   (Intercept)         X
# 1   0.7033685 0.7794396
# 2   0.7110790 0.9662543
# 3   0.7139412 0.6947599
# 4   0.7383144 0.8239593
# 5   0.7010285 0.8233850
# 6   0.6855417 0.7371017

v0 <- mean(se[,1]^2) + var(ranef(model)$ID[,1])
sqrt(v0) #切片の変量効果の標準偏差
#> sqrt(v0)
#[1] 1.038639 #真値は1.0(VarCorrの出力結果1.0386と一致)

v1 <- mean(se[,2]^2) + var(ranef(model)$ID[,2])
sqrt(v1) #傾きの変量効果の標準偏差
# > sqrt(v1)
# [1] 1.352343 #真値は1.5(VarCorrの出力結果1.3524と一致)

 完全に一致!ということで解決です。個々の変量効果の点推定値を使う時にはこのことを知っておいたほうが良いと思います。例えば論文等を書くときにranef()の分散を変量効果の分散として報告してしまうといった誤りを犯さないように気を付けましょう。最尤推定を使って変量効果を推定すると,推定値の背後にある標本分布の存在をついつい忘れてしまいがちですが,ベイズ推定ならMCMCサンプルのばらつきという形で推定の不確かさを自然に確かめることができるので,混合効果モデルはいっそベイズでやってしまったほうが安全かもしれません。推定の正確さという点でもベイズに利がありそうです(特にサンプルサイズが小さいとき)。という訳でEnjoy!

posted by mutopsy at 20:28 | Rに関するTips