Rによる一般化線型モデル（GLM）

20252 ワード

テキストリンク

はじめに

本記事では、Rによる一般化線型モデル解析を紹介する。線型回帰、ロジスティック回帰、ポアソン回帰を行う。入門的な記事で記されている内容に加え、係数ベクトルによる算出、対比検定、offset項を用いたポアソン回帰による率比推定を記載した。

【参考文献】

一般線型モデルと一般化線型モデル
Package
一般線型モデル（LM）
- Cervical Dystonia longitudinal dataset
- Variables
- データの読み込み
- 線型回帰分析
- 係数ベクトルによる算出
  - 例1: treat_c2群の16週目のtwstrs
  - 例2: treat_c2群の治療効果
- 対比検定
一般化線型モデル（GLM）
- Byar & Greene prostate cancer data
- Variables
- データの読み込み
- ロジスティック回帰
  - 対比検定
- ポアソン回帰（率比の推定）
まとめ
- （メモ）今後、追加する可能性のある解析

一般線型モデルと一般化線型モデル

一般線型モデルでは以下のようなモデルを考える。

E(Y_i) = \mu_i = \bm{x_i^T\beta} \\ Y_i \sim N(\mu_i, \sigma^2)

$Y_i$ は各対象者でのアウトカム、 $x_i$ は各対象者での説明変数とする。単回帰分析、重回帰分析、ANOVA、ANCOVAは全てこのモデルで書き表せられる。一般線型モデルをLinear Model（LM）と呼ぶ。ちなみに一般線型モデルをGeneral linear modelと記載することもあるが、次に述べる一般化線型モデル Generalized linear modelと区別しづらいため本記事内ではLMで統一する。

LMを正規分布以外の分布（指数型分布族）に拡張したモデルが一般化線型モデルであり、Generalized linear model（GLM）と呼ばれる。GLMでは以下のようなモデルを考える。

E[g(Y_i)] = \mu_i = \bm{x_i^T\beta} \\

LMとの違いとして $g()$ という変換を考えることができ、この関数 $g()$ を連結関数やリンク関数と呼ぶ。リンク関数が恒等関数、つまり $g(Y_i)=Y_i$ のときLMと等しくなる。GLMの解析では、 $Y_i$ (厳密には誤差)が従う分布とリンク関数を考えることになる。正規分布であれば恒等関数、2項分布であればlogit、ポアソン分布であればlogがリンク関数としてしばしば利用されるが、必ずしも分布とリンク関数は1対1対応するものではない。とはいえ、標準的な統計ソフトウェアであれば分布を指定すればデフォルトでこれらのリンク関数が指定される。

Rを利用することで簡便に推定できるため本記事では述べないが、LMとGLMにおける $\bm{\beta}$ の推定方法は「一般化線形モデル入門原著第2版」に詳しい。

ここからは実際にRを利用してLMとGLMを行う。

Package

利用するPackageは以下の通り。LMやGLMを行うだけであれば特にPackageを追加しなくても可能だが、データハンドリングや係数ベクトルを用いた算出を行うために以下のPackageを利用した。

#package
library(tidyverse)
library(dplyr)
library(multcomp)

一般線型モデル（LM）

LMでは「Rによる経時データ解析」で用いたデータセットと同一のものを利用する。データの記述はそちらを参考にしていただきたい。

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が高いほど、障害の程度が重いことを意味する。

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

Variables

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

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

データの読み込み

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

th1 <- cdystonia %>% 
  mutate(treat_c = case_when(
    treat == "Placebo" ~ "0",
    treat == "5000U" ~ "1",
    treat == "10000U" ~ "2" )) %>% 
  mutate(treat_c = as.factor(treat_c))

th2 <- th1 %>% filter(week==16)

線型回帰分析

16週目のtwstrsを結果変数として、治療と性別、年齢をモデルに含めた回帰分析を行った。治療と性別の交互作用項、治療と年齢の交互作用項もモデルに含めた。特に分布やリンク関数を指定しない場合は、正規分布と恒等関数が用いられる。

Coefficientsを見れば、各推定値がわかる。解釈について少し言及すると、ageの0.4238はtreat_c0群で年齢が1つあがるとスコアが0.4238を上がることを意味する。そのため年齢が10違えば、4.238だけスコアは変化する。またtreat_c1群での年齢の変化は単にageを見るだけでは不十分で、age + treat_c1:age の-0.098となり、年齢が上がるにつれてスコアが減少することになる。臨床的な解釈はわからないが、交互作用項をモデルに含めると各推定値の意味がわかりづらくなるので注意が必要である。各群のスコアや治療効果を算出する方法は次のセクションで述べる。

