変量効果モデルとAICと目で見る尤度関数の学習記録

【学習動機】

前回, AICについて少し書いた.

moratoriamuo.hatenablog.com

その続きで, 線形混合効果モデル(LMM)で変量効果(ランダム効果)をモデルに入れたときって, AICが使えるのか, つまりは, モデルの尤度関数が正規分布で近似できそうなのか否かということが気になって, Rで人工データ作って実験してみることにした.

【学習内容】

概要を書く. グループごとの変量効果と誤差項の和を真の分布として, 人工データを生成する. 3つのモデルをあてはめる.

1. 切片のみで誤差項を正規分布とした線形モデル(つまりは, ただの正規分布)

2. グループごとの固定効果と誤差項を正規分布とした線形モデル

3. グループごとに正規分布からの変量効果と誤差項を正規分布とした線形モデル(真の分布と同じ構造)

推定されたパラメータを確認し, AICも確認する.

次に, 同じようにデータを1000回発生させて, それぞれのモデルのAICを確認する.

最後に, 3つ目のモデル, 変量効果モデルの尤度関数の形状を確認する.

では, 人工データの説明をする.

f:id:moratoriamuo271:20191012122507p:plain

グループ数がM, 各グループにそれぞれ3ケースのバランスしたデータを生成する. biがグループの変量効果, εiがそれらで捉えられない誤差項.

使用するパッケージと乱数の種指定.

library(tidyverse)
library(ggplot2)
library(lme4)
library(lmerTest)
set.seed(271)

真の分布のパラメータの設定.
σ=2, σb=2

sigma <- 2
sigma_b <- 2

グループ数は12. それぞれ3ケース. グループ番号もつけておく.

M <- 12
Mi <- rep(3,M)
gindex <- rep(x = 1:M,each=3)

データの生成.

b <- numeric(M)
e <- numeric(M)
y <- numeric(0)

for(i in 1:M){
  b[i] <- rnorm(1,0,sigma_b)
  y <- c(y,b[i]*rep(1,Mi[i]) + rnorm(Mi[i],0,sigma))
}

df <- data.frame(
  y = y,
  x1 = rep(1,sum(Mi)),
  g = gindex %>% as.factor()
)

yのヒストグラムのチェック.

hist(df$y)

ggplot(df) + geom_histogram(aes(y,fill=g))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

切片のみの単純な線形モデルと, 切片にグループの固定効果をつけた線形モデルのあてはめ.
あてはめた後, 残差をグループごとに色付けしてチェック.

model1 <- lm(y~1,data=df)
summary(model1)

## 
## Call:
## lm(formula = y ~ 1, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5753 -1.9004 -0.1703  2.3216  6.1391 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.1077     0.4769  -0.226    0.823
## 
## Residual standard error: 2.861 on 35 degrees of freedom

ggplot(df %>% mutate(res=residuals(model1))) + geom_histogram(aes(res,fill=g))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

model2 <- lm(y~0+g,data=df)
summary(model2)

## 
## Call:
## lm(formula = y ~ 0 + g, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0633 -0.9106 -0.3025  0.8729  3.3338 
## 
## Coefficients:
##     Estimate Std. Error t value Pr(>|t|)    
## g1   -1.6948     0.9199  -1.842 0.077811 .  
## g2    4.9049     0.9199   5.332  1.8e-05 ***
## g3    1.0250     0.9199   1.114 0.276207    
## g4   -3.8096     0.9199  -4.141 0.000368 ***
## g5   -3.4402     0.9199  -3.740 0.001014 ** 
## g6    1.6255     0.9199   1.767 0.089932 .  
## g7    1.4748     0.9199   1.603 0.121978    
## g8   -2.6787     0.9199  -2.912 0.007643 ** 
## g9   -0.5467     0.9199  -0.594 0.557850    
## g10   1.3370     0.9199   1.453 0.159067    
## g11  -1.5896     0.9199  -1.728 0.096826 .  
## g12   2.0996     0.9199   2.282 0.031624 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.593 on 24 degrees of freedom
## Multiple R-squared:  0.7877, Adjusted R-squared:  0.6815 
## F-statistic:  7.42 on 12 and 24 DF,  p-value: 1.751e-05

