Rによる一般化線型モデル(GLM)
はじめに
本記事では、Rによる一般化線型モデル解析を紹介する。線型回帰、ロジスティック回帰、ポアソン回帰を行う。入門的な記事で記されている内容に加え、係数ベクトルによる算出、対比検定、offset項を用いたポアソン回帰による率比推定を記載した。
【参考文献】
目次
- 一般線型モデルと一般化線型モデル
- Package
- 一般線型モデル(LM)
- Cervical Dystonia longitudinal dataset
- Variables
- データの読み込み
- 線型回帰分析
- 係数ベクトルによる算出
- 例1: treat_c2群の16週目のtwstrs
- 例2: treat_c2群の治療効果
- 対比検定
- 一般化線型モデル(GLM)
- Byar & Greene prostate cancer data
- Variables
- データの読み込み
- ロジスティック回帰
- 対比検定
- ポアソン回帰(率比の推定)
- まとめ
- (メモ)今後、追加する可能性のある解析
一般線型モデルと一般化線型モデル
一般線型モデルでは以下のようなモデルを考える。
LMを正規分布以外の分布(指数型分布族)に拡張したモデルが一般化線型モデルであり、Generalized linear model(GLM)と呼ばれる。GLMでは以下のようなモデルを考える。
LMとの違いとして
Rを利用することで簡便に推定できるため本記事では述べないが、LMとGLMにおける
ここからは実際に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に掛ける形で処理する(右辺)ことになる。
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項に対数時間を指定する。以下のような式変形をすると、
対数率比が推定されるため、ロジスティック回帰と同様に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
Discussion
初歩的な質問で恐縮なのですが質問させていただきたいです.
質問は,Cervical Dystonia longitudinal dataset の読み込み方法がわからないというものです.
方法としてはリンク先から,cdystonia.savをダウンロードし,それをRで使えるように読み込むという方法でよろしいでしょうか.
私の環境だけかもしれませんが,cdystonia.savのダウンロード後,作成したディレクトリ下でforeignパッケージやhavenパッケージを用いてsavファイルを読み込もうとしても,エラー (file has unsupported features.) が出てしまい読み込めません.ディレクトリの指定先が間違っているというようなエラーではないかと思っています.savファイルはどのように読み込めばよいのでしょうか.
ちなみに,Rはv4.2.0です.
私は統計専門ではなく,農学系の修士学生ですが,GLMについて勉強したいと思い本ページを見させていただきました.
返信いただけますと幸いです.よろしくお願いいたします.
ぞんじん様
コメントありがとうございます。こちらsavではなく、dtaから読み込むことは可能でしょうか?僕も試してみたのですが、なぜかsavだとうまくRで読み込めず、dtaであれば読み込めました。
サンプルコード
library(foreign)
th1 <- read.dta("path/cdystonia.dta")
返信ありがとうございます.
dtaから読み込むことができました! m(_ _)m
勉強になりました.ありがとうございます.