# LM
lm1 <- glm(twstrs ~ treat_c + sex + treat_c:sex + age + treat_c:age, data = th2)

summary(lm1)

## 
## Call:
## glm(formula = twstrs ~ treat_c + sex + treat_c:sex + age + treat_c:age, 
##     data = th2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -27.279   -7.069    1.396    7.256   27.964  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    21.3994     9.4547   2.263   0.0259 *
## treat_c1       26.8725    13.2492   2.028   0.0453 *
## treat_c2       26.3271    13.2773   1.983   0.0502 .
## sexM           -3.3703     4.0698  -0.828   0.4096  
## age             0.4238     0.1630   2.600   0.0108 *
## treat_c1:sexM   7.6661     5.6428   1.359   0.1775  
## treat_c2:sexM   0.3914     6.1604   0.064   0.9495  
## treat_c1:age   -0.5218     0.2278  -2.290   0.0242 *
## treat_c2:age   -0.3909     0.2320  -1.685   0.0952 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 132.3536)
## 
##     Null deviance: 14727  on 104  degrees of freedom
## Residual deviance: 12706  on  96  degrees of freedom
## AIC: 821.54
## 
## Number of Fisher Scoring iterations: 2

また信頼区間を推定したい場合にはconfint()を用いればよく、ここでは95%信頼区間を出力している。

confint(lm1)

## Waiting for profiling to be done...

##                     2.5 %      97.5 %
## (Intercept)     2.8686251 39.93020797
## treat_c1        0.9045071 52.84054792
## treat_c2        0.3040626 52.35011590
## sexM          -11.3469279  4.60628566
## age             0.1042815  0.74340968
## treat_c1:sexM  -3.3937045 18.72581962
## treat_c2:sexM -11.6828603 12.46561872
## treat_c1:age   -0.9682283 -0.07528326
## treat_c2:age   -0.8455945  0.06377826

係数ベクトルによる算出

一般的な入門向けの説明は以上であり、ここからは難易度が上がるように感じられるかもしれないが、交互作用項をモデルに含めた場合の各推定値の解釈はやや難しくなる。例えば、treat_c2の26.3は0歳の女性におけるtreat_c2 - treat_c0 の値であり意味を持たない（そもそも0歳まで解釈できるモデルではない可能性が高い）。

このように場合、回帰係数から興味のある推定量を算出する操作が必要になることがある。そのときには、各推定値をどのように足し合わせるのかは目的に応じて自ら係数ベクトルを指定する必要がある。

例えば、以下では各Treat群の16週目のtwstrsの値を求めることを考えている。各Treat群の16週目のtwstrsの値は共変量の分布、ここでは性別と年齢に依存するが男性と女性に同様の重みをつけ、年齢は今回の集団の平均値を与えた。

先ほどの結果の出力をみると、
intercept, treat_c1, treat_c2, sexM, age, treat_c1:sexM, treat_c2:sexM, treat_c1:age, treat_c2:age
の順になっており、この順に合わせてmatrixを指定している。例えば、treat_c2ではinterceptとtreat_c2に1を、sexMとtreat_c2:sexMに1/2を、ageとtreat_c2:ageにage_meanを、その他には0を指定している。multcomp::glhtを用いて計算を行うが、全パラメータに自分で値を与える必要がある点に注意が必要である。

age_mean <- mean(th2$age)

treat_c2 <- matrix(c(1, 0, 1, 1/2, age_mean, 0, 1/2, 0, age_mean), 1)
treat_c1 <- matrix(c(1, 1, 0, 1/2, age_mean, 1/2, 0, age_mean, 0), 1)
treat_c0 <- matrix(c(1, 0, 0, 1/2, age_mean, 0, 0, 0, 0), 1)

例1: treat_c2群の16週目のtwstrs

treat_c2群の16週目のtwstrsを算出するときに、係数ベクトルの指定が終われば、multcomp::glhtで計算でき、48.063であることがわかる。

effect_treat_c2 <- glht(lm1, linfct = treat_c2)
summary(effect_treat_c2)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = twstrs ~ treat_c + sex + treat_c:sex + age + treat_c:age, 
##     data = th2)
## 
## Linear Hypotheses:
##        Estimate Std. Error z value Pr(>|z|)    
## 1 == 0   48.063      2.308   20.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

