R語言卡方檢驗大全

本文首發(fā)于公眾號：醫(yī)學(xué)和生信筆記

醫(yī)學(xué)和生信筆記，專注R語言在臨床醫(yī)學(xué)中的使用确虱，R語言數(shù)據(jù)分析和可視化。主要分享R語言做醫(yī)學(xué)統(tǒng)計學(xué)替裆、meta分析校辩、網(wǎng)絡(luò)藥理學(xué)窘问、臨床預(yù)測模型、機(jī)器學(xué)習(xí)宜咒、生物信息學(xué)等惠赫。

卡方檢驗/列聯(lián)表資料的卡方檢驗在臨床中非常常見！

因為最近又有一批臨床數(shù)據(jù)要進(jìn)行統(tǒng)計故黑，所以趁機(jī)把卡方檢驗的R語言實現(xiàn)再重新梳理一遍儿咱。

這篇文章涵蓋了孫振球，徐勇勇《醫(yī)學(xué)統(tǒng)計學(xué)》第4版 卡方檢驗章節(jié) 中的 所有內(nèi)容场晶。課本電子版和配套數(shù)據(jù)已上傳到QQ群混埠，需要的朋友加群下載即可。

image.png

不同類型卡方檢驗的選擇四格表資料的卡方檢驗方法1方法2配對四格表資料的卡方檢驗四格表資料的 Fisher 確切概率法行 x 列表資料的卡方檢驗多個樣本率的比較樣本構(gòu)成比的比較雙向無序分類資料的關(guān)聯(lián)性檢驗雙向有序分組資料的線性趨勢檢驗多個樣本率間的多重比較Cochran-Mantel-Haenszel 卡方統(tǒng)計量檢驗頻數(shù)分布擬合優(yōu)度卡方檢驗

不同類型卡方檢驗的選擇

課本中關(guān)于四格表資料的卡方檢驗的方法選擇以及R x C表資料的檢驗方法選擇做了非常好的總結(jié)峰搪，在這里一并和大家分享一下：

四格表資料的方法選擇：

  • 當(dāng) n(樣本量)≥40 且所有的T(期望頻數(shù))≥5時岔冀，用χ2檢驗的基本公式或四格表資料之χ2檢驗的專用公式；當(dāng)P ≈ α?xí)r概耻，改用四格表資料的 Fisher 確切概率法使套；
  • 當(dāng) n≥40 但有 1≤T<5 時，用四格表資料χ2檢驗的校正公式鞠柄，或改用四格表資料的 Fisher 確切概率法侦高。
  • 當(dāng) n<40，或 T<1時厌杜，用四格表資料的 Fisher 確切概率法奉呛。

R×C表資料的分類及其檢驗方法的選擇：

R×C表資料可以分為雙向無序、單向有序夯尽、雙向有序?qū)傩韵嗤碗p向有序?qū)傩圆煌?類瞧壮。

  1. 雙向無序R×C表資料 R×C表資料中兩個分類變量皆為無序分類變量對于該類資料，若研究目的為多個樣本率（或構(gòu)成比）的比較匙握，可用行×列表資料的χ2檢驗：若研究目的為分析兩個分類變量之間有無關(guān)聯(lián)性以及關(guān)系的密切程度時泳秀，可用行×列表資料的χ2檢驗以及Pearson列聯(lián)系數(shù)進(jìn)行分析押框。
  2. 單向有序R×C表資料 有兩種形式凉翻。一種是R×C表資料中的分組變量（如年齡）是有序的默赂，而指標(biāo)變量（如傳染病的類型）是無序的。其研究目的通常是分析不同年齡組各種傳染病的構(gòu)成情況蛾娶，此種單向有序R×C表資料可用行×列表資料的χ2檢驗進(jìn)行分析灯谣。另一種情況是R×C表資料中的分組變量
    (如療法)為無序的，而指標(biāo)變量（如療效按等級分組）是有序的蛔琅。其研究目的為比較不同療法的療效胎许，此種單向有序R×C表資料宜用秩轉(zhuǎn)換的非參數(shù)檢驗進(jìn)行分析。
  3. 雙向有序?qū)傩韵嗤腞×C表資料 R×C表資料中的兩個分類變量皆為有序且屬性相同。實際上是配對四格表資料的擴(kuò)展呐萨，即水平數(shù)≥3的配伍資料杀饵，如用兩種檢測方法同時對同一批樣品的測定結(jié)果。其研究目的通常是分析兩種檢測方法的一致性谬擦，此時宜用一致性檢驗或稱Kappa檢驗切距；也可用特殊模型分
    析方法（可用SAS軟件）。
  4. 雙向有序?qū)傩圆煌腞×C表資料 R×C表資料中兩個分類變量皆為有序的惨远，但屬性不同谜悟。對于該類資料，若研究目的為分析不同年齡組患者療效之間有無差別時北秽，可把它視為單向有序R×C表資料葡幸，選用秩轉(zhuǎn)換的非參數(shù)檢驗；若研究目的為分析兩個有序分類變量間是否存在相關(guān)關(guān)系贺氓，宜用等級相關(guān)分析：若研究目的為分析兩個有序分類變量間是否存在線性變化趨勢蔚叨，宜用前述的雙向有序分組資料的線性趨勢檢驗(test for linear trend)。

