Zenn
🐝

Rによる経時データ解析

2022/03/18に公開
R
#
統計
#
生物統計
#
経時データ解析
tech

はじめに

本記事では、医療分野における経時データ解析でしばしば用いられる変量効果モデル、MMRM、一般化推定方程式のRでの実装を紹介する。解析できそうなpackageを調べてはみたものの、分散の推定方法や自由度、積分の計算方法はあまり調べていない。また生物統計の初学者よりは一定の知識がある者を対象とし数式やRの使用方法に関する記述を省いたためわかりづらい箇所があるかもしれないが、ご容赦いただければ幸いである。

本記事では、説明のために一般線型モデルと一般化線型モデルで区別している。しかし尤度に基づく手法とは異なり、GEEは分布に仮定を置かない手法でありこのような区別は適していない可能性がある点は始めに述べておく。

経時データ解析でよい教科書やチュートリアル論文を知らないので、ご存じの方がいれば教えて頂きたい。変量効果についての「医学統計のための線型混合モデル―SASによるアプローチ」は絶版のため、英語に苦手意識がなければ原著を読まれるといいかもしれない。

【参考文献】

目次

  • Package
  • 線形モデル（Linear model）
    • Cervical Dystonia longitudinal dataset
    • Variables
    • データの前処理
    • データの要約
    • 線型混合効果モデル(LMM)
    • Mixed-effects models for repeated measures（MMRM）
    • 一般化推定方程式(GEE)
  • 一般化線型モデル（Generalized linear model）
    • Diggle, Liang, and Zeger (1994)
    • データの要約
    • 一般化線型混合効果モデル（GLMM）
    • 一般化線型混合効果モデルでのMMRM
    • 一般化推定方程式（GEE）
  • まとめ
    • （メモ）今後、追加する可能性のある解析

Package

利用するPackageは以下の通り。

#package
library(tidyverse)
library(reshape2)
library(gtsummary)
library(ggplot2)
library(lme4)
library(lmerTest)
library(geepack)
library(nlme)
library(glmmTMB)

線形モデル（Linear model）

Cervical Dystonia longitudinal dataset

2002年Springer出版のCharles S. Davisによる"Statistical Methods for the Analysis of Repeated Measurements"の表6.9のデータセットを利用した。多施設ランダム化比較試験で、頸部ジストニア患者に対するB型ボツリヌス毒素の効果が検討された。割付はPlacebo群、5000 units、10000 unitsの3群であり、ベースライン（0週目）および治療開始後2、4、8、12、16週目にToronto Western Spasmodic Torticollis Rating Scale (TWSTRS)が測定された。TWSTRSが高いほど、障害の程度が重いことを意味する。

本記事の解析方針として、Placebo群と投薬群の2群間の比較を行い、解析には　4, 8, 12, 16週目のtwstrsのスコアを用いて、治療効果として16週時点でのスコアの群間差を検討する。

データはVanderbilt Biostatistics Datasetsの" Cervical Dystonia longitudinal dataset"からダウンロードが可能である（2022-3-3時点）

Variables

・データに含まれている変数

変数名 概要
week 経過した週
site 施設識別コード
id 患者コード
treat 割付群（placebo, 5000 units or 10000 units）
age 年齢
sex 性別
twstrs twstrsのスコア

データの前処理

cdystonia <- read.table("path/cdystonia.csv",
                       sep=",",
                       header=TRUE,
                       na.strings=c(''),
                       fill=TRUE)

th1 <- cdystonia %>% 
  filter(week!=2) %>% 
  mutate(site_id = id) %>% 
  mutate(tmp_id = paste(site, id, sep="_")) %>% 
  mutate(treat_c = case_when(
    treat == "Placebo" ~ 0,
    treat == "5000U" ~ 1,
    treat == "10000U" ~ 1 )) %>% 
  mutate(treat_c = as.factor(treat_c)) %>% 
  mutate(week_c = week/4) %>% 
  mutate(site_c = as.factor(site))

