お久しぶりの投稿です。

本日のテーマ

過去に書いた自分の記事に突っ込みを入れます。

過去に、lmer/glmerで信頼区間を算出する方法、としてconfint.mermod関数を紹介しました。
https://qiita.com/Jwellyfish/items/4b667a4dd1e0da18ee9f

現時点では (おそらくその当時も) もっと簡単な方法が存在していますので、そのご紹介です。
併せて、効果量 (effect size/risk ratio) と信頼区間をプロットする方法も紹介します。

データはHSAURパッケージのBtheBを使うことにします。

概要

lmer()/glmer()で作成されたmerModオブジェクトにbroom::tidy()を噛ませることで、tibbleとして出力してくれます。summary(lmer/glmer) ではp値および信頼区間を算出してくれませんが、broom::tidy()の引数conf.int = TRUEを設定することで信頼区間を出力してくれます。broom::tidy()で出力した効果量と信頼区間をforestplot::forestplot()に引き渡すことで、簡単に効果量をプロットできます。

準備

まずはBtheBを準備しましょう。
この辺の前処理は、REQUIER02の土屋さんのご発表をtidyverseで再現したものです。

> library(HSAUR)
> data <- BtheB %>%
  tibble::rownames_to_column(var = "ID") # adding rownames as a variable "ID"

> head(data)
  ID drug length treatment bdi.pre bdi.2m bdi.4m bdi.6m bdi.8m
1  1   No    >6m       TAU      29      2      2     NA     NA
2  2  Yes    >6m     BtheB      32     16     24     17     20
3  3  Yes    <6m       TAU      25     20     NA     NA     NA
4  4   No    >6m     BtheB      21     17     16     10      9
5  5  Yes    >6m     BtheB      26     23     NA     NA     NA
6  6  Yes    <6m     BtheB       7      0      0      0      0

BtheBデータは横持ちなので、縦持ちに変換します。

data1 <- data %>%
  gather(key = "month",
         value = "bdi",
         `bdi.2m`,`bdi.4m`,`bdi.6m`,`bdi.8m`
         )

横持時の変数名を時系列変数に変換します。

data1$month <- ifelse(data1$month == "bdi.2m", 2, data1$month)
data1$month <- ifelse(data1$month == "bdi.4m", 4, data1$month)
data1$month <- ifelse(data1$month == "bdi.6m", 6, data1$month)
data1$month <- ifelse(data1$month == "bdi.8m", 8, data1$month)
 
subset(data1, ID %in% c("1", "2", "3"))
 
> subset(data1, ID %in% c("1", "2", "3"))
    ID drug length treatment bdi.pre month bdi
1    1   No    >6m       TAU      29     2   2
2    2  Yes    >6m     BtheB      32     2  16
3    3  Yes    <6m       TAU      25     2  20
101  1   No    >6m       TAU      29     4   2
102  2  Yes    >6m     BtheB      32     4  24
103  3  Yes    <6m       TAU      25     4  NA
201  1   No    >6m       TAU      29     6  NA
202  2  Yes    >6m     BtheB      32     6  17
203  3  Yes    <6m       TAU      25     6  NA
301  1   No    >6m       TAU      29     8  NA
302  2  Yes    >6m     BtheB      32     8  20
303  3  Yes    <6m       TAU      25     8  NA

IDのrandom-slopeモデルでLMM (Linear Mixed Model)

> model01 <- lmer(bdi ~ treatment + month + bdi.pre + drug + length + (1|ID), data = data1)
> summary(model01)
Linear mixed model fit by REML ['lmerMod']
Formula: bdi ~ treatment + month + bdi.pre + drug + length + (1 | ID)
   Data: data1

REML criterion at convergence: 1859.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7135 -0.4974 -0.0827  0.4084  3.7246 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 51.41    7.170   
 Residual             25.52    5.052   
Number of obs: 280, groups:  ID, 97

Fixed effects:
               Estimate Std. Error t value
(Intercept)     4.46830    2.26253   1.975
treatmentBtheB -2.35589    1.70968  -1.378
month4         -1.42557    0.81954  -1.739
month6         -2.64780    0.89391  -2.962
month8         -4.37403    0.92965  -4.705
bdi.pre         0.63909    0.07968   8.021
drugYes        -2.79416    1.76710  -1.581
length>6m       0.23614    1.67690   0.141