例2: treat_c2群の治療効果

treat_c2の治療効果として、treat_c2 - treat_c0を考える場合も同様である。初めから差に該当するベクトルを指定してもよいが、ここでは先ほど指定した各群の差分として求めた。

treat_c2_c0 <- treat_c2 - treat_c0
effect_treat_c2_c0 <- glht(lm1, linfct = treat_c2_c0)
summary(effect_treat_c2_c0)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = twstrs ~ treat_c + sex + treat_c:sex + age + treat_c:age, 
##     data = th2)
## 
## Linear Hypotheses:
##        Estimate Std. Error z value Pr(>|z|)
## 1 == 0    4.852      3.071    1.58    0.114
## (Adjusted p values reported -- single-step method)

confint(effect_treat_c2_c0)

## 
##   Simultaneous Confidence Intervals
## 
## Fit: glm(formula = twstrs ~ treat_c + sex + treat_c:sex + age + treat_c:age, 
##     data = th2)
## 
## Quantile = 1.96
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##        Estimate lwr     upr    
## 1 == 0  4.8516  -1.1666 10.8697

対比検定

対比検定も同じようにできる。介入や曝露に用量反応性等の関連がある場合には、その関連を考慮したほうが検出力は向上する。ここでは等間隔の線型な関連が用量反応性があると仮定したときの対比検定を実施した。そのためtreat_c0, treat_c1, treat_c2に（-1, 0, 1）を当てている。この数をn倍したとしても、Estimateは変化するがz-value, p値は変化しない。

#対比検定
trt_contrast <- 1*treat_c2 + 0*treat_c1 - 1*treat_c0
trt_contrast2 <- glht(lm1, linfct = trt_contrast)
summary(trt_contrast2)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = twstrs ~ treat_c + sex + treat_c:sex + age + treat_c:age, 
##     data = th2)
## 
## Linear Hypotheses:
##        Estimate Std. Error z value Pr(>|z|)
## 1 == 0    4.852      3.071    1.58    0.114
## (Adjusted p values reported -- single-step method)

一般化線型モデル（GLM）

GLMでは「Rによる生存時間解析」で用いたデータセットと同一のものを利用する。データの記述はそちらを参考にしていただきたい。GLMの解析は基本的にはLMと同じように扱える。そのため被る内容は扱っていない。

Byar & Greene prostate cancer data

前立腺がん治療のランダム化比較試験（Byar and Green, 1980）のデータを利用した。前立腺がんのStage3、Stage4の患者を対象にエストロゲン治療群（0.2 mg, 1.0 mg, 5.0 mg）とPlacebo群の4群に割付が行われた。データはVanderbilt Biostatistics Datasetsの"Byar & Greene prostate cancer data"からダウンロードが可能である（2022-2-25時点）。

Variables

・データに含まれている変数のうち利用した変数

変数名	概要
stage	Stage(3 or 4)
rx	割付群（placebo, 0.2 mg estrogen, 1.0 mg estrogen or 5.0 mg estrogen）
dtime	追跡期間（月）
status	alive もしくは死因

・新たに作成した変数

変数名	概要
rx_c	rxを数値化
status2	生存、前立腺がんによる死亡、心疾患による死亡、その他に分類
status3-5	status2のダミー変数（status5が前立腺がんによる死亡）
stage_c	StageをFactorにした変数
dtime_y	追跡期間（年）
log_dtime_y	追跡期間（年）を対数変換した変数

データの読み込み

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

th1 <- prostate %>% 
  mutate(rx_c = case_when(
    rx == "placebo" ~ 0,
    rx == "0.2 mg estrogen" ~ 1,
    rx == "1.0 mg estrogen" ~ 2,
    rx == "5.0 mg estrogen" ~ 3,
  )) %>% 
  mutate(rx_c=as.factor(rx_c)) %>% 
  mutate(status2 = case_when(
    status == "alive" ~ 0,
    status == "dead - prostatic ca" ~ 1,
    status %in% c("dead - heart or vascular","dead - cerebrovascular") ~ 2,
    TRUE ~ 3
  )) %>% 
  mutate(status3 =
           ifelse(status2==0,0,
                  ifelse(status2==2,1,0))) %>% 
  mutate(status4 = 
           ifelse(status2==0,0,
                  ifelse(status2==3,1,0))) %>% 
  mutate(status5 = 
           ifelse(status2==1,1,0)) %>% 
  mutate(stage_c=as.factor(stage)) %>% 
  mutate(dtime_y=dtime/12) %>% 
  mutate(log_dtime_y=log(dtime/12))