th2 <- th1 %>% select(tmp_id, week_c, week, twstrs)

base <- th1 %>% names() %>% dput()

## c("week", "site", "id", "treat", "age", "sex", "twstrs", "site_id", 
## "tmp_id", "treat_c", "week_c", "site_c")

base <- base[c(-1,-3,-7,-11)]

ID <- c(1:nrow(filter(th1, week==0)))

visit <- data.frame()
for (i in ID) {
  visit <- rbind(visit, c(0:4))
} 

th3 <- cbind(ID, select(filter(th1, week==0), base), visit)

th4 <- melt(th3, c("ID",base),value.name="week_c")

th5 <- th4[order(th4$ID, th4$week_c),]
rownames(th5) <- 1:nrow(th5)

th6 <- left_join(th5, th2, by = c("tmp_id", "week_c"),all=FALSE)

dat <- th6 %>% filter(week_c != 0) %>% 
  mutate(week_c = as.factor(week_c))

データの要約

背景因子

性別、年齢、施設ごとの人数を表1に要約した。年齢は中央値と四分位範囲、カテゴリカル変数は人数と割合を算出した。

# データの要約(背景因子)
th10 <- th6 %>% 
  filter(week == 0) %>% 
  dplyr::select(treat,sex,age,site_c)

table1 <- tbl_summary(
  th10,
  by = treat,
  missing = "ifany") %>% 
  add_n() %>% 
  modify_header(label = "**Variable**")

twstrs

twstrsの推移を示した。平均値を折れ線グラフで表すと、治療群では4週目ではPlacebo群に対して改善が見られるが、徐々に差は縮まり16週時点ではPlacebo群のほうが治療群よりもtwstrsの値が小さくなった。

# データの要約(アウトカムの推移)
th11 <- summarize(group_by(th6,treat,week), 
          n = n(), 
          mean = mean(twstrs), 
          sd = sd(twstrs))

## `summarise()` has grouped output by 'treat'. You can override using the
## `.groups` argument.

g <- ggplot(th11, aes(x=week, y=mean, color=treat)) + 
  geom_line()


plot(g)

線型混合効果モデル（LMM）

lmerを利用して、線型混合効果モデルによる解析を行った。3つのモデルで解析を行った。16週目時点の群間差を算出するために係数から計算することを避けるために切片項を推定しないモデルで解析を行った。

LMM1では患者（ID）のみを変量効果に含めた解析を行い、LMM2では患者（ID）と施設（site_c）を変量効果に加えた解析を行った。また、LMM3ではIDと、連続量であるweekとの交互作用を変量効果とした解析を行った。

モデルの指定の仕方だが、(1 + week | ID)のように "|"のまえに変量効果に対するデザイン行列を示すことになる。よく書かれる式で言えば、固定効果とそのデザイン行列を\boldsymbol{\beta, X}として変量効果とそのデザイン行列を\boldsymbol{Z, b}としたときの\boldsymbol{Z}にあたる。交互作用項のみを指定したい場合には(0 + week | ID)とすればよく、単に(week | ID)とすると切片項がデフォルトで含まれてしまうことがある。

16週目の群間差の結はFixed effectsのtreat_c1:week_c4を見ればよい。

\boldsymbol {Y = X\beta + Zb + \epsilon}

LMM1

# LMM
LMM1 <- lmer(twstrs ~ 0 + treat_c:week_c + week_c + (1 | ID), data = dat)
summary(LMM1, ddf = "Kenward-Roger")

## Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
## lmerModLmerTest]
## Formula: twstrs ~ 0 + treat_c:week_c + week_c + (1 | ID)
##    Data: dat
## 
## REML criterion at convergence: 2946.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3647 -0.5171 -0.0571  0.5202  2.9196 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 126.92   11.266  
##  Residual              35.12    5.926  
## Number of obs: 419, groups:  ID, 108
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## week_c1           39.3429     2.1517 148.2164  18.285   <2e-16 ***
## week_c2           41.4000     2.1517 148.2164  19.241   <2e-16 ***
## week_c3           42.3199     2.1616 150.6078  19.578   <2e-16 ***
## week_c4           43.4964     2.1616 150.6078  20.123   <2e-16 ***
## treat_c1:week_c1  -3.3384     2.6206 148.9052  -1.274    0.205    
## treat_c1:week_c2  -2.2814     2.6247 149.7124  -0.869    0.386    
## treat_c1:week_c3   0.8969     2.6310 150.9595   0.341    0.734    
## treat_c1:week_c4   3.5940     2.6293 150.6000   1.367    0.174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             wek_c1 wek_c2 wek_c3 wek_c4 t_1:_1 t_1:_2 t_1:_3
## week_c2      0.783                                          
## week_c3      0.780  0.780                                   
## week_c4      0.780  0.780  0.779                            
## trt_c1:wk_1 -0.821 -0.643 -0.640 -0.640                     
## trt_c1:wk_2 -0.642 -0.820 -0.639 -0.639  0.780              
## trt_c1:wk_3 -0.641 -0.641 -0.822 -0.640  0.778  0.778       
## trt_c1:wk_4 -0.641 -0.641 -0.640 -0.822  0.779  0.778  0.778

LMM2

LMM2 <- lmer(twstrs ~ 0 + treat_c:week_c + week_c + (1 | ID) + (1 | site_c), data = dat)
summary(LMM2, ddf = "Kenward-Roger")

## Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
## lmerModLmerTest]
## Formula: twstrs ~ 0 + treat_c:week_c + week_c + (1 | ID) + (1 | site_c)
##    Data: dat
## 
## REML criterion at convergence: 2936.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4075 -0.5296 -0.0427  0.5110  2.8693 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 104.08   10.202  
##  site_c   (Intercept)  22.54    4.748  
##  Residual              35.15    5.929  
## Number of obs: 419, groups:  ID, 108; site_c, 9
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)    
## week_c1            39.388      2.560  25.052  15.384 2.85e-14 ***
## week_c2            41.445      2.560  25.052  16.187 8.90e-15 ***
## week_c3            42.342      2.569  25.374  16.484 4.54e-15 ***
## week_c4            43.518      2.569  25.374  16.942 2.40e-15 ***
## treat_c1:week_c1   -3.108      2.434 146.500  -1.277    0.204    
## treat_c1:week_c2   -2.024      2.438 147.390  -0.830    0.408    
## treat_c1:week_c3    1.162      2.445 148.747   0.475    0.635    
## treat_c1:week_c4    3.856      2.443 148.418   1.578    0.117    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             wek_c1 wek_c2 wek_c3 wek_c4 t_1:_1 t_1:_2 t_1:_3
## week_c2      0.846                                          
## week_c3      0.844  0.844                                   
## week_c4      0.844  0.844  0.843                            
## trt_c1:wk_1 -0.639 -0.478 -0.476 -0.476                     
## trt_c1:wk_2 -0.477 -0.638 -0.475 -0.475  0.745              
## trt_c1:wk_3 -0.476 -0.476 -0.641 -0.476  0.742  0.742       
## trt_c1:wk_4 -0.476 -0.476 -0.476 -0.641  0.743  0.743  0.743

LMM3

LMM3 <- lmer(twstrs ~ 0 + treat_c:week_c + week_c + (1 + week | ID), data = dat)
summary(LMM3, ddf = "Kenward-Roger")

## Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
## lmerModLmerTest]
## Formula: twstrs ~ 0 + treat_c:week_c + week_c + (1 + week | ID)
##    Data: dat
## 
## REML criterion at convergence: 2908.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7250 -0.4252 -0.0126  0.4439  2.5680 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID       (Intercept) 208.8740 14.4525       
##           week          0.4709  0.6862  -0.64
##  Residual              22.5179  4.7453       
## Number of obs: 419, groups:  ID, 108
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## week_c1           39.3429     2.3196 113.1383  16.961   <2e-16 ***
## week_c2           41.4000     2.1403 128.8737  19.344   <2e-16 ***
## week_c3           42.2575     2.0610 132.5187  20.504   <2e-16 ***
## week_c4           43.3685     2.0801 115.8341  20.850   <2e-16 ***
## treat_c1:week_c1  -3.4093     2.8248 113.5869  -1.207    0.230    
## treat_c1:week_c2  -2.2768     2.6085 129.8079  -0.873    0.384    
## treat_c1:week_c3   0.9875     2.5080 132.7076   0.394    0.694    
## treat_c1:week_c4   3.7578     2.5298 115.7659   1.485    0.140    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             wek_c1 wek_c2 wek_c3 wek_c4 t_1:_1 t_1:_2 t_1:_3
## week_c2      0.852                                          
## week_c3      0.779  0.826                                   
## week_c4      0.666  0.753  0.820                            
## trt_c1:wk_1 -0.821 -0.700 -0.639 -0.547                     
## trt_c1:wk_2 -0.699 -0.821 -0.678 -0.618  0.849              
## trt_c1:wk_3 -0.640 -0.679 -0.822 -0.674  0.777  0.825       
## trt_c1:wk_4 -0.548 -0.619 -0.674 -0.822  0.666  0.753  0.819

Mixed-effects models for repeated measures（MMRM）

MMRMは変量効果を指定するのではなく周辺化分散構造を指定する。UnstructuredやAR(1),CSなどがよく使われ、以下の例では無構造（unstructured）を指定している。

# MMRM（UN）
mmrm <- gls(twstrs ~ 0 + treat_c:week_c + week_c,
            na.action=na.exclude,
            dat, 
            correlation = corSymm(form = ~ 1 | ID))

summary(mmrm)

## Generalized least squares fit by REML
##   Model: twstrs ~ 0 + treat_c:week_c + week_c 
##   Data: dat 
##        AIC      BIC    logLik
##   2925.232 2985.511 -1447.616
## 
## Correlation Structure: General
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##  Correlation: 
##   1     2     3    
## 2 0.865            
## 3 0.760 0.822      
## 4 0.680 0.702 0.855
## 
## Coefficients:
##                     Value Std.Error   t-value p-value
## week_c1          39.34286  2.148729 18.309830  0.0000
## week_c2          41.40000  2.148729 19.267207  0.0000
## week_c3          42.27246  2.158638 19.582932  0.0000
## week_c4          43.39562  2.164042 20.053041  0.0000
## treat_c1:week_c1 -3.38877  2.615565 -1.295618  0.1958
## treat_c1:week_c2 -2.39302  2.618597 -0.913854  0.3613
## treat_c1:week_c3  0.96862  2.625798  0.368885  0.7124
## treat_c1:week_c4  3.72225  2.629249  1.415709  0.1576
## 
##  Correlation: 
##                  wek_c1 wek_c2 wek_c3 wek_c4 t_1:_1 t_1:_2 t_1:_3
## week_c2           0.865                                          
## week_c3           0.756  0.819                                   
## week_c4           0.676  0.697  0.852                            
## treat_c1:week_c1 -0.822 -0.710 -0.621 -0.555                     
## treat_c1:week_c2 -0.710 -0.821 -0.672 -0.572  0.863              
## treat_c1:week_c3 -0.622 -0.673 -0.822 -0.701  0.758  0.820       
## treat_c1:week_c4 -0.556 -0.573 -0.701 -0.823  0.679  0.699  0.851
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.85339211 -0.72695440  0.06939305  0.73356757  2.75690514 
## 
## Residual standard error: 12.71205 
## Degrees of freedom: 419 total; 411 residual

一般化推定方程式(GEE)

線型混合効果モデルやMMRMは尤度ベースの解析だが、推定方程式を利用するGEEも経時データ解析では利用されることがある。GEEではMMRM同様に分散の構造を指定するが一般にロバスト分散（サンドイッチ分散）を利用する。ロバスト分散であれば、分散構造を誤って指定していても時点数に比べてサンプルサイズが十分に大きければ問題にならない。