ggplot(df %>% mutate(res=residuals(model2))) + geom_histogram(aes(res,fill=g))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

混合効果モデルを推定. 最尤推定とREMLで. 今回は固定効果がないので分散成分の推定が一緒になる.

model3 <- lmer(y~0+(1|g),data=df,REML = FALSE)
summary(model3)

## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: y ~ 0 + (1 | g)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##    163.8    166.9    -79.9    159.8       34 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3706 -0.6620 -0.2384  0.4543  2.1791 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  g        (Intercept) 5.432    2.331   
##  Residual             2.539    1.593   
## Number of obs: 36, groups:  g, 12

model3_REML <- lmer(y~0+(1|g),data=df,REML = TRUE)
summary(model3_REML)

## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: y ~ 0 + (1 | g)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##    163.8    166.9    -79.9    159.8       34 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3706 -0.6620 -0.2384  0.4543  2.1791 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  g        (Intercept) 5.432    2.331   
##  Residual             2.539    1.593   
## Number of obs: 36, groups:  g, 12

AIC(model1,model2,model3)

##        df      AIC
## model1  2 180.8405
## model2 13 147.1056
## model3  2 163.7519

AICを見るとモデル2の固定効果のモデルが一番小さく出ている. 本来であれば, モデル3が選ばれてほしいところだが, これは偶然なのか違うのか. 検証するために, 1000回発生させてAICを算出してみる.

3つのモデルのAICのヒストグラムを描く

vAIC <- matrix(0, nrow = 1000, ncol = 3)

for(j in 1:1000){
  Mi <- rep(3,M)
  gindex <- rep(x = 1:M,each=3)
  
  b <- numeric(M)
  e <- numeric(M)
  y <- numeric(0)
  
  for(i in 1:M){
    b[i] <- rnorm(1,0,sigma_b)
    y <- c(y,b[i]*rep(1,Mi[i]) + rnorm(Mi[i],0,sigma))
  }
  
  df <- data.frame(
    y = y,
    x1 = rep(1,sum(Mi)),
    g = gindex %>% as.factor()
  )
  
  model1 <- lm(y~1,data=df)
  model2 <- lm(y~0+g,data=df)
  model3 <- lmer(y~0+(1|g),data=df,REML = FALSE)
  vAIC[j,] <- AIC(model1,model2,model3)$AIC
}

## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular

dfAIC <-data.frame(
  AIC1 = vAIC[,1],
  AIC2 = vAIC[,2],
  AIC3 = vAIC[,3]
)

ggplot(data=dfAIC) + geom_histogram(aes(x=AIC1),alpha=0.3,fill="red") +
  geom_histogram(aes(x=AIC2),alpha=0.3,fill="blue") +
  geom_histogram(aes(x=AIC3),alpha=0.3,fill="green") +
  theme_bw() + labs(x="AIC")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

