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

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

2019年01月04日

StanやJAGSの計算結果からパラメータのMAP推定値とHDIをまとめて得るためのRの自作関数

 こんにちは。mutopsyです。この記事では,StanやJAGSのfitオブジェクトを渡すとパラメータのMAP推定値とHDI (最高密度区間) をまとめて返してくれる自作関数を紹介します。

 StanやJAGSでパラメータの推定結果の要約を出力すると,通常は事後分布のEAP (i.e., 平均) や分位点の値が出力されますが,事後分布のMAP推定値やHDIが知りたい場合もありますよね。そんなわけで,そういう関数を作ってみました (もう既にあるかもしれないけど知ったこっちゃない)。自作関数といっても,HDIの計算にはHDIntervalというパッケージを利用しますが。関数は以下の通り。

Rコード01: MCMCサンプルからMAP推定値とHDIを計算する関数
library(HDInterval)
maphdi <- function(fit,credMass=0.99,include=TRUE,pars=NA){
  if(credMass<0|credMass>1) print("警告:credMassの値が正しくありません。")
  Class <- class(fit)[1]
  if(Class=="stanfit"){
    e <- as.data.frame(rstan::extract(fit))
  }else if(Class=="rjags"){
    e <- as.data.frame(fit$BUGSoutput$sims.list)
  }else if(Class=="data.frame"){
    e <- fit
  }else{
    print("エラー:未対応なオブジェクトです。")
    stop()
  }

  if(is.na(pars[1])) pars <- colnames(e)
  if (include==TRUE){
    e <- e[,pars]
  }else{
    e <- e[,!colnames(e) %in% pars]
  }
  if(length(pars)==1){
    MAP <- density(e)$x[which.max(density(e)$y)]
    HDI95 <- hdi(e)
    HDIAny <- hdi(e,credMass = credMass)
    Result <- data.frame(MAP = MAP, Lower95 = HDI95[1], Upper95 = HDI95[2],
                         LowerAny = HDIAny[1],UpperAny = HDIAny[2])
    rownames(Result) <- pars
  }else{
    MAP <- apply(e,2,function(z) density(z)$x[which.max(density(z)$y)])
    HDI95 <- apply(e,2,HDInterval::hdi)
    HDIAny <- apply(e,2,function(x) HDInterval::hdi(x,credMass = credMass))
    Result <- data.frame(MAP = MAP, Lower95 = HDI95[1,], Upper95 = HDI95[2,],
                         LowerAny = HDIAny[1,],UpperAny = HDIAny[2,])
  }

  colnames(Result)[4:5] <- c(
    paste("Lower",round(credMass*100,2),sep=""),
    paste("Upper",round(credMass*100,2),sep="")
  )
  return(Result)
  }

 試しにRからStanのfitオブジェクトを渡してみます。データとモデルは何でもよいので適当なものを渡します。

Rコード02: 関数を使ってみる
>fit.eg <- sampling(stanmodel, data=data, seed=123,
                 chains=4, iter=10500, warmup=500, thin=1)
> maphdi(fit.eg)
              MAP    Lower95    Upper95    Lower99    Upper99
beta.1  1.8502036  1.3540767  2.3442492  1.1987086  2.5187496
beta.2  0.6478270  0.5188389  0.7805164  0.4767612  0.8295815
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 こんな風に,Stanのfitオブジェクトをmaphdi()に引数として渡すだけで各パラメータのMAP推定値と95% HDI・99% HDIを計算して出力することができます。Stan2.18, rstan2.18.2, Rtools3.5, R3.5.2, Windows10で動作確認済み。JAGSにも対応しています (JAGS4.3.0, R2jags0.5-7で確認済み)。

 また,fitオブジェクトの代わりに,MCMCサンプルが格納されたデータフレームを渡してもMAP推定値とHDIを返してくれます。

Rコード03: MCMCサンプルを直接渡す
> e <- as.data.frame(rstan::extract(fit.eg))
> head(e)
    beta.1    beta.2    beta.3     beta.4     sigma     lp__
1 2.194799 0.5785987 0.6255656 -0.4086454 0.3094092 97.46736
2 1.726557 0.6609587 0.7787883 -0.7231942 0.3058624 97.31513
3 2.053161 0.6021407 0.6534262 -0.3992853 0.3425830 96.40984
4 1.941408 0.6503287 0.6399921 -0.4178822 0.2976153 97.93040
5 2.518750 0.4448252 0.7676119 -0.7428853 0.3209268 91.90027
6 1.758346 0.6890341 0.6522970 -0.4046868 0.2985565 97.78386
> maphdi(e)
              MAP    Lower95    Upper95    Lower99    Upper99