本記事では線型混合効果モデルと一般化推定方程式の使い分けについてはあまり言及しないが、

  • 欠測の有無
  • 欠測の機序（MCAR, MAR, MNAR）や構造（monotone, non-motonone）
  • 周辺化の有無による解釈の違い

などを考慮するべきだと私は考えている。

下の例では分散構造に独立とAR(1)を指定している。

独立

# GEE
GEE1 <- geese(twstrs ~ 0 + treat_c:week_c + week_c , 
            id = ID,
            waves = week_c,
            data = dat, 
            corstr = "independence", 
            family=gaussian)

summary(GEE1)

## 
## Call:
## geese(formula = twstrs ~ 0 + treat_c:week_c + week_c, id = ID, 
##     waves = week_c, data = dat, family = gaussian, corstr = "independence")
## 
## Mean Model:
##  Mean Link:                 identity 
##  Variance to Mean Relation: gaussian 
## 
##  Coefficients:
##                   estimate   san.se        wald         p
## week_c1          39.342857 1.970369 398.6910355 0.0000000
## week_c2          41.400000 2.253777 337.4262343 0.0000000
## week_c3          41.735294 2.099311 395.2329968 0.0000000
## week_c4          42.911765 2.286324 352.2710985 0.0000000
## treat_c1:week_c1 -3.399195 2.552241   1.7738178 0.1829103
## treat_c1:week_c2 -2.400000 2.779117   0.7457766 0.3878167
## treat_c1:week_c3  1.750420 2.539738   0.4750144 0.4906893
## treat_c1:week_c4  4.017813 2.622268   2.3476069 0.1254755
## 
## Scale Model:
##  Scale Link:                identity 
## 
##  Estimated Scale Parameters:
##             estimate   san.se     wald p
## (Intercept) 158.2185 16.53285 91.58384 0
## 
## Correlation Model:
##  Correlation Structure:     independence 
## 
## Returned Error Value:    0 
## Number of clusters:   108   Maximum cluster size: 4

AR(1)

GEE2 <- geese(twstrs ~ 0 + treat_c:week_c + week_c , 
              id = ID,
              waves = week_c,
              data = dat, 
              corstr = "ar1", 
              family=gaussian)
summary(GEE2)

## 
## Call:
## geese(formula = twstrs ~ 0 + treat_c:week_c + week_c, id = ID, 
##     waves = week_c, data = dat, family = gaussian, corstr = "ar1")
## 
## Mean Model:
##  Mean Link:                 identity 
##  Variance to Mean Relation: gaussian 
## 
##  Coefficients:
##                    estimate   san.se        wald         p
## week_c1          39.3428571 1.970369 398.6910355 0.0000000
## week_c2          41.4000000 2.253777 337.4262343 0.0000000
## week_c3          42.2547413 2.119839 397.3242111 0.0000000
## week_c4          43.3571075 2.298908 355.6948015 0.0000000
## treat_c1:week_c1 -3.3891082 2.533731   1.7891621 0.1810283
## treat_c1:week_c2 -2.3176727 2.749425   0.7105922 0.3992473
## treat_c1:week_c3  0.9578033 2.535352   0.1427171 0.7055947
## treat_c1:week_c4  3.7624867 2.638691   2.0331660 0.1538997
## 
## Scale Model:
##  Scale Link:                identity 
## 
##  Estimated Scale Parameters:
##             estimate   san.se     wald p
## (Intercept) 158.2762 16.54296 91.53868 0
## 
## Correlation Model:
##  Correlation Structure:     ar1 
##  Correlation Link:          identity 
## 
##  Estimated Correlation Parameters:
##        estimate     san.se     wald p
## alpha 0.8573399 0.02786535 946.6234 0
## 
## Returned Error Value:    0 
## Number of clusters:   108   Maximum cluster size: 4

一般化線型モデル（Generalized linear model）

正規分布以外の分布を考えるGLMでの経時データ解析を行う。てんかん発作の回数をアウトカムにした"seizure"というデータセットを用いてポアソン分布のあてはめを行う。

