疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

Rでglmを用いてオッズ比を算出するコマンド

今回はRを用いた場合のオッズ比の求め方を記述しておく。
(オッズ比が{e^{回帰式の推定値}}で求めることができることを知っていれば、簡単に理解できる。)

はじめに

医学研究を行う上で、なんらかの指標で群分けしたグループでのある事象のおこりやすさ(オッズ)を比較することがある(オッズ比)。今回はRを用いた場合の簡単なORの求め方を記述する。
なお、模擬データは本記事の末尾に記述するので参照されたい(ランダムで乱数を発生させているので、同じ結果にならないことに注意。)

簡単にデータのプレビューをする。

head(Ndtnew)

#         age  gender   height   weight      BMI alcohol ID BMI01
# 837 44.54117 Female 152.4315 42.94982 18.48467     low  1     0
# 730 37.54635 Female 153.3675 48.91365 20.79522     low  2     0
# 290 33.72811   Male 155.6905 58.27677 24.04204    high  3     0
# 24  52.81553   Male 164.7295 77.22198 28.45760  middle  4     1
# 573 50.11663 Female 155.1740 36.71924 15.24950     low  5     0
# 582 47.65373 Female 141.8624 47.11554 23.41152     low  6     0

データフレームには年齢、性別、身長、体重、BMI、アルコールの摂取頻度、ID、BMIの01データ(BMI >25 →1, BMI <=25 →0)の項目が含まれている。

今回の解析ではアルコール低摂取群(low)に対するアルコール中程度摂取群(middle)、アルコール中程度摂取群(high)の性別、年齢で補正した肥満のORsを求める。

早速、ORの算出を行う。

1. まずはglmを用いて、ロジスティック回帰分析を行う。

obesity1.glm <- glm(BMI01~alcohol+age+gender, family=binomial,data=Ndtnew)

2. その結果を表示する。

# Call:
# glm(formula = BMI01 ~ alcohol + age + gender, family = binomial, 
#     data = Ndtnew)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -1.2800  -0.0900  -0.0840  -0.0741   3.3821  
# 
# Coefficients:
#                Estimate Std. Error z value Pr(>|z|)    
# (Intercept)    0.329512   0.500709   0.658    0.510    
# alcoholmiddle -0.248181   0.215683  -1.151    0.250    
# alcoholhigh    0.135045   0.289193   0.467    0.641    
# age           -0.007546   0.010460  -0.721    0.471    
# genderFemale  -5.553169   0.721649  -7.695 1.41e-14 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
#     Null deviance: 1000.80  on 999  degrees of freedom
# Residual deviance:  604.55  on 995  degrees of freedom
# AIC: 614.55
# 
# Number of Fisher Scoring iterations: 8

ここメイン
3. この結果を用いて、ORの算出をする。

ModelName <- obesity1.glm
alpha <- 0.05
x <- summary(ModelName)
y <- confint(ModelName, level=1-alpha)
OR <- exp(x$coefficients[,1])
LowerCL <- exp(y[,1])
UpperCL <- exp(y[,2])
pvalue <- summary(obesity1.glm)$coefficients[,"Pr(>|z|)"]
z <- cbind(OR, LowerCL, UpperCL, pvalue)
round(z, 5)

結果はこんな感じになる。

#                   OR LowerCL UpperCL  pvalue
# (Intercept)   1.39029 0.52132 3.72268 0.51048
# alcoholmiddle 0.78022 0.51060 1.19023 0.24987
# alcoholhigh   1.14459 0.64927 2.02241 0.64052
# age           0.99248 0.97227 1.01302 0.47068
# genderFemale  0.00388 0.00063 0.01247 0.00000

まとめ

ちなみにこのサンプルデータではアルコールの摂取頻度の間に有意な関連は見られなかった。
本来であれば、glmにアルコール摂取の連続量(0,1,2など)の変数を組み込んだ、傾向性の検討も実施すべきと考えられる。

サンプルデータ(Ndtnew)の作成

age <- runif(420, 30, 63)
gender <- rep("Male", 420)
a <- rep(c("L", "M", "H"), times = c(150, 200, 70))
alcohol <- as.factor(sample(a, length(a)))
height <- rnorm(420, 165.4, 6.2)
weight <- rnorm(420, 67.2, 8.9)
BMI <- weight/((height/100)^2)
NdtMen <- data.frame(age, gender, height, weight, BMI, alcohol)
age <- runif(580, 37, 68)
gender <- rep("Female", 580)
a <- rep(c("L", "M", "H"), times = c(400, 150, 30))
alcohol <- as.factor(sample(a, length(a)))
height <- rnorm(580, 154.4, 6.2)
weight <- rnorm(580, 43.4, 4.5)
BMI <- weight/((height/100)^2)
NdtFemale <- data.frame(age, gender, height, weight, BMI, alcohol)
Ndtnew <- rbind(NdtMen, NdtFemale)
ID <- 1:1000
ID <- sample(ID,1000,replace=FALSE)
Ndtnew <- transform(Ndtnew, ID=ID)
sortlist <- order(Ndtnew$ID)
Ndtnew <- Ndtnew[sortlist, ]
Ndtnew$BMI01[Ndtnew$BMI < 25] <- 0
Ndtnew$BMI01[Ndtnew$BMI >= 25] <- 1
Ndtnew$alcohol <- factor(Ndtnew$alcohol, level=c("L", "M", "H"), label=c("low", "middle", "high"))


20160202
RF