Relative Risk Regression
Relative Risk Regression
2022-06-29
Relative Risk Regression
결과값이 삶/죽음, 예/아니오 와 같이 이항
반응(binary response)인 경우 가장 널리 사용되는 분석 방법은 logistic
regression이다.
\[\textrm{logit}(\mathbb{E}[y]) = \beta_0
+ \beta_1 x_1 + \beta_2 x_2 + \cdots\]
이름이 로지스틱인 이유는 위의 식을 다음과 같이 변형할 수 있기
때문이다.
\[\mathbb{E}[y] = \textrm{logistic}(
\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots)\]
로지스틱 회귀의 링크 함수는 logit()이며, Relative Risk
Regression의 링크 함수는 다음과 같이 log()이다.
\[\log (\mathbb{E}[y]) = \beta_0 + \beta_1
x_1 + \beta_2 x_2 + \cdots\]
이 회귀식이 Relative Risk Regression(a.k.a log-binomial
regression)이라 불리는 이유를 알아 보자.
Risk, Relative Risk
Risk는 확률의 다른 말이다. 예를 들어 “오늘 벼락 맞을 확률”은 “오늘
벼락 맞을 Risk”로 바꿀 수 있다.
Relative Risk 또는 Risk Ratio(RR)은 상대위험도이다. 두 위험(확률)을
비교하기 위해 두 위험의 비를 계산한다. 예를 들어 “그냥 서 있을 때
벼락맞을 확률”과 “우산을 들고 서 있을 때 벼락맞을 확률”을 비교하기 위해
“우산을 들고 서 있을 때 벼락맞을 확률”을 “그냥 서 있을 때 벼락맞을
확률”으로 나눠주면 “우산을 들고 서 있을 때 벼락맞을 확률”이 “그냥 서
있을 때”보다 몇 배 더 많은지 알 수 있다. 이를 Relative
Risk 또는 Risk Ratio라고 한다.
Relative Risk Regression의 계수는 Relative Risk를 나타낸다.
위의 식에서 \(x_1\) 을 우산을 들고
서 있음을 나타내는 지시변수라고 하자(우산을 들고 서 있다면 \(x_1=1\) , 그렇지 않다면 \(x_1 = 0\) ). 이 경우 \(x_1=1\) 일 때 \(y\) ( \(y_1\) )와 \(x_1=0\) 일 때의 \(y\) ( \(y_0\) )를 비교해보면,
\[{\log}(y_1) = \beta_0 +
\beta_1\]
\[\Rightarrow y_1 = \exp(\beta_0 +
\beta_1)\]
\[\Rightarrow y_1 = \exp(\beta_0)
\exp(\beta_1)\]
\[y_0 = \exp(\beta_0)\]
위의 식에서 \(y_1\) (우산을 들고
있을 때 벼락맞을 확률) 에 \(y_0\)
(우산을 들고 있지 않을 때 벼락맞을 확률)를 나눠주면 \(y_1/y_0 = \exp(\beta_1)\) 으로 \(\beta_1\) 은 두 위험(확률)의 비를 나타내고
있음을 확인할 수 있다.
MLE 추정의 어려움
Relative Risk Regression은 모수에 제약이 있다. 예를 들어 \(y_0\) 과 \(y_1\) 은 모두 확률을 나타내기 때문에 0보다
작을 수 없고, 1보다 클 수 없다. 따라서 \[{\log}(y_1) = \beta_0 + \beta_1\] 의
우측은 양수가 될 수 없다. 그리고 \(\beta_0 +
beta_1\) 이 0과 매우 가까운 경우처럼, 모수가 제약 조건의 경계면에
위치하게 되면 모수의 분산을 추정하기가 까다롭다.
R에서 Relative Risk Regression
R에서 Relative Risk Regression을 하는 방법은 glm()
함수를 사용할 수도 있지만, logbin 패키지를 사용하면 좀더
안정적으로 추정을 할 수 있다고 한다. logbin::logbin()
함수의 method=를 사용하여 여러 가지 추정방법 중 하나를
선택하여 추정할 수도 있다.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(ggplot2)
library(logbin)
# from https://github.com/mdonoghoe/logbin
require(glm2, quietly = TRUE)
data(heart)
start.p <- sum(heart$Deaths) / sum(heart$Patients)
t.glm <- system.time(
fit.glm <- logbin(cbind(Deaths, Patients-Deaths) ~ factor(AgeGroup) + factor(Severity) +
factor(Delay) + factor(Region), data = heart,
start = c(log(start.p), -rep(1e-4, 8)), method = "glm", maxit = 10000)
)
## Warning: step size truncated due to divergence
## Warning: step size truncated due to divergence
## ...위의 코드는 logbin 패키지의 저자 github에서 가져온
테스트용 코드이다. 나머지 코드를 모두 실행하여 서로 다른 추정 방법을
사용해본 결과는 다음과 같다.
t.glm2 <- system.time(fit.glm2 <- update(fit.glm, method='glm2'))## Warning: step size truncated due to divergence
## Warning: step size truncated due to divergence
## Warning: step size truncated due to divergence
## Warning: step size truncated due to divergence## Warning: step size truncated due to increasing deviance## Warning: step size truncated due to divergence## Warning: step size truncated due to increasing deviance
## Warning: step size truncated due to increasing deviance
## Warning: step size truncated due to increasing deviance
## Warning: step size truncated due to increasing deviance
## Warning: step size truncated due to increasing deviancet.cem <- system.time(fit.cem <- update(fit.glm, method = "cem"))## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "packedMatrix" of class "replValueSp"; definition not updated## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "packedMatrix" of class "mMatrix"; definition not updated## Loading required package: doParallel## Loading required package: foreach## Loading required package: iterators## Loading required package: parallel## Loading required package: numDeriv## Loading required package: quantreg## Loading required package: SparseM##
## Attaching package: 'SparseM'## The following object is masked from 'package:base':
##
## backsolve##
## Attaching package: 'turboEM'## The following objects are masked from 'package:numDeriv':
##
## grad, hessiant.em <- system.time(fit.em <- update(fit.glm, method = "em"))
t.cem.acc <- system.time(fit.cem.acc <- update(fit.cem, accelerate = "squarem"))
t.em.acc <- system.time(fit.em.acc <- update(fit.em, accelerate = "squarem"))
objs = list("glm"=fit.glm,
"glm2"=fit.glm2,
"cem"=fit.cem,
"em"=fit.em,
"cem.acc" = fit.cem.acc,
"em.acc" = fit.em.acc)
params = c('converged', "loglik", "iter")
to_dataframe = function(objs, params) {
#param = params[1]
#obj[[param]]
dat = data.frame(model=names(objs))
for (param in params) {
dat[[param]] = sapply(objs,
function(x)
x[[param]])
}
return(dat)
}
dat = to_dataframe(objs, params)
dat$time = c(t.glm['elapsed'],
t.glm2['elapsed'],
t.cem['elapsed'],
t.em['elapsed'],
t.cem.acc['elapsed'],
t.em.acc['elapsed'])
print(dat)## model converged loglik iter time
## 1 glm FALSE -186.7366 10000 1.66
## 2 glm2 TRUE -179.9016 14 0.00
## 3 cem TRUE -179.9016 223196, 8451 47.31
## 4 em TRUE -179.9016 6492 3.13
## 5 cem.acc TRUE -179.9016 4215, 114 4.05
## 6 em.acc TRUE -179.9016 81 0.09저자는 cem이 가장 안정적으로 추정한다고 했으나 속도는
가장 느렸다. 가장 빠른 방법은 저자가 비교 대상에서 제외한
glm2이었다.
Relative-Risk vs Logistic
그렇다면 동일하게 이항반응을 모형화하는 logistic regression과는 어떤
차이가 있을까? Relative Risk는 \(\log(p)\) ( \(p\) : 확률)에서 선형이며, logistic은 \(\textrm{logit}(p)\) 에서 선형적이다.
예를 들어 다음의 모형을 보자.
\[\textrm{logit}(p) = \beta_0 + \beta_1
x_1 + \beta_2 x_2\]
\(x_1\) 과 \(x_2\) 을 모두 이항 변수라고 할 때 \(p|x_1=1, x_2=1\) 은 \(\textrm{logit}(p)\) 가 \(\beta_0 + \beta_1 + \beta_2\) 가 되어야
한다. 반면 아래의 모형에서는 \(p|x_1=1,
x_2=1\) 은 \(\log(p)\) 가 \(\beta_0 + \beta_1 + \beta_2\) 가 되어야
한다.
\[\log(p) = \beta_0 + \beta_1 x_1 +
\beta_2 x_2\]
예를 들어 \(p_{x_1=0, x_2=0} = p_{00} =
0.3\) , \(p_{10} = 0.4\) , \(p_{01}=0.5\) 인 경우를 생각해보자.
상호작용이 없는 로짓 모형을 사용한다면 \(\beta_0\) , \(\beta_1\), \(\beta_2\) 는 각각 -0.8472979, 0.4418328,
0.8472979이고, \(\textrm{logit}(p_{11}) =
-0.8472979 +0.4418328 + 0.8472979\) 이 되어야 한다. 계산해보면
\(p_{11} = 0.6086957\) 이 된다.
반면 상호작용이 없는 로그-바이노미얼(Relative Risk) 모형을
사용한다면, \(\beta_0\) , \(\beta_1\), \(\beta_2\) 는 각각 -1.203973, 0.2876821,
0.5108256이고, \(\textrm{log}(p_{11}) =
-1.203973 + 0.2876821 + 0.5108256\) 이 되어야 한다. 계산해보면
\(p_{11} = 0.6666665\) 이 된다.
따라서 상호작용이란 상대적이라고 말할 수 있다. logit 상에서
상호작용이 없다는 것과 log 상에서 상호작용이 없다는 것이 같지 않다.
relative risk model에서 상호작용이 없다면 \(x_2=1\) 의 효과가 \(x_1=0\) 에서도 \(x_1=1\) 에서도 동일하게 relative
risk=exp(0.51) (1.6배)로 나타난다. 하지만 동일한 결과가 logit 상에서는
상호작용이 존재한다!
Relative-Risk vs Logistic 2: Variance Estimation
몇몇 논문은 Relative-Risk의 모수 se가 더 정확하기 때문에 Logistic
regression 대신 Relative Risk Regression을 쓰길 권장한다. 과연 그럴까?
실험을 해보자.
dat = data.frame(y= c(rep(T,4), rep(F,16)))
# logistic regression
fit1 = glm(y ~ 1, family=binomial,
data=dat)
# relative-risk regression
fit2 = logbin(y ~ 1,
data = dat,
method = "cem",
maxit = 10000)
fit1$converged## [1] TRUEfit2$converged## [1] TRUE추정확률은 모두 같다.
logistic = function(x) 1/(1+exp(-x))
logistic(coef(fit1))## (Intercept)
## 0.2exp(coef(fit2))## (Intercept)
## 0.2하지만 95%-신뢰구간은 다르다.
logistic(confint(fit1))## Waiting for profiling to be done...## 2.5 % 97.5 %
## 0.0668305 0.4053389exp(confint(fit2))## 2.5 % 97.5 %
## (Intercept) 0.08324556 0.4805061prop = prop.test(4, 16)
binom = binom.test(4, 16)
prop$conf.int## [1] 0.08332776 0.52591618
## attr(,"conf.level")
## [1] 0.95binom$conf.int## [1] 0.07266204 0.52377082
## attr(,"conf.level")
## [1] 0.95dat = data.frame(p = c('p0.025', 'p0.975'),
logistic = logistic(confint(fit1)),
rr = exp(confint(fit2)['(Intercept)', ]),
prop = prop$conf.int,
binom = binom$conf.int)## Waiting for profiling to be done...dat %>%
gather('key', 'value', logistic:binom) %>%
spread(key='p', value='value') %>%
mutate(key=factor(key, levels=c('logistic', 'rr', 'prop', 'binom'))) %>%
ggplot(aes(x=key)) +
geom_errorbar(aes(ymin=p0.025, ymax=p0.975), width=0.2) +
ylim(c(0,1)) +
geom_hline(yintercept=0.2)## Warning: attributes are not identical across measure variables;
## they will be dropped
Leave a comment