beta.1  1.8502036  1.3540767  2.3442492  1.1987086  2.5187496
beta.2  0.6478270  0.5188389  0.7805164  0.4767612  0.8295815
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 引数のcredMassに0〜1の範囲の値を与える事で,99% HDIの代わりに任意の確度のHDIを出力することもできます。95% HDIは必ず出力されるようになっています。

Rコード04: 任意のHDIを出力
> maphdi(e,credMass = 0.50)
              MAP    Lower95    Upper95    Lower50    Upper50
beta.1  1.8502036  1.3540767  2.3442492  1.7011924  2.0385575
beta.2  0.6478270  0.5188389  0.7805164  0.6060299  0.6964437
beta.3  0.7103451  0.5996639  0.8201469  0.6729364  0.7490473
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.6429855 -0.4713128
sigma   0.3135227  0.2814639  0.3548575  0.3030061  0.3280369
lp__   97.7211912 93.6796922 99.1662181 96.7857008 98.5912001

 引数のparsにパラメータの名前を与えれば,特定のパラメータの結果のみを出力することもできます。また,include = FALSE とすると,指定したパラメータ以外の結果を出力することもできます。

Rコード05: 一部のパラメータのみ計算

> maphdi(fit.eg,pars=c("beta.1","beta.2"))
            MAP   Lower95   Upper95   Lower99   Upper99
beta.1 1.850204 1.3540767 2.3442492 1.1987086 2.5187496
beta.2 0.647827 0.5188389 0.7805164 0.4767612 0.8295815
> maphdi(fit.eg,pars=c("beta.3"))
             MAP   Lower95   Upper95   Lower99   Upper99
beta.3 0.7103451 0.5996639 0.8201469 0.5627496 0.8524817
> maphdi(fit.eg,include=FALSE,pars=c("beta.1","beta.2"))
              MAP    Lower95    Upper95    Lower99    Upper99
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 以上です。Enjoy!

posted by mutopsy at 19:10 | Rに関するTips

2019年01月01日

既知の確率質量関数からMAP推定値と最高密度区間 (HDI) を計算するためのRの自作関数

 あけましておめでとうございます。mutopsyです。年末にノロウイルスらしきものに感染してしまいましたがぶじ元気に2019年を迎えることができました。今年も何卒よろしくお願いいたします。新年早々ではありますが,必要に駆られてRの自作関数を作ったので共有します。

 既に分かっている確率質量関数を使ってMAP (最大事後確率) 推定値とHDI (highest density interval; 最高密度区間) を求める関数です。実現値として離散値kを返す既知の確率質量関数 (e.g., 二項分布) で事後分布が表現できる場合に,そのような事後分布のMAPとHDIを計算するのに使えます。ちなみに連続値を返す確率密度関数の場合でも,離散化すればこの関数を使って近似的にMAP推定値とHDIを求めることができるはずです。

 2019年1月2日追記:MAPとなる値やHDIの範囲が2つ以上ある場合 (例えば2峰分布など) にも対応できるように関数を修正しました。結果の出力の仕方も微修正。