ロジスティック回帰

ロジスティック回帰を行った。LMとは異なり、family=binomial(link = "logit")を指定している。また推定される値は対数オッズ比のため、exp()でオッズ比にしている。stage_c4の対数オッズ比1.19とオッズ比3.29は扱いがやや異なり、以下の式からもわかるように対数オッズ比の場合はその他のEstimateに足す形で処理する（左辺）が、オッズ比の場合はその他のEstimateに掛ける形で処理する（右辺）ことになる。

e^{{\beta_1}+{\beta_2}} = e^{\beta_1} \cdot e^{\beta_2}

2値のイベントに対してロジスティック回帰ではオッズ比が推定される。リスク比やリスク差を推定する回帰分析については、「Rによる計数データ（カテゴリカルデータ）解析」の回帰分析によるリスク比、オッズ比、リスク差をご参考いただきたい。

# Logistic
glm1 <- glm(status5 ~ rx_c + stage_c + rx_c:stage_c ,
            data = th1,
            family=binomial(link = "logit"))

summary_glm1 <- summary(glm1)
OR <- exp(glm1[["coefficients"]])
OR_95cl <- exp(confint(glm1))

## Waiting for profiling to be done...

cbind(summary_glm1[["coefficients"]], OR, OR_95cl)

##                  Estimate Std. Error    z value     Pr(>|z|)        OR
## (Intercept)    -1.4552872  0.2968084 -4.9031198 9.432643e-07 0.2333333
## rx_c1          -0.3852624  0.4516284 -0.8530517 3.936307e-01 0.6802721
## rx_c2          -1.3631110  0.5941473 -2.2942309 2.177724e-02 0.2558635
## rx_c3          -0.7576857  0.4965656 -1.5258522 1.270467e-01 0.4687500
## stage_c4        1.1895841  0.4060872  2.9293807 3.396382e-03 3.2857143
## rx_c1:stage_c4  1.1722625  0.6038720  1.9412432 5.222878e-02 3.2292906
## rx_c2:stage_c4  1.0691984  0.7130175  1.4995401 1.337336e-01 2.9130435
## rx_c3:stage_c4  0.4927606  0.6346658  0.7764095 4.375072e-01 1.6368286
##                     2.5 %     97.5 %
## (Intercept)    0.12529909  0.4048369
## rx_c1          0.27383783  1.6372406
## rx_c2          0.06957166  0.7579342
## rx_c3          0.16768236  1.2072610
## stage_c4       1.50011789  7.4288712
## rx_c1:stage_c4 0.99751076 10.7370592
## rx_c2:stage_c4 0.75981264 12.9834588
## rx_c3:stage_c4 0.47935964  5.8631471

対比検定

対比検定の考え方はLMとGLMでそれほど変わらないが、ここでは4群になっているのでその点について触れる。K1は等間隔の用量反応性を仮定しており順に（-3, -1, 1, 3）を当てて、K2は2次の反応性を仮定しており順に（1, -1, -1, 1）を当てている。今回は2次の関連性を想定する妥当な理由はないが、対比検定では1次以外にも2次あるいは3次などの関連性を指定できるため、その1例として行った。

# int c1 c2 c3 stage4 c1 c2 c3
rx_c3 <- matrix(c(1, 0, 0, 1, 1/2, 0, 0, 1/2), 1)
rx_c2 <- matrix(c(1, 0, 1, 0, 1/2, 0, 1/2, 0), 1)
rx_c1 <- matrix(c(1, 1, 0, 0, 1/2, 1/2, 0, 0), 1)
rx_c0 <- matrix(c(1, 0, 0, 0, 1/2, 0, 0, 0), 1)

K1 <- -3 * rx_c0 + -1 * rx_c1 + 1 * rx_c2 + 3 *rx_c3
K2 <- 1 * rx_c0 + -1 * rx_c1 -1 * rx_c2 + 1 *rx_c3

# 1次
contrast1 <- glht(glm1, linfct = K1)
summary(contrast1)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = status5 ~ rx_c + stage_c + rx_c:stage_c, family = binomial(link = "logit"), 
##     data = th1)
## 
## Linear Hypotheses:
##        Estimate Std. Error z value Pr(>|z|)  
## 1 == 0   -2.563      1.021  -2.511    0.012 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

