應(yīng)用統(tǒng)計學(xué)與R語言實現(xiàn)筆記（番外篇四）——bookdown使用與OR值計算

本期是之前做的應(yīng)用統(tǒng)計學(xué)與R語言實現(xiàn)筆記的番外篇四，本期主要關(guān)注兩個問題，一個是重新利用R的bookdown包創(chuàng)建新的電子書，另一個是計算公共衛(wèi)生當(dāng)中一個比較常見的指標(biāo)OR值。

1 bookdown使用

bookdown是謝益輝之前開發(fā)的R語言包，可以基于rmarkdown快速生成在線電子書，并且可以輸出pdf和epub。具體的使用方法可以參見官方文檔。

https://bookdown.org/yihui/bookdown/

這里由于中文在輸出pdf中容易有bug，因此中文圖書推薦使用謝益輝提供的模板進(jìn)行修改。同時可以參考李東風(fēng)的這本中文使用指南輔助進(jìn)行。

https://www.math.pku.edu.cn/teachers/lidf/docs/Rbook/html/_Rbook/

這里提供一些使用過程中的tips經(jīng)驗。

  • Latex公式在在$與公式間不要空格。

  • 插入圖片建議使用r的function。例子如下：


knitr::include_graphics("fig/fig1.jpg")

其中include_graphics括號后面為圖片路徑。同時在R code的設(shè)置里，在設(shè)置圖片大小時推薦使用out.width和out.height參數(shù)，設(shè)置為'100%'，這樣圖片可以根據(jù)排版情況進(jìn)行自適應(yīng)。

  • 如果想加入r代碼塊而不想運行僅作為展示的話，需要在R code的設(shè)置里設(shè)置為r eval=FALSE, echo = T。

  • 目前這個版本的封面圖片設(shè)置參數(shù)cover-image只能生成epub里的封面，pdf無法添加。

  • 默認(rèn)模板會生成圖和表目錄，不需要的可以在index.rmd的輸出設(shè)置里把lot和lof設(shè)置為false。

  • 默認(rèn)模板pdf里有一句話“獻(xiàn)給……呃，愛誰誰吧”，需要在模板的latex文件夾下的before_body.tex里去掉。

  • 默認(rèn)模板設(shè)置是B5的紙張大小，邊距設(shè)置也是左右不對稱。這個是在index.rmd的輸出設(shè)置里。實際上也是latex的設(shè)置?？筛鶕?jù)自己喜好做調(diào)整。


geometry: [b5paper, tmargin=2.5cm, bmargin=2.5cm, lmargin=3.5cm, rmargin=2.5cm]
  • 用github托管的話，可以在bookdown.yml文件里設(shè)置輸出文件夾參數(shù)，在最后一行添加參數(shù)（output_dir: "docs"）。然后在github的Pages設(shè)置對應(yīng)的根目錄。同時需要在R里輸入如下命令。讓網(wǎng)頁不使用默認(rèn)jekyll主題。

file.create('docs/.nojekyll')

最后奉上最新的bookdown在線電子書地址：

http://gisersqdai.top/Note-of-Applied-Statistics-with-R-Book/

2 公式更正

在修改的過程里，我發(fā)現(xiàn)了來自BruceZhaoR同學(xué)18年的一條pr，雖然不知道什么原因我一直沒注意到這條pr。這里鄭重向這位同學(xué)道歉，非常感謝他的指正。他指出在原本第三章描述性統(tǒng)計里的樣本方差與標(biāo)準(zhǔn)差公式里有誤。并給出了wiki上的參考公式。

image

wiki:https://en.wikipedia.org/wiki/Standard_deviation#Corrected_sample_standard_deviation

具體錯誤這里也說明下。原公式如下：

樣本方差：

s^2=\frac{\sum_{i=1}^N (x_i-\mu)^2}{n-1} 或 s^2=\frac{\sum_{i=1}^k (M_i-\mu)^2f_i}{n-1}

樣本標(biāo)準(zhǔn)差：

s=\sqrt {\frac{\sum_{i=1}^N (x_i-\mu)^2}{n-1}} 或 s=\sqrt{\frac{\sum_{i=1}^k (M_i-\mu)^2f_i}{n-1}}

這里用\mu是不對的，\mu雖然可以指代統(tǒng)計學(xué)中的均值，但是\mu是代表總體均值。而嚴(yán)格來說，樣本均值通常只是近似總體均值，因此必須作區(qū)分，故常用\bar x來做為樣本均值。故修改后公式為

樣本方差：