Rコード01: 既知の確率質量関数からMAP推定値とHDIを計算する関数
hdimass <- function(prob.k,start.k=0,credMass=0.95){
  #警告メッセージ
  if(sum(prob.k)>1) print("警告:確率の和が1を越えています。")
  if(credMass<0|credMass>1) print("警告:credMassの値が正しくありません。")
  
  #HDIの計算の下準備
  credMass1m <- 1 - credMass
  accum<-1:length(prob.k)
  for(k in 1:length(prob.k)){
    if(k ==1){
      accum[k] <- sort(prob.k)[k]
    }else{
      accum[k] <- accum[k-1] + sort(prob.k)[k]
    }
  }
  Interval <- grep(TRUE,prob.k>max(sort(prob.k)[accum<=(min(abs(accum-credMass1m))+credMass1m)]))
  
  #HDIが複数ある場合やHDIが存在しない場合 (全範囲がHDI) を考慮
  if (length(Interval) <= 1){
    Check <- 0
    Interval <- 1:length(prob.k)
  } else{
    Check <- Interval[2:length(Interval)]-Interval[1:(length(Interval)-1)]-1
  }
  if(sum(Check)==0){
    Lower <- min(Interval) + (start.k-1)
    Upper <- max(Interval) + (start.k-1)
  } else{
    Lower <- c(min(Interval),Interval[grep(TRUE,Check!=0)+1]) + (start.k-1)
    Upper <- c(Interval[grep(TRUE,Check!=0)],max(Interval)) + (start.k-1)
  }
  
  #実際にその範囲の値となる確率
  ActualProb <- sum(prob.k[Interval])
  
  #MAPの計算
  MAP <- grep(TRUE,prob.k == max(prob.k)) + (start.k-1)
  
  #結果をリストにまとめて返す
  Result <- list(MAP=MAP,Lower=Lower,Upper=Upper,ActualProb=ActualProb)
  return(Result)
}

 第1引数 (prob.k) としてそれぞれのkが実現する確率のベクトルを,第2引数 (start.k) としてkの最小値 (デフォルトはstart.k=0) を与えることで,MAPと95% HDIの上限・下限を格納したリストが返り値として得られます。また,実際にHDIの範囲内の値をとる確率も返してくれます (離散分布なので95% HDIであっても厳密には95%にはなりませんが,95%に最も近い値になるはずです)。第1引数 (prob.k) には現実的な長さのベクトルを与えれば十分です。例えば N = 20,θ = 0.5 の二項分布であれば,長さが20もあれば総和がほぼ1になるので計算上は問題ありません (以下の例では一貫して長さ500のベクトルを与えていますがもっと少なくても結果は変わりません)。ちなみに第3引数であるcreadMassに(0,1)の範囲の値を与えれば (例えばcredMass = 0.99),異なる確度のHDIを出力することもできます (99% HDIなど)。

 計算の中身を簡単に説明します。まず最初にfor文で確率質量 (prob.k)を昇順にソートし,順番に累積確率を計算しています。次に,累積確率が閾値 (95% HDIなら5%)に最も近い点を探し,各々のkがHDIに含まれるか否かを論理値ベクトルで表現します。この論理値ベクトルを使って,(credMass)% HDIに含まれるkの値を抽出し,Intervalという内部の変数に格納しています。HDIが1つしかない場合は,Intervalの最小値と最大値ををHDIの下限・上限とみなしています。HDIが2つ以上ある場合は,Interval内で数値が飛び飛びになっている部分を見つけてその変化点をHDIの境界とみなしています。MAP推定値は最大事後確率を実現するkの値を求めることで計算しています。

 実際にこの関数を使って,さまざまなパラメータを持つ二項分布のMAPとHDIを計算してみました。

Rコード02: パラメータと引数を変えながら二項分布で試してみた結果
> #(1) 成功確率0.5,試行回数20の二項分布 (95% HDI)
> hdimass(dbinom(0:500,20,0.5),start.k = 0)
$MAP
[1] 10

$Lower
[1] 6

$Upper
[1] 14

$ActualProb
[1] 0.9586105

> #(2) 成功確率0.2,試行回数20の二項分布 (95% HDI)
> hdimass(dbinom(0:500,20,0.2),start.k = 0)
$MAP
[1] 4

$Lower
[1] 1

$Upper
[1] 7

$ActualProb
[1] 0.9563281

> #(3) 成功確率0.6,試行回数100の二項分布 (95% HDI)
> hdimass(dbinom(0:500,100,0.6),start.k = 0)
$MAP
[1] 60

$Lower
[1] 51

$Upper
[1] 69

$ActualProb
[1] 0.948118

> #(4) 成功確率0.6,試行回数100の二項分布 (99% HDI)
> hdimass(dbinom(0:500,100,0.6),start.k = 0,credMass = 0.99)
$MAP
[1] 60

$Lower
[1] 48

$Upper
[1] 72

$ActualProb
[1] 0.9896389

> #(5) 成功確率0.5,試行回数20,k>=7の切断二項分布 (95% HDI)
> hdimass(dbinom(7:500,20,0.5)/sum(dbinom(7:500,20,0.5)),start.k = 7,credMass = 0.95)
$MAP
[1] 10

$Lower
[1] 7

$Upper
[1] 13

$ActualProb
[1] 0.9388129

> #(6) MAPとHDIが2つずつある混合二項分布 (95% HDI)
> hdimass(dbinom(0:500,100,0.3) * 0.5 + dbinom(0:500,100,0.7) * 0.5)
$MAP
[1] 30 70

$Lower
[1] 22 61

$Upper
[1] 39 78

$ActualProb
[1] 0.9501802
#↑95% HDIは[22,39]または[61,78]の範囲

 いい感じに計算できていると思います。

 最後に,この関数を作った経緯を少しだけ説明します。僕が今行っている分析で離散パラメータの推定を行っているのですが,MCMCの結果からMAP推定値とHDIを十分高い精度で計算するためには発生させるMCMCサンプルの数 (JAGSのn.iterの値) をかなり増やさなければならないことが分かりました。これは現実的ではなかったので,数値積分を利用して事後分布を近似的に計算し,その事後分布からMAP推定値とHDIを計算することにしました。その副産物として今回の関数ができたという訳です。Rを使って事後分布を近似的に求める方法についてもいつか記事で解説するかもしれません。

 という訳で,需要があるかどうかも怪しい誰得関数の紹介から2019年が始まりました。Enjoy!

posted by mutopsy at 16:52 | Rに関するTips