四格表資料的卡方檢驗

使用課本例7-1的數(shù)據(jù)辙培。

首先是構(gòu)造數(shù)據(jù)蔑水，本次數(shù)據(jù)自己從書上摘錄。扬蕊。

ID<-seq(1,200)
treat<-c(rep("treated",104),rep("placebo",96))
treat<- factor(treat)
impro<-c(rep("marked",99),rep("none",5),rep("marked",75),rep("none",21))
impro<-as.factor(impro)
data1<-data.frame(ID,treat,impro)

str(data1)
## 'data.frame':    200 obs. of  3 variables:
##  $ ID   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ treat: Factor w/ 2 levels "placebo","treated": 2 2 2 2 2 2 2 2 2 2 ...
##  $ impro: Factor w/ 2 levels "marked","none": 1 1 1 1 1 1 1 1 1 1 ...

數(shù)據(jù)一共3列搀别，第一列是id，第二列是治療方法尾抑，第三列是等級（有效和無效）歇父。

簡單看下各組情況：

table(data1$treat,data1$impro)
##          
##           marked none
##   placebo     75   21
##   treated     99    5

做卡方檢驗有2種方法，分別演示：

方法1

直接使用 gmodels包里面的 CrossTable()函數(shù)再愈，非常強(qiáng)大榜苫，直接給出所有結(jié)果，和SPSS差不多翎冲。

library(gmodels)

CrossTable(data1$treat, data1$impro, digits = 4, expected = T, chisq = T, fisher = T, mcnemar = T, format = "SPSS")
## 
##    Cell Contents
## |-------------------------|
## |                   Count |
## |         Expected Values |
## | Chi-square contribution |
## |             Row Percent |
## |          Column Percent |
## |           Total Percent |
## |-------------------------|
## 
## Total Observations in Table:  200 
## 
##              | data1$impro 
##  data1$treat |   marked  |     none  | Row Total | 
## -------------|-----------|-----------|-----------|
##      placebo |       75  |       21  |       96  | 
##              |  83.5200  |  12.4800  |           | 
##              |   0.8691  |   5.8165  |           | 
##              |  78.1250% |  21.8750% |  48.0000% | 
##              |  43.1034% |  80.7692% |           | 
##              |  37.5000% |  10.5000% |           | 
## -------------|-----------|-----------|-----------|
##      treated |       99  |        5  |      104  | 
##              |  90.4800  |  13.5200  |           | 
##              |   0.8023  |   5.3691  |           | 
##              |  95.1923% |   4.8077% |  52.0000% | 
##              |  56.8966% |  19.2308% |           | 
##              |  49.5000% |   2.5000% |           | 
## -------------|-----------|-----------|-----------|
## Column Total |      174  |       26  |      200  | 
##              |  87.0000% |  13.0000% |           | 
## -------------|-----------|-----------|-----------|
## 
##  
## Statistics for All Table Factors
## 
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  12.85707     d.f. =  1     p =  0.0003362066 
## 
## Pearson's Chi-squared test with Yates' continuity correction 
## ------------------------------------------------------------
## Chi^2 =  11.3923     d.f. =  1     p =  0.0007374901 
## 
##  
## McNemar's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  50.7     d.f. =  1     p =  1.076196e-12 
## 
## McNemar's Chi-squared test with continuity correction 
## ------------------------------------------------------------
## Chi^2 =  49.40833     d.f. =  1     p =  2.078608e-12 
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio:  0.1818332 
## 
## Alternative hypothesis: true odds ratio is not equal to 1
## p =  0.0005286933 
## 95% confidence interval:  0.05117986 0.5256375 
## 
## Alternative hypothesis: true odds ratio is less than 1
## p =  0.0002823226 
## 95% confidence interval:  0 0.4569031 
## 
## Alternative hypothesis: true odds ratio is greater than 1
## p =  0.9999541 
## 95% confidence interval:  0.06281418 Inf 
## 
## 
##  
##        Minimum expected frequency: 12.48