Diggle, Liang, and Zeger (1994)

geepackに含まれている"seizure"というデータセットを利用する。てんかん発作の回数を経時的に記録したデータである。

## Diggle, Liang, and Zeger (1994) pp166-168, compare Table 8.10
seiz.l <- reshape(seizure,
                  varying=list(c("base","y1", "y2", "y3", "y4")),
                  v.names="y", 
                  times=0:4, 
                  direction="long")

th20 <- seiz.l[order(seiz.l$id, seiz.l$time),]
seiz.l <- seiz.l[order(seiz.l$id, seiz.l$time),] %>% 
  filter(time != 0)

データの要約

てんかん発作の回数についてベースラインと各時点での平均値から折れ線グラフを作成した。

th21 <- summarize(group_by(th20,trt,time), 
                  n = n(), 
                  mean = mean(y), 
                  sd = sd(y))


g2 <- ggplot(th21, aes(x=time, y=mean, color=factor(trt))) + 
  geom_line()

plot(g2)

一般化線型混合効果モデル（GLMM）

GLMMはglmerで行える。変量効果の指定の仕方は線型混合効果モデルの時と同様である。ここではポアソン分布をあてはめており、Link関数はLogである。出力される結果をx倍の単位には直していないので、解釈しやすくするためにはFixed effectsの推定値をexp()で直す必要がある。

## GLMM
GLMM <- glmer(y ~ 0 + trt:factor(time) + factor(time) + (1 | id), 
      data = seiz.l,
      family=poisson)

summary(GLMM)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: y ~ 0 + trt:factor(time) + factor(time) + (1 | id)
##    Data: seiz.l
## 
##      AIC      BIC   logLik deviance df.resid 
##   1407.9   1439.0   -694.9   1389.9      227 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3890 -0.7862 -0.1086  0.6203  6.6750 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 0.8713   0.9334  
## Number of obs: 236, groups:  id, 59
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## factor(time)1       1.8574     0.1891   9.822   <2e-16 ***
## factor(time)2       1.7358     0.1904   9.117   <2e-16 ***
## factor(time)3       1.7944     0.1898   9.456   <2e-16 ***
## factor(time)4       1.7007     0.1908   8.914   <2e-16 ***
## trt:factor(time)1  -0.2957     0.2628  -1.125    0.260    
## trt:factor(time)2  -0.1930     0.2638  -0.732    0.464    
## trt:factor(time)3  -0.2867     0.2636  -1.088    0.277    
## trt:factor(time)4  -0.3802     0.2659  -1.430    0.153    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             fct()1 fct()2 fct()3 fct()4 tr:()1 tr:()2 tr:()3
## factor(tm)2  0.889                                          
## factor(tm)3  0.892  0.886                                   
## factor(tm)4  0.887  0.881  0.884                            
## trt:fctr()1 -0.718 -0.638 -0.640 -0.637                     
## trt:fctr()2 -0.640 -0.720 -0.638 -0.634  0.888              
## trt:fctr()3 -0.640 -0.636 -0.718 -0.635  0.889  0.885       
## trt:fctr()4 -0.635 -0.631 -0.633 -0.716  0.881  0.878  0.878

一般化線型混合効果モデルでのMMRM

glmmTMBで、正規分布以外でもMMRMが行える。以下では分散構造は無構造でポアソン分布を当てはめて解析している。

## GLMM MMRM
# https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html
MMRM2 <- glmmTMB(y ~ 0 + trt:factor(time) + factor(time) + us(0 + factor(time)| id),
              data = seiz.l,
              family=poisson)

summary(MMRM2)