# 2次
contrast2 <- glht(glm1, linfct = K2)
summary(contrast2)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = status5 ~ rx_c + stage_c + rx_c:stage_c, family = binomial(link = "logit"), 
##     data = th1)
## 
## Linear Hypotheses:
##        Estimate Std. Error z value Pr(>|z|)
## 1 == 0   0.1163     0.4863   0.239    0.811
## (Adjusted p values reported -- single-step method)

ポアソン回帰（率比の推定）

ポアソン回帰を行った。ポアソン回帰は疫学などで発生率（/人時間）の比を推定するときに用いられることがある。そのときにはOffset項に対数時間を指定する。以下のような式変形をすると、 $\bm{\beta}$ が対数率比として解釈できることがわかる。

E[log(Y_i)] = \bm{x_i^T\beta} + log(T_i)\\ E[log(Y_i)-log(T_i)] = \bm{x_i^T\beta} \\ E[log(\frac{Y_i}{T_i})] = \bm{x_i^T\beta} \\

対数率比が推定されるため、ロジスティック回帰と同様にexp()で率比にしている。

th2 <- th1 %>% filter(dtime_y!=0)

glm2<- glm(status5 ~ rx_c + stage_c + rx_c:stage_c ,
           data = th2,
           family=poisson(link = "log"),
           offset = log_dtime_y)

summary_glm2 <- summary(glm2)
RR <- exp(glm2[["coefficients"]])
RR_95cl <- exp(confint(glm2))

## Waiting for profiling to be done...

cbind(summary_glm1[["coefficients"]], RR, RR_95cl)

##                  Estimate Std. Error    z value     Pr(>|z|)         RR
## (Intercept)    -1.4552872  0.2968084 -4.9031198 9.432643e-07 0.05414115
## rx_c1          -0.3852624  0.4516284 -0.8530517 3.936307e-01 0.83733607
## rx_c2          -1.3631110  0.5941473 -2.2942309 2.177724e-02 0.29860944
## rx_c3          -0.7576857  0.4965656 -1.5258522 1.270467e-01 0.53236875
## stage_c4        1.1895841  0.4060872  2.9293807 3.396382e-03 3.68870167
## rx_c1:stage_c4  1.1722625  0.6038720  1.9412432 5.222878e-02 1.30702094
## rx_c2:stage_c4  1.0691984  0.7130175  1.4995401 1.337336e-01 1.94665411
## rx_c3:stage_c4  0.4927606  0.6346658  0.7764095 4.375072e-01 1.20584587
##                     2.5 %     97.5 %
## (Intercept)    0.03049911 0.08764829
## rx_c1          0.36097703 1.87140007
## rx_c2          0.08464247 0.83242477
## rx_c3          0.18851837 1.32699945
## stage_c4       1.92077836 7.34503934
## rx_c1:stage_c4 0.49723191 3.54146221
## rx_c2:stage_c4 0.58511543 7.76813429
## rx_c3:stage_c4 0.40200019 3.95166618

まとめ

本記事では、Rによる一般化線型モデル解析として線型回帰、ロジスティック回帰、ポアソン回帰を行った。これらの解析は簡単にRで実装できるが、推定された値を解釈するためにはある程度の一般化線型モデルへの理解が求められる。本記事が多少でもお役に立てば幸いである。

また記事内に誤りが見つかった場合には、コメント等でご連絡下さい。

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

解析ではないが、multcomp::glhtのSE、95%信頼区間の推定方法について調べて追記（付録）

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

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

Author And Source

この問題について(Rによる一般化線型モデル（GLM）), 我々は、より多くの情報をここで見つけました https://zenn.dev/masahiro_kondo/articles/ea502d7594e48b

著者帰属：元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。

Collection and Share based on the CC protocol

Dockerでprotocを行って開発を楽にする

Tekton Pipelineを触ってみた on GKE Autopilot

Rによる一般化線型モデル（GLM）

はじめに

目次

一般線型モデルと一般化線型モデル

Package

一般線型モデル（LM）

Cervical Dystonia longitudinal dataset

Variables

データの読み込み

線型回帰分析

係数ベクトルによる算出

例1: treat_c2群の16週目のtwstrs

例2: treat_c2群の治療効果

対比検定

一般化線型モデル（GLM）

Byar & Greene prostate cancer data

Variables

データの読み込み

ロジスティック回帰

対比検定

ポアソン回帰（率比の推定）

まとめ

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

Author And Source