可以看到這個函數(shù)直接給出所有結(jié)果单刁，根據(jù)需要自己選擇合適的即可。

本例符合pearson卡方府适，卡方值為12.85707，p<0.01肺樟，和課本一致檐春。

方法2

先把數(shù)據(jù)變成2x2列聯(lián)表，然后用 chisq.test函數(shù)做.

mytable <- table(data1$treat,data1$impro)

mytable
##          
##           marked none
##   placebo     75   21
##   treated     99    5

chisq.test(mytable,correct = F) # 和SPSS一樣
## 
##  Pearson's Chi-squared test
## 
## data:  mytable
## X-squared = 12.857, df = 1, p-value = 0.0003362

這個結(jié)果和課本也是一致的么伯，和SPSS算出來的也是一樣的疟暖。

四格表資料卡方檢驗的專用公式/四格表資料卡方檢驗的校正公式/配對四格表資料的卡方檢驗/四格表資料的Fisher精確概率法，都可以用方法1可直接解決。

下面使用R語言自帶的chisq.test()函數(shù)進(jìn)行演示俐巴。

使用課本例7-2的數(shù)據(jù)骨望，這是一個連續(xù)校正卡方檢驗。

per <- matrix(c(46,6,18,8),
              nrow = 2, byrow = T,
              dimnames = list(group = c("胞磷膽堿","神經(jīng)節(jié)苷脂"),
                              effect = c("有效","無效")
                              )
              )

per
##             effect
## group        有效 無效
##   胞磷膽堿     46    6
##   神經(jīng)節(jié)苷脂   18    8

進(jìn)行連續(xù)校正的卡方檢驗：

chisq.test(per, correct = T)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  per
## X-squared = 3.1448, df = 1, p-value = 0.07617

配對四格表資料的卡方檢驗

使用課本例7-3的數(shù)據(jù)欣舵。

ana <- matrix(c(11,12,2,33), nrow = 2, byrow = T,
              dimnames = list(免疫熒光 = c("陽性","陰性"),
                              乳膠凝集 = c("陽性","陰性")
                              )
              )

ana
##         乳膠凝集
## 免疫熒光 陽性 陰性
##     陽性   11   12
##     陰性    2   33

配對四個表資料需要用McNemar檢驗：

mcnemar.test(ana, correct = T)
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  ana
## McNemar's chi-squared = 5.7857, df = 1, p-value = 0.01616

四格表資料的 Fisher 確切概率法

使用課本例7-4的數(shù)據(jù)擎鸠。

hbv <- matrix(c(4,18,5,6), nrow = 2, byrow = T,
              dimnames = list(組別 = c("預(yù)防注射組","非預(yù)防組"),
                              效果 = c("陽性","陰性")
                              )
              )
hbv
##             效果
## 組別       陽性 陰性
##   預(yù)防注射組    4   18
##   非預(yù)防組      5    6

進(jìn)行 Fisher 檢驗：

fisher.test(hbv)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  hbv
## p-value = 0.121
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.03974151 1.76726409
## sample estimates:
## odds ratio 
##  0.2791061

P值為0.121，和課本一樣缘圈。

行 x 列表資料的卡方檢驗

行 x 列表資料的卡方檢驗有很多種情況劣光，不是所有的列聯(lián)表資料都可以直接用卡方檢驗，大家要注意甄別糟把！方法選擇可以參考本篇開頭部分绢涡。

多個樣本率的比較

使用課本例7-6的數(shù)據(jù)。

首先是構(gòu)造數(shù)據(jù)遣疯，本次數(shù)據(jù)直接讀取雄可，也可以自己手動摘錄。

df <- foreign::read.spss("./例07-06物理藥物外用膏用.sav", to.data.frame = T,reencode = "utf-8")
## re-encoding from utf-8