##  Family: poisson  ( log )
## Formula:          
## y ~ 0 + trt:factor(time) + factor(time) + us(0 + factor(time) |      id)
## Data: seiz.l
## 
##      AIC      BIC   logLik deviance df.resid 
##   1332.6   1395.0   -648.3   1296.6      218 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name          Variance Std.Dev. Corr           
##  id     factor(time)1 1.0242   1.0120                  
##         factor(time)2 0.7859   0.8865   0.87           
##         factor(time)3 1.2642   1.1244   0.81 0.89      
##         factor(time)4 0.7630   0.8735   0.91 0.95 0.92 
## Number of obs: 236, groups:  id, 59
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## factor(time)1       1.8184     0.2087   8.713  < 2e-16 ***
## factor(time)2       1.7689     0.1879   9.416  < 2e-16 ***
## factor(time)3       1.5586     0.2334   6.679  2.4e-11 ***
## factor(time)4       1.7612     0.1836   9.591  < 2e-16 ***
## factor(time)1:trt  -0.4177     0.2908  -1.436    0.151    
## factor(time)2:trt  -0.1264     0.2563  -0.493    0.622    
## factor(time)3:trt  -0.2336     0.3216  -0.726    0.468    
## factor(time)4:trt  -0.3815     0.2551  -1.496    0.135    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

一般化推定方程式（GEE）

一般化推定方程式であれば、正規分布以外でも正規分布のときと同様にプログラムすればよい。

## GEE
GEE3 <- geese(y ~ 0 + trt:factor(time) + factor(time), 
              id = id,data=seiz.l, corstr="exch", 
              waves = time,
              family=poisson)
summary(GEE3)

## 
## Call:
## geese(formula = y ~ 0 + trt:factor(time) + factor(time), id = id, 
##     waves = time, data = seiz.l, family = poisson, corstr = "exch")
## 
## Mean Model:
##  Mean Link:                 log 
##  Variance to Mean Relation: poisson 
## 
##  Coefficients:
##                      estimate    san.se         wald            p
## factor(time)1      2.23613999 0.2010414 1.237163e+02 0.000000e+00
## factor(time)2      2.11453286 0.1828579 1.337215e+02 0.000000e+00
## factor(time)3      2.17312703 0.3099234 4.916561e+01 2.352341e-12
## factor(time)4      2.07944154 0.1765511 1.387244e+02 0.000000e+00
## trt:factor(time)1 -0.08663089 0.4260137 4.135219e-02 8.388596e-01
## trt:factor(time)2  0.01600034 0.3088349 2.684144e-03 9.586811e-01
## trt:factor(time)3 -0.07768514 0.4327254 3.222933e-02 8.575253e-01
## trt:factor(time)4 -0.17109449 0.3437487 2.477365e-01 6.186734e-01
## 
## Scale Model:
##  Scale Link:                identity 
## 
##  Estimated Scale Parameters:
##             estimate  san.se     wald          p
## (Intercept) 18.26834 8.81939 4.290628 0.03832301
## 
## Correlation Model:
##  Correlation Structure:     exch 
##  Correlation Link:          identity 
## 
##  Estimated Correlation Parameters:
##        estimate     san.se     wald p
## alpha 0.8054581 0.09123654 77.93781 0
## 
## Returned Error Value:    0 
## Number of clusters:   59   Maximum cluster size: 4

まとめ

経時データについて、混合効果モデル、MMRM、一般化推定方程式を線型モデルと一般化線型モデルの場合での解析を紹介した。筆者の意見としては、変量効果や周辺化分散、推定方程式はどれもやや難解な点があり、経時データやクラスターデータのいかなる場面でも適用できるような解析は存在しないように思う。常に目的やデータの性質を考慮する必要があるため、これらの解析を行う際には統計の専門家に相談することを特に推奨したい分野である。もし本記事内で誤りを見つけられた方は、ご連絡頂ければ幸いである。

（メモ）今後、追加する可能性のある解析

  • 各手法の分散推定、計算アルゴリズム、検定時の自由度の算出方法

サポートして頂けるとモチベーションに繋がりますのでぜひ宜しくお願いします。データ解析や臨床研究でのご相談があれば、お気軽にTwiiterもしくはメールにてご連絡下さい。

作成者：Masahiro Kondo
作成日：2022/3/18
連絡先：m.kondo1042(at)gmail.com

Discussion