Correlation of Fixed Effects:
            (Intr) trtmBB month4 month6 month8 bdi.pr drugYs
tretmntBthB -0.395                                          
month4      -0.141  0.021                                   
month6      -0.138  0.025  0.440                            
month8      -0.125  0.019  0.423  0.432                     
bdi.pre     -0.691  0.122  0.012  0.025  0.018              
drugYes     -0.076 -0.323 -0.016 -0.023 -0.023 -0.237       
length>6m   -0.247  0.001 -0.034 -0.039 -0.039 -0.241  0.158
>

broom()を用いてtidyな感じで固定効果を抽出

> library(broom)
> es.model01 <- tidy(model01, effects = "fixed", conf.int = TRUE)
> es.model01
            term   estimate  std.error  statistic    conf.low  conf.high
1    (Intercept)  4.4683042 2.26253164  1.9749135  0.03382372  8.9027848
2 treatmentBtheB -2.3558889 1.70967729 -1.3779729 -5.70679483  0.9950170
3         month4 -1.4255691 0.81954421 -1.7394658 -3.03184629  0.1807080
4         month6 -2.6477953 0.89390922 -2.9620405 -4.39982513 -0.8957654
5         month8 -4.3740325 0.92964553 -4.7050541 -6.19610425 -2.5519607
6        bdi.pre  0.6390882 0.07967549  8.0211392  0.48292710  0.7952493
7        drugYes -2.7941559 1.76710183 -1.5812082 -6.25761183  0.6693001
8      length>6m  0.2361410 1.67689698  0.1408202 -3.05051673  3.5227986
>

大本のmodel01は線形混合モデルの結果です。
broom::tidy()を用いて綺麗にまとめたのがex.model01。引数conf.intをTRUEとすることで信頼区間を表示できます。信頼区間の水準はconf.levelで。effects = "fixed"を指定すると、固定効果のみを抽出できます。ランダム効果が効いてない場合にeffectsをデフォルト指定するとここでエラーを吐きます。
なお、glmer(family = binomial)、すなわちロジスティック回帰で一般化線形混合モデルを用いる場合、tidy関数のexponentiateは無効化されている模様 (glmではまだ有効？)。この辺、もう一工夫が必要そう。

broom()でまとめたlmer/glmerの結果をforestplotで出力

# plotting modeled effects
plot.es <- forestplot(
  labeltext = es.model01$term, mean = es.model01$estimate, lower = es.model01$conf.low, upper = es.model01$conf.high, # 変数名、効果量、信頼区間 (下限)、信頼区間 (上限)
  is.summary = FALSE, # メタアナっぽく表を同時表示するならTRUEに
  
  fn.ci_norm = fpDrawCircleCI, # マーカーの形状を指定
  boxsize = 0.3, # 主にマーカーの大きさを指定
  col = fpColors(box = "royalblue", lines = "gray", zero = "darkred"), # グラフ要素の色指定。
  colgap = unit(6, "mm"), # 変数間の間隔
  
  # grid = structure(c(-1, 1, 2, 3), gp = gpar(lty = 2, col = "#CCCCFF")), # グリッドの設定。お好みで。
  lineheight = unit(1, "cm"), # default: auto
  zero = 0, lwd.zero = 2, # X軸の0の設定。対数軸ならzero = 1、そうでなければzero = 0 (default)
  
  xtics = seq(from = -6, to = 9, by = 3), # X軸の軸ラベルの設定
  xlab = "BDI scores",
  txt_gp = fpTxtGp(label = gpar(fontfamily = "", cex = 1),
                   ticks = gpar(fontfamily = "", cex = 1),
                   xlab  = gpar(fontfamily = "", cex = 1)), # 文字のサイズ設定
  
  mar = unit(rep(5, times = 4), "mm") # グラフの四隅からのマージンの設定
)

ここでプロットされるグラフはggplot2系ではないので、jpeg()/png()とdev.off()で挟んであげましょう。
出力されたグラフがこちら。

forestplot()のis.summary以降は描画の細かい設定で、効果量や信頼区間の限界値はlabeltextの行で指定済。推定結果をtibbleで返してくれるbroom::tidy()が活きてます。