s^2=\frac{\sum_{i=1}^N (x_i-\bar x)^2}{n-1} 或 s^2=\frac{\sum_{i=1}^k (M_i-\bar x)^2f_i}{n-1}

樣本標(biāo)準(zhǔn)差：

s=\sqrt {\frac{\sum_{i=1}^N (x_i-\bar x)^2}{n-1}} 或 s=\sqrt{\frac{\sum_{i=1}^k (M_i-\bar x)^2f_i}{n-1}}

3 OR值計算

由于我目前主要從事健康地理學(xué)方面的研究，最近碰上了一個基礎(chǔ)的OR值計算問題。首先OR值的全稱是odds ratio值，這是公共衛(wèi)生領(lǐng)域的一個專業(yè)名詞。這里給出Encyclopedia of Public Health的定義。

The odds ratio (OR) provides a measure of the strength of relationship between two variables, most commonly an exposure and a dichotomous outcome. It is most commonly used in a case control study where it is defined as “the ratio of the odds of being exposed in the group with the outcome to the odds of being exposed in the group without the outcome.”

This concept can be extended to a situation with multiple levels of exposure (e.g., low, moderate, or high exposure to an environmental containment). One exposure level is assigned as the “reference” level. For each of the remaining exposure levels, one divides the odds of that exposure level in the outcome positive group (compared with the

reference level) by the odds of that exposure level in the outcome negative group.

The OR ranges in value from 0 to infinity. Values close to 1.0 indicate no relationship between the exposure and the outcome. Values less than 1.0 suggest a protective effect, while values greater than 1.0 suggest a causative or adverse effect of exposure.

這里簡單翻譯一下，OR值是用來度量兩個變量之間關(guān)系強度的指標(biāo)，常見于暴露水平與二分的健康結(jié)局變量。最常用在案例對照研究。OR被定義為“組中暴露患病幾率與暴露未患病幾率的比值”。這個概念通?？梢酝卣沟蕉嗨奖┞吨笜?biāo)，通常定義某一類別的暴露水平為參考水平，對于剩余的暴露水平，則除以該參考水平的暴露幾率用來進(jìn)行比較。這里要先給出odds的定義。odds，稱為幾率、比值、比數(shù)，是指某事件發(fā)生的可能性(概率)與不發(fā)生的可能性（概率）之比。用p表示事件發(fā)生的概率，則：odds = p/(1-p)。OR值的公共衛(wèi)生意義如下：范圍從0到無窮大。接近1.0的值表示暴露與結(jié)果之間沒有關(guān)系。小于1.0的值表示保護(hù)作用，而大于1.0的值表示暴露的致病性或不利影響。

針對一個標(biāo)準(zhǔn)的2x2流病表格（如下）。實際上OR值計算如下：

暴露時患病幾率=\frac{暴露時患病病例數(shù)}{未暴露時患病病例數(shù)}=\frac{a}{c}

暴露時未患病幾率=\frac{暴露時未患病病例數(shù)}{未暴露時未患病病例數(shù)}=\fracu0z1t8os

OR = \frac{a/c}{b/d}=\frac{ad}{bc}

| | Outcome +ve | Outcome -ve |

| :-------------------------- | :---------- | :---------------|

| Exposed | a | b |

| No exposed | c | d |

這里選用一個ICU的數(shù)據(jù)來進(jìn)行說明，這個數(shù)據(jù)來源于David W. Hosmer等人出版的applied logistic regression一書中的數(shù)據(jù)。獲取途徑可以通過安裝這個書的r包，命令如下。


install.packages('aplore3')

安裝完成以后，載入數(shù)據(jù)做個初步探索。


library(aplore3)

data(icu)

head(icu)
image

為了簡單化，我們這里定義健康結(jié)局變量為status，即數(shù)據(jù)中的sta（Lived或者Died），感興趣的自變量為age。首先繪制一下圖。由于status是個二分變量。所以圖就變成了如下的樣子。

image

你是否覺得很熟悉？其實這就是logistic regression的典型數(shù)據(jù)。那么根據(jù)age的數(shù)據(jù)，我們做一個處理，統(tǒng)計不同年齡段的死亡率，以10歲為分界線。我們可以得到如下的圖。

image