str(df)
## 'data.frame':    6 obs. of  3 variables:
##  $ .?.: Factor w/ 3 levels ".....?...","?........",..: 1 1 2 2 3 3
##  $ ..Ч: Factor w/ 2 levels "..Ч","..Ч_duplicated_2": 1 2 1 2 1 2
##  $ f  : num  199 7 164 18 118 26
##  - attr(*, "variable.labels")= Named chr(0) 
##   ..- attr(*, "names")= chr(0)

數(shù)據(jù)一共3列缠犀，第1列是療法数苫，第2列是有效無效，第3列是頻數(shù).

進(jìn)行 行 x 列表資料的卡方檢驗夭坪，首先要對數(shù)據(jù)格式轉(zhuǎn)換一下文判，變成 table或者 矩陣：

M <- matrix(df$f,nrow = 3,byrow = T,
            dimnames = list(trt = c("物理", "藥物", "外用"),
                            effect = c("有效","無效")))

M
##       effect
## trt    有效 無效
##   物理  199    7
##   藥物  164   18
##   外用  118   26

這里教大家一個可視化列聯(lián)表資料非常好用的馬賽克圖：

mosaicplot(M)
image.png

進(jìn)行 行 x 列表資料的卡方檢驗：

kf <- chisq.test(M, correct = F)

kf
## 
##  Pearson's Chi-squared test
## 
## data:  M
## X-squared = 21.038, df = 2, p-value = 2.702e-05

多個樣本率的比較也可以使用以下函數(shù)進(jìn)行檢驗：

# 只適用于兩列的，類似于 有效/無效 這種室梅！
prop.test(M, correct = TRUE)
## 
##  3-sample test for equality of proportions without continuity correction
## 
## data:  M
## X-squared = 21.038, df = 2, p-value = 2.702e-05
## alternative hypothesis: two.sided
## sample estimates:
##    prop 1    prop 2    prop 3 
## 0.9660194 0.9010989 0.8194444

可以看到兩種結(jié)果是一樣的戏仓，和課本一致的！

樣本構(gòu)成比的比較

使用課本例7-7的數(shù)據(jù)亡鼠。

ace <- matrix(c(42,48,21,30,72,36),nrow = 2,byrow = T,
              dimnames = list(dn = c("dn組","非dn組"),
                              idi = c("dd","id","ii")
                              )
              )
ace
##         idi
## dn       dd id ii
##   dn組   42 48 21
##   非dn組 30 72 36

進(jìn)行卡方檢驗：

