モデル選択できないGLMM

GLMMで解析する際、lme4::lmer()（あるいはlme4::glmer()）を使うことは多い。しかし、これらの関数はデフォルト設定ではモデル選択できない。なぜなら、デフォルトでREML = TRUE（制限付き最尤法「Restricted Maximum Likelihood」）となっているため。

制限付き最尤法（REML=TRUE）はランダム効果の分散推定を不偏にすることを目的とした推定方法のため、固定効果を入れ替えるようなモデル間比較には使えない。推定方法を最尤法（REML=FALSE）と明示せずにモデル選択すると、信頼できない結果が返ってくる。ちなみに、影響は極めて大きい。

pacman::p_load(tidyverse,
               lme4,
               MuMIn)

options(na.action = "na.fail")

## produce dummy data
set.seed(123)
sleepstudy_aug <- sleepstudy %>%
  mutate(
    Age = rnorm(n(), 
                mean = 35, 
                sd = 10),
    Caffeine = sample(c("Low", "Medium", "High"),
                      n(),
                      replace = TRUE),
    Exercise = rnorm(n(), 
                     mean = 3, 
                     sd = 1)
  )

# REML = TRUE by default; Restricted Maximum Likelihood
m0 <- lmer(Reaction ~ Days + Age + Caffeine + Exercise + (1 | Subject),
           data = sleepstudy_aug)

# REML = FALSE; Maximum Likelihood
m1 <- lmer(Reaction ~ Days + Age + Caffeine + Exercise + (1 | Subject),
           data = sleepstudy_aug,
           REML = FALSE)

それぞれモデル選択すると（当然だが）結果は全く異なり、結論に直接影響する可能性が極めて高い。

REML = TRUEとした場合（不適切）¹：

(ms0 <- dredge(m0,
               rank = "AIC") %>% 
    filter(delta < 4))

## Warning in dredge(m0, rank = "AIC"): comparing models fitted by REML

## Fixed term is "(Intercept)"

## Global model call: lmer(formula = Reaction ~ Days + Age + Caffeine + Exercise + 
##     (1 | Subject), data = sleepstudy_aug)
## ---
## Model selection table 
##      (Intrc)     Age Caffn  Days Exrcs df   logLik    AIC delta weight
## [1,]   237.1             + 10.43 4.969  7 -883.375 1780.7  0.00  0.683
## [2,]   245.9 -0.2443     + 10.45 4.866  8 -883.364 1782.7  1.98  0.254
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Subject

REML = FALSEとした場合（適切）：

(ms1 <- dredge(m1,
               rank = "AIC") %>% 
    filter(delta < 4))

## Fixed term is "(Intercept)"

## Global model call: lmer(formula = Reaction ~ Days + Age + Caffeine + Exercise + 
##     (1 | Subject), data = sleepstudy_aug, REML = FALSE)
## ---
## Model selection table 
##      (Intrc)     Age Caffn  Days Exrcs df   logLik    AIC delta weight
## [1,]   238.0               10.35 4.587  5 -895.025 1800.0  0.00  0.357
## [2,]   246.5 -0.2353       10.36 4.488  6 -894.585 1801.2  1.12  0.204
## [3,]   251.4               10.47        4 -897.039 1802.1  2.03  0.130
## [4,]   237.1             + 10.43 4.962  7 -894.178 1802.4  2.31  0.113
## [5,]   260.4 -0.2582       10.47        5 -896.523 1803.0  3.00  0.080
## [6,]   245.9 -0.2438     + 10.45 4.861  8 -893.701 1803.4  3.35  0.067
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Subject

REML = TRUEの場合、Caffeineがすべてのトップモデル（\(\Delta AIC < 2.0\)）に含まれるのに対し、REML = FALSEとすると全く含まれなくなる。実際の論文で起きた場合、結論が真逆になる可能性すらある。怖い。

複数マッチのエラーがでない

RのデフォルトデータであるirisにはSpecies列にsetosa, versicolor, virginica の３種が含まれる。その中の２種、setosa, versicolorをデータフレームから抜き出したいとする。このような複数マッチを指定したい場合、%in%を使う。

df_sub <- iris %>% 
  filter(Species %in% c("setosa",
                        "versicolor"))

table(df_sub$Species)

## 
##     setosa versicolor  virginica 
##         50         50          0

もともと、各種50点の観察値が含まれるデータである。Speciesの頻度をとるとsetosa, versicolorが50となっており、正しく機能しているとわかる。

しかし、ありがちなミスとして、%in%を使うべきところに==（単一の値・文字列との符合を指示する際の演算子）を使ってしまうことがある。このミスは全く異なる結果を返すにも関わらず、エラーおよび警告メッセージが一切でない。試してみる。

df_sub <- iris %>% 
  filter(Species == c("setosa",
                      "versicolor"))

table(df_sub$Species)

## 
##     setosa versicolor  virginica 
##         25         25          0

なぜか、setosa, versicolor からそれぞれ25行ずつ選ばれている。本来50行ずつ選ばれるはず。この25行がどのような基準で選ばれているのかはわからない。このエラーに気付かないと、知らぬ間に不適切なデータの間引きが起き、意味のないデータ解析を進めてしまう。怖い。

“sample” surprise

シミュレーションなどをしていると気が付く、sample関数の（意味の分からない）挙動。あるベクターからいくつかの要素をランダムに取り出したい場合、sample関数のお世話になることが多い。例えば：

# random five draws from 1 to 10
sample(1:10, size = 5)

## [1] 9 5 7 1 2

この挙動から、最初の引数には「くじ」（この場合、1から10までの数字）を作ってやればよい、と理解する。なので、sample(10, size = 1)とすれば、かならず10が返ってくるはず、と思っていると痛い目にあう。

# random one draw from 10...?
set.seed(123)
sample(10, size = 1)

## [1] 3

なぜに3！？となる。

これは、最初の引数が単一値の場合、1:10を与えた場合と同じ挙動をするためのようだ。これはsample関数のHelpページに”Surprise”として紹介されるとともに、より安全な関数の定義が提案されている。

# safer function, return 10 every time
resample <- function(x, ...) x[sample.int(length(x), ...)]
resample(10, size = 1)

## [1] 10

最近、自分でシミュレーションするときは、必ずresample関数を定義してからパッケージを組むようにしている。怖い。