那么我們突然發(fā)現(xiàn)，這個散點是有線性趨勢的。假設(shè)我們采用線性回歸來做分析，即假定有：pr(death)=\beta_0+\beta_1(age)，不就可以擬合了嗎？但是我們又會發(fā)現(xiàn)一個問題。那就是這里的y(pr(death))是有現(xiàn)實意義的實數(shù)，也就是它的值域必須在（0,1）中。然而等式右邊實際上是可以取任何值的（根據(jù)\beta_0 , \beta_1, age），因此這個線性方程即使求解出來，預(yù)測值通常會超過實際的值域。所以為了解決這個問題，logistics regression就提出了。首先是定義了logit函數(shù)為：

logit(p)=log(\frac{p}{1-p})

p=pr(death)

那么這個logit函數(shù)的現(xiàn)實意義是事件發(fā)生幾率的對數(shù)。那么同時模型就變成了：

log(\frac{p}{1-p})=\beta_0+\beta_1(age)

這時我們就會發(fā)現(xiàn)p的值域是在（0,1），而logit(p)的值域則是[ - \infty, + \infty ]

image

那這個時候我們就可以用線性回歸方法求解系數(shù)了，因此logistic regression也被稱為廣義線性回歸的一類。

那么我們再來看一個更特殊的情況，就是前面說的2x2聯(lián)表的情況。假定自變量是個分類變量。這里選用icu數(shù)據(jù)里的type來分析（健康結(jié)局變量不變）。也就是說方程如下：

logit(p)=\beta_0+\beta_1(I_{type})

2x2聯(lián)表則為：

| | Lived | Died |

| :-------------------------- | :---------- | :---------------|

| elective admission | a | b |

| emeregency admission | c | d |

那么這時候I_{type}=0時是elective admission，I_{type}=1時是emeregency admission。因此我們可以得到對應(yīng)的y值。也就是elective adminssion的logit(p)為\beta_0。而emergency admission的logit(p)為\beta_0+\beta_1。那么根據(jù)logit函數(shù)的定義，我們就有如下的式子：

對elective adminssion的odds：

odds_{ele}=\frac{p}{1-p}=\frac{a}=e^{\beta_0}

對emergency adminssion的odds：

odds_{eme}=\frac{p}{1-p}=\frac{c}u0z1t8os=e^{\beta_0+\beta_1}

那么所以這個OR值就可以計算：

OR =\frac{ad}{bc} =odds_{ele}/odds_{eme}=\frac{a}/\frac{c}u0z1t8os=e^{\beta_0+\beta_1}/e^{\beta_0}=e^{\beta_0+\beta_1-\beta_0}=e^{\beta_1}

也就是說，實際上OR值就是e的logistics regression中的回歸系數(shù)次方。因此通常在公共衛(wèi)生研究中求OR值，第一步就是先做logistics regression然后接著進(jìn)行計算。對應(yīng)其實也可以計算OR的95%置信區(qū)間以及p值(Explaining Odds Ratios)。這塊這里就不詳述了。我這里還是采用icu數(shù)據(jù)做個樣例分析，展示三種方式（第一種不借助其他包，第二個使用epiDisplay，第三個是用questionr）求OR。


## 1 without any packages

modellogit <- glm(sta~type, data = icu, family = binomial)

ORDF <- data.frame(exp(cbind(OR = coef(modellogit), confint(modellogit))))

ORDF



## 2 Using epiDisplay

install.packages('epiDisplay')

library(epiDisplay)

modellogit <- glm(sta~type, data = icu, family = binomial)

ORDF <- logistic.display(modellogit)

ORDF



## 3 Using questionr

install.packages('questionr')

library(questionr)

modellogit <- glm(sta~type, data = icu, family = binomial)

ORDF <- odds.ratio(modellogit)

ORDF

目前個人推薦第三種，能把p值一起算出來，這里要注意R里面默認(rèn)factor的第一個因子作為參考組（baseline），如需要設(shè)置不同的參考組?？梢杂萌缦碌暮瘮?shù)。

image

icu$type <- relevel(icu$type, ref = "emergency")
image

最后本次的代碼也都是在之前的github項目上。歡迎大家使用。最后再放一下兩個項目地址：

Note-of-Applied-Statistics-with-R

Note-of-Applied-Statistics-with-R-Book

參考鏈接:

  1. Logistic regression for a Yes/No outcome
  1. Confused with the reference level in logistic regression in R
  1. Logistic regression - defining reference level in R
  1. Reference category and interpreting regression coefficients in R
  1. Confidence Intervals for RRs, ORs in R
  1. How to calculate the p.value of an odds ratio in R?
  1. Calculating Odds Ratio in R
  1. Odds ratios and logistic regression: further examples of their use and interpretation
  1. BMI 541/699 Lecture 22
  1. aplore3
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容