chisq.test(ace, correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  ace
## X-squared = 7.9127, df = 2, p-value = 0.01913

卡方值為7.91赏殃，和課本一致。

雙向無序分類資料的關(guān)聯(lián)性檢驗

使用課本例例7-8的數(shù)據(jù)间涵。

blood <- matrix(c(431,490,902,388,410,800,495,587,950,137,179,32),
                nrow = 4,byrow = T,
                dimnames = list(abo = c("o","a","b","ab"),
                                mn = c("m","n","mn")
                                )
                )
blood
##     mn
## abo    m   n  mn
##   o  431 490 902
##   a  388 410 800
##   b  495 587 950
##   ab 137 179  32

進(jìn)行 關(guān)聯(lián)性檢驗：

chisq.test(blood,correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  blood
## X-squared = 213.16, df = 6, p-value < 2.2e-16

計算列聯(lián)系數(shù)：

library(vcd)
## Loading required package: grid

assocstats(blood)
##                     X^2 df P(> X^2)
## Likelihood Ratio 248.14  6        0
## Pearson          213.16  6        0
## 
## Phi-Coefficient   : NA 
## Contingency Coeff.: 0.188 
## Cramer's V        : 0.136

Pearson列聯(lián)系數(shù)是0.188仁热，和課本一樣。

雙向有序分組資料的線性趨勢檢驗

使用課本例7-9的數(shù)據(jù)勾哩。

ather <- matrix(c(70,22,4,2,27,24,9,3,16,23,13,7,9,20,15,14),
                nrow = 4,byrow = T,
                dimnames = list(age = c("20~","30~","40~","≥50"),
                                level = c("-","+","++","+++")
                                )
                )
ather
##      level
## age    -  + ++ +++
##   20~ 70 22  4   2
##   30~ 27 24  9   3
##   40~ 16 23 13   7
##   ≥50  9 20 15  14

進(jìn)行卡方檢驗：

chisq.test(ather)
## 
##  Pearson's Chi-squared test
## 
## data:  ather
## X-squared = 71.432, df = 9, p-value = 7.97e-12

但這種情況下做這種卡方檢驗并不能說明什么問題抗蠢，課本是看兩者之間有沒有線性趨勢，我們可以直接用lm()函數(shù)做思劳，把a(bǔ)ge作為自變量迅矛，把level作為因變量即可，由于沒有原始數(shù)據(jù)潜叛，這里就不演示了秽褒。

多個樣本率間的多重比較

主要有卡方分割法壶硅、Scheffe可信區(qū)間法、SNK法等销斟，這里主要演示卡方分割法庐椒。

其實非常簡單，就是把多個組手動拆分為多個 兩個組蚂踊，分別進(jìn)行卡方檢驗约谈，和P值比較，只不過這里的P值不再是0.05悴势，而是和組數(shù)（比較次數(shù)）有關(guān)窗宇。

使用例7-10的數(shù)據(jù)。

df <- foreign::read.spss("./例07-06物理藥物外用膏用.sav", to.data.frame = T,reencode = "utf-8")

M <- matrix(df$f,nrow = 3,byrow = T,
            dimnames = list(trt = c("物理", "藥物", "外用"),
                            effect = c("有效","無效")))

M

Show in New Window
re-encoding from utf-8
      effect
trt    有效 無效
  物理  199    7
  藥物  164   18
  外用  118   26

手動拆分特纤，兩兩比較军俊，直接取子集即可：

# 物理治療組和藥物治療組的卡方檢驗
chisq.test(M[1:2,], correct = F)

    Pearson's Chi-squared test

data:  M[1:2, ]
X-squared = 6.756, df = 1, p-value = 0.009343

# 物理治療組和外用膏藥組的卡方檢驗
chisq.test(M[c(1,3),], correct = F)

    Pearson's Chi-squared test

data:  M[c(1, 3), ]
X-squared = 21.323, df = 1, p-value = 3.881e-06

# 藥物治療組和外用膏藥組的卡方檢驗
chisq.test(M[2:3,], correct = F)

    Pearson's Chi-squared test

data:  M[2:3, ]
X-squared = 4.591, df = 1, p-value = 0.03214

可以看到和課本是一樣的。

這時的 P' = P / (K * (K - 1) / 2 + 1)捧存，K是組數(shù)粪躬，一般情況下P=0.05，所以P' = 0.05/(3*(3-1)/2+1) = 0.0125昔穴，上面3個卡方分析的P值和0.0125比較即可镰官！

Cochran-Mantel-Haenszel 卡方統(tǒng)計量檢驗

中文名又叫行均分檢驗，常用于按照某個變量進(jìn)行分層后的檢驗吗货，這個方法課本上說用于檢驗兩個有序分類變量是否存在線性相關(guān)泳唠，但實際上用途很廣泛，比如因變量是有序變量的單向有序列聯(lián)表宙搬，也可以用笨腥。

使用課本例7-12的數(shù)據(jù)。

這個數(shù)據(jù)有3個變量勇垛，首先是年齡脖母，根據(jù)年齡分成兩層，然后是是否心肌梗死和是否口服避孕藥闲孤，我們可以直接把這個數(shù)據(jù)錄入成3維array的形式：

myo <- array(c(17,47,
               121,944,
               12,158,
               14,663),
             dim = c(2,2,2),
             dimnames = list(心肌梗死 = c("病例","對照"),
                             口服避孕藥 = c("是","否"),
                             年齡分層 = c("<40歲","≥40歲")
                             )
             )
myo
## , , 年齡分層 = <40歲
## 
##         口服避孕藥
## 心肌梗死 是  否
##     病例 17 121
##     對照 47 944
## 
## , , 年齡分層 = ≥40歲
## 
##         口服避孕藥
## 心肌梗死  是  否
##     病例  12  14
##     對照 158 663

這樣就能直接進(jìn)行Cochran-Mantel-Haenszel檢驗了谆级，這個檢驗的函數(shù)是R語言自帶的，不需要另外的包：

mantelhaen.test(myo,correct = F)
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  myo
## Mantel-Haenszel X-squared = 24.184, df = 1, p-value = 8.755e-07
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.930775 4.933840
## sample estimates:
## common odds ratio 
##          3.086444

這樣就得到結(jié)果了讼积，這個結(jié)果和課本也是一樣的肥照。

課本還介紹了Breslow-Day對各層的效應(yīng)值進(jìn)行齊性檢驗，這個檢驗可以通過DescTools包實現(xiàn)：

library(DescTools)
## Registered S3 method overwritten by 'DescTools':
##   method         from 
##   reorder.factor gdata

BreslowDayTest(myo)
## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  myo
## X-squared = 0.23409, df = 1, p-value = 0.6285

結(jié)果也是和課本一模一樣勤众。
如果你的是原始數(shù)據(jù)的形式建峭，也是很簡單的，我們構(gòu)造一個和上面數(shù)據(jù)一樣的原始數(shù)據(jù)：

myoo <- data.frame(年齡分層 = c(rep("<40歲",1129),rep("≥40歲",847)),
                   心肌梗死 = c(rep("病例",64),rep("對照",1065),
                                rep("病例",170),rep("對照",677)
                                ),
                   口服避孕藥 = c(rep("是",17),rep("否",47),
                                  rep("是",121),rep("否",944),
                                  rep("是",12),rep("否",158),
                                  rep("是",14),rep("否",663)
                                  )
                   )

# 分類變量變?yōu)橐蜃有?myoo$年齡分層 <- factor(myoo$年齡分層,levels = c("<40歲","≥40歲"))
myoo$心肌梗死 <- factor(myoo$心肌梗死, levels = c("病例","對照"))
myoo$口服避孕藥 <- factor(myoo$口服避孕藥, levels = c("是","否"))

head(myoo)
##   年齡分層 心肌梗死 口服避孕藥
## 1    <40歲     病例         是
## 2    <40歲     病例         是
## 3    <40歲     病例         是
## 4    <40歲     病例         是
## 5    <40歲     病例         是
## 6    <40歲     病例         是

用xtabs查看數(shù)據(jù)决摧，結(jié)果和我們的array的形式是一樣的：

myoo.tab <- xtabs(~口服避孕藥+心肌梗死+年齡分層,data = myoo)
myoo.tab
## , , 年齡分層 = <40歲
## 
##           心肌梗死
## 口服避孕藥 病例 對照
##         是   17  121
##         否   47  944
## 
## , , 年齡分層 = ≥40歲
## 
##           心肌梗死
## 口服避孕藥 病例 對照
##         是   12   14
##         否  158  663

這樣就可以直接進(jìn)行Cochran-Mantel-Haenszel檢驗了：

mantelhaen.test(myoo.tab, correct = F)
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  myoo.tab
## Mantel-Haenszel X-squared = 24.184, df = 1, p-value = 8.755e-07
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.930775 4.933840
## sample estimates:
## common odds ratio 
##          3.086444

結(jié)果也是一樣的。

還可用woolf法檢驗不同分層之間的效應(yīng)值有沒有統(tǒng)計學(xué)顯著性，通過使用?mantelhaen.test查看幫助文檔掌桩，作者直接給了一個現(xiàn)成的計算方法：

woolf <- function(x) {
  x <- x + 1 / 2
  k <- dim(x)[3]
  or <- apply(x, 3, function(x) (x[1,1]*x[2,2])/(x[1,2]*x[2,1]))
  w <-  apply(x, 3, function(x) 1 / sum(1 / x))
  1 - pchisq(sum(w * (log(or) - weighted.mean(log(or), w)) ^ 2), k - 1)
}

woolf(myoo.tab)
## [1] 0.6400154

直接給出了P值边锁。

頻數(shù)分布擬合優(yōu)度卡方檢驗

使用課本例7-13的數(shù)據(jù)。

R語言做卡方擬合優(yōu)度檢驗非常簡單波岛，關(guān)鍵是概率的計算茅坛，這里我們直接用課本中的概率。

x <- c(26,51,75,63,38,17,9)
p <- c(0.0854,0.2102,0.2585,0.2120,0.1304,0.0641,0.0394)

chisq.test(x=x, p =p)
## 
##  Chi-squared test for given probabilities
## 
## data:  x
## X-squared = 2.0377, df = 6, p-value = 0.9162

結(jié)果和課本非常接近则拷。

這里貼一個網(wǎng)絡(luò)教程的概率計算方法：

x<-0:6
y<-c(26,51,75,63,38,17,9)
mean<-mean(rep(x,y))
q<-ppois(x,mean)
n<-length(y)
p<-c()
p[1]<-q[1]
p[n]<-1-q[n-1]
for(i in 2:(n-1))
  p[i]<-q[i]-q[i-1]
chisq.test(y, p=p,correct=F)

    Chi-squared test for given probabilities

data:  y
X-squared = 2.0569, df = 6, p-value = 0.9144

結(jié)果和課本非常接近贡蓖。

本文首發(fā)于公眾號：醫(yī)學(xué)和生信筆記

醫(yī)學(xué)和生信筆記，專注R語言在臨床醫(yī)學(xué)中的使用煌茬，R語言數(shù)據(jù)分析和可視化斥铺。主要分享R語言做醫(yī)學(xué)統(tǒng)計學(xué)、meta分析坛善、網(wǎng)絡(luò)藥理學(xué)晾蜘、臨床預(yù)測模型、機(jī)器學(xué)習(xí)眠屎、生物信息學(xué)等剔交。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末，一起剝皮案震驚了整個濱河市改衩，隨后出現(xiàn)的幾起案子岖常，更是在濱河造成了極大的恐慌，老刑警劉巖葫督，帶你破解...
    沈念sama閱讀 206,013評論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件竭鞍，死亡現(xiàn)場離奇詭異，居然都是意外死亡候衍，警方通過查閱死者的電腦和手機(jī)笼蛛，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來蛉鹿，“玉大人滨砍，你說我怎么就攤上這事⊙欤” “怎么了惋戏？”我有些...
    開封第一講書人閱讀 152,370評論 0贊 342
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵，是天一觀的道長他膳。 經(jīng)常有香客問我响逢，道長，這世上最難降的妖魔是什么棕孙？ 我笑而不...
    開封第一講書人閱讀 55,168評論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任舔亭，我火速辦了婚禮些膨，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘钦铺。我一直安慰自己订雾，他們只是感情好，可當(dāng)我...
    茶點故事閱讀 64,153評論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布矛洞。 她就那樣靜靜地躺著洼哎，像睡著了一般。 火紅的嫁衣襯著肌膚如雪沼本。 梳的紋絲不亂的頭發(fā)上噩峦，一...
    開封第一講書人閱讀 48,954評論 1贊 283
  • 城市分裂傳說
    那天，我揣著相機(jī)與錄音抽兆，去河邊找鬼识补。 笑死，一個胖子當(dāng)著我的面吹牛郊丛，可吹牛的內(nèi)容都是我干的李请。 我是一名探鬼主播，決...
    沈念sama閱讀 38,271評論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼厉熟，長吁一口氣：“原來是場噩夢啊……” “哼导盅！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起揍瑟，我...
    開封第一講書人閱讀 36,916評論 0贊 259
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤白翻，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后绢片，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體滤馍，經(jīng)...
    沈念sama閱讀 43,382評論 1贊 300
  • ?護(hù)林員之死
    正文 獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,877評論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年底循，在試婚紗的時候發(fā)現(xiàn)自己被綠了巢株。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點故事閱讀 37,989評論 1贊 333
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡熙涤，死狀恐怖阁苞，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情祠挫，我是刑警寧澤那槽，帶...
    沈念sama閱讀 33,624評論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布，位于F島的核電站等舔，受9級特大地震影響骚灸，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜慌植，卻給世界環(huán)境...
    茶點故事閱讀 39,209評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一甚牲、第九天 我趴在偏房一處隱蔽的房頂上張望义郑。 院中可真熱鬧，春花似錦鳖藕、人聲如沸魔慷。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,199評論 0贊 19
  • 一樁弒父案著恩，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至蜻展，卻和暖如春喉誊，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背纵顾。 一陣腳步聲響...
    開封第一講書人閱讀 31,418評論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工伍茄， 沒想到剛下飛機(jī)就差點兒被人妖公主榨干…… 1. 我叫王不留，地道東北人施逾。 一個月前我還...
    沈念sama閱讀 45,401評論 2贊 352
  • 代替公主和親
    正文 我出身青樓敷矫，卻偏偏與公主長得像，于是被迫代替她去往敵國和親汉额。 傳聞我的和親對象是個殘疾皇子曹仗，可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 42,700評論 2贊 345

推薦閱讀更多精彩內(nèi)容