boundary (singular) fit: see ?isSingular とメッセージが何度か出ていて, lmerの推定が上手くいってないところもある気がするが今回は無視.. [https://www.researchgate.net/post/What_does_singular_fit_mean_in_Mixed_Models]

モデル3とモデル2のAICの差をヒストグラムで描き, 1000個中何個がモデル3を選ぶのかをチェックする.

ggplot(data=dfAIC) + geom_histogram(aes(x=AIC3-AIC2),colour="grey") +
  geom_vline(xintercept = 0, colour="red") +theme_bw() +
  geom_text(aes(x=-7,y=80,label=sum(AIC3<AIC2))) + geom_text(aes(x=-4,y=80,label="/1000"))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

どうもモデル2の方が選ばれやすいみたいだ. この原因が何なのか考えていて, 一つそれっぽいと思ったのが, 渡辺澄夫先生の『階層ベイズ法とWAIC』[http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/h_b.html]に書いてあるようなことだ. つまり, 何を予測しているかが, モデル2とモデル3では違っていてAICが比較できるものになっていないのかもしれない. これについては, まだ自分の中で理解と整理ができていないので, stanの学習を進めつつ, 分かり次第, また記事として書いていくつもりだ.

原因として, もう一つ考えられたのは, モデル3の尤度関数が正規分布で近似できそうにないからAICの適用が上手くいっていないという可能性だ. これを調べるため, 以下ではモデル3の尤度関数の形を描くことにする. (というか, これをやりたいがためにパラメータを2つに抑えていた. それより多いと可視化ができないから.)

尤度関数の描き方については, 以下の記事を参考にした.

abrahamcow.hatenablog.com

この記事へのコメントをここに書いてしまいますが, 混合正規分布の乱数発生の関数に記載ミス?があると思われます. 記事記載のコードでは, 半分半分で各正規分布からの乱数が発生しているだけになっています.

rmixnorm <- function(n,a,b){
    n1 <- sum(runif(n) < a)
    c(rnorm(n-n1,0,1),rnorm(n1,b,1))
}

などとすると上手くいくはずです.

さて, 尤度関数を使える形になるよう導出する.

f:id:moratoriamuo271:20191012122703p:plain

# par = c(sigma,sigma_b)
# df = cbind.data.frame(y,g)
dlmmnorm <- function(par,df,log=TRUE){
  a <- par[1]
  b <- par[2]
  N <- length(df$y)
  gs <- unique(df$g)
  ngs <- length(gs)
  element_sum <- numeric(ngs)
  for(i in gs){
    dfi <- df %>% filter(g==i) %>% mutate(y2=y^2)
    ni <- length(dfi$y)
    A <- ni*b^2 + a^2
    B <- sum(dfi$y2)
    C <- sum(dfi$y)^2
    element_sum[i] <- log(1 / sqrt(2*pi*a^2)^ni * sqrt(a^2/A)) - B/(2*a^2) + b^2*C/(2*A*a^2)
  }
  ll <- sum(element_sum)
  if(log){
    ll
  }else{
    exp(ll)
  }
}

op <- optim(c(2.5,0.5),dlmmnorm,control = list(fnscale=-1), df=df, method = "BFGS",
            hessian = TRUE)
op

## $par
## [1] 1.724198 2.580073
## 
## $value
## [1] -82.95735
## 
## $counts
## function gradient 
##       32       12 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
##             [,1]       [,2]
## [1,] -16.2874199 -0.6085378
## [2,]  -0.6085378 -2.7317524

#a <- seq(0.1,4.1,length.out = 80)
#b <- seq(0.1,4.1,length.out = 80)

a <- seq(-3.9,4.1,length.out = 80)
b <- seq(-3.9,4.1,length.out = 80)


parms <-expand.grid(a,b)
L1 <- apply(parms,1,dlmmnorm,df=df,log=FALSE)
datL1 <- as.data.frame(parms) %>% 
  set_names(c("a","b")) %>% 
  mutate(L=L1)

ggplot(datL1,aes(x=a,y=b))+
  geom_contour(aes(z=L),size=1)+
  theme_bw()

最適化の計算を回していて何度か負の値に落ち込んだので, 負の方向まで確認してみたら, 対称性をもって尤度関数の峰ができていた. 標準偏差の尤度関数としてグラフを書いているから, そりゃあそうかと納得. 正の領域で制約をかければ, 尤度関数自体は正規分布で近似できそうな気がするという結論で良さそう.

ということは, AICのモデル選択の問題点は先に述べた『階層ベイズ法とWAIC』に書かれているようなことに起因しているということなのかもしれない.

【追記：2019/10/19】

長島健悟先生が線形混合効果モデルの同時・周辺尤度関数の導出をメモしているのを発見した. 計算する人の参考になるかもしれないので, リンクを貼っておく.

https://nshi.jp/contents/stat/lme/

【学習予定】

ベイズ推測とWAICの理解を深めていきたい. stanをもっと気軽に使えるように練習していきたいところ.

霞と側杖を食らう

ほしいものです。なにかいただけるとしあわせです。[https://www.amazon.jp/hz/wishlist/ls/2EIEFTV4IKSIJ?ref_=wl_share]

変量効果モデルとAICと目で見る尤度関数の学習記録

【学習動機】

【学習内容】

【学習予定】