統(tǒng)計(jì)學(xué)第七章 卡方檢驗(yàn)【R語言實(shí)現(xiàn)】

知識(shí)點(diǎn)

  • 卡方分布(Chi-square distribution)
  • 四格表資料的卡方檢驗(yàn)
  • 配對(duì)四格表的卡方檢驗(yàn)
  • 四格表資料的Fisher確切概率法
  • 行*列資料的卡方檢驗(yàn)
  • 多個(gè)樣本率間的多重比較
  • 頻數(shù)分布擬合優(yōu)度的卡方檢驗(yàn)

1. 卡方分布

當(dāng)自由度趨于無窮大時(shí)，圖形趨于正態(tài)分布

# 產(chǎn)生服從卡方分布的觀測(cè)數(shù)為1000的樣本
df_n <- seq(1, 5)
chisq_data <- function(n){
  x <- seq(-0.2, 16, length=300)
  prob <- dchisq(x, df=n[1, 1])
  return(data.frame(x=x, prob=prob, df=n[1,1]))
}
require(plyr, quietly=TRUE)
require(ggplot2, quietly=TRUE)
data <- ddply(data.frame(n=df_n), .(n), chisq_data)
head(data)

##   n           x     prob df
## 1 1 -0.20000000 0.000000  1
## 2 1 -0.14581940 0.000000  1
## 3 1 -0.09163880 0.000000  1
## 4 1 -0.03745819 0.000000  1
## 5 1  0.01672241 3.059352  1
## 6 1  0.07090301 1.446043  1

tail(data)

##      n        x        prob df
## 1495 5 15.72910 0.003186504  5
## 1496 5 15.78328 0.003117378  5
## 1497 5 15.83746 0.003049697  5
## 1498 5 15.89164 0.002983433  5
## 1499 5 15.94582 0.002918558  5
## 1500 5 16.00000 0.002855045  5

ggplot(data, aes(x=x, y=prob, color=factor(df), group=df))+geom_line(lwd=1)+scale_y_continuous(limits=c(0, 0.7))
  1. 其定義為：若n個(gè)相互獨(dú)立的隨機(jī)變量ξ?、ξ?、……、ξn ，均服從標(biāo)準(zhǔn)正態(tài)分布（也稱獨(dú)立同分布于 標(biāo)準(zhǔn)正態(tài)分布），則這n個(gè)服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量的平方和服從卡方分布
  2. 可加性：兩個(gè)服從卡方分布的獨(dú)立隨機(jī)變量相加服從自由度為兩自由度之和的卡方分布
  3. 卡方檢驗(yàn)的基本思想：由于在假設(shè)符合某種情況的前提下，樣本實(shí)際值偏離理論值的偏差服從正態(tài)分布，其均值為理論值，方差也為理論值？？？（有點(diǎn)疑惑）

2. 四格表資料的卡方檢驗(yàn)

書上P98例7-2：表格為

組別 有效 無效 合計(jì)
胞磷膽堿組 46 6 52
神經(jīng)節(jié)苷酯組 18 8 26
合計(jì) 64 14 78

H0：兩種藥物療效相同
H1：有效率不等

table7_2 <- matrix(c(46, 18, 6, 8), nrow=2, ncol=2)
chisq.test(table7_2)

## Warning in chisq.test(table7_2): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table7_2
## X-squared = 3.1448, df = 1, p-value = 0.07617

得到warning "Chi-squared approximation may be incorrect"
因?yàn)楸砀裰杏蠺<5, 此時(shí)可以采用校正【自動(dòng)校正】或者fisher.test()
可以用以下代碼查看理論值

chisq.test(table7_2)$expected

## Warning in chisq.test(table7_2): Chi-squared approximation may be incorrect

##          [,1]     [,2]
## [1,] 42.66667 9.333333
## [2,] 21.33333 4.666667

參考：
http://r.789695.n4.nabble.com/In-chisq-test-x-Chi-squared-approximation-may-be-incorrect-td845040.html

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

  • 超幾何分布
    從一個(gè)有限總體中進(jìn)行不放回抽樣，設(shè)N件產(chǎn)品，有M件不合格品，若從中不放回地隨機(jī)抽取n件，則其中含有的不合格品件數(shù)X服從超幾何分布，記為X~h(n, N, M)
    P100 例7-4
    組別 感染 未感染 合計(jì)
    預(yù)防組 4 18 22
    非預(yù)防組 5 6 11
    合計(jì) 9 24 33

假設(shè)兩組（預(yù)防組和非預(yù)防組）的感染率都是9：33【零假設(shè)】，則邊緣值固定的情況下，相當(dāng)于在總數(shù)33的所有個(gè)體中【有9個(gè)感染的，24個(gè)未感染的】，取22個(gè)值作為有效組，在這22個(gè)值中，記感染的人數(shù)為X，則X~h(22, 9, 24)。
H0: 兩個(gè)組無查別

# x代表取出來的白球數(shù)， k代表取的次數(shù)，m代表總白球數(shù)，n代表總黑球數(shù)
sum(dhyper(x=0:9, k=22, m=9, n=24))

## [1] 1

p_values <- dhyper(x=0:9, k=22, m=9, n=24)
p <- sum(p_values[p_values<=dhyper(x=4, k=22, m=9, n=24)])
print(p)

## [1] 0.1210448

fisher.test(matrix(c(4, 5, 18, 6), nrow=2, ncol=2))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(4, 5, 18, 6), nrow = 2, ncol = 2)
## 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

4. 配對(duì)四格表的卡方檢驗(yàn)

  • 在次處的配對(duì)即同一樣本接受兩種處理
  • 用結(jié)果不一致的兩種情況計(jì)算統(tǒng)計(jì)量卡方（分b+c>=40和b+c<40的情況）
  • 稱為McNemar卡方檢驗(yàn)
mat <- matrix(c(11, 2, 12, 33), nrow=2)
# 方法1：
if (mat[1,2]+mat[2,1]>=40) {
  x_sq <- (mat[1,2]-mat[2,1])^2/(mat[1,2]+mat[2,1])
}else { x_sq <- (abs(mat[1,2]-mat[2,1])-1)^2/(mat[1,2]+mat[2,1]) }
cat("chi-squareed = ", x_sq, 
    "    p-value = ", 
    pchisq(x_sq, df=1, lower.tail=FALSE), 
    sep="")

## chi-squareed = 5.785714    p-value = 0.01615693

# 方法2：
mcnemar.test(mat)

## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  mat
## McNemar's chi-squared = 5.7857, df = 1, p-value = 0.01616

參考：
https://stat.ethz.ch/R-manual/R-devel/library/stats/html/mcnemar.test.html

5. 行*列資料的卡方檢驗(yàn)

  • 多個(gè)樣本率的比較（同四格表代碼）

  • 樣本構(gòu)成比的比較（同四格表代碼）

  • 雙向無序分類變量的關(guān)聯(lián)性檢驗(yàn)
    方法1：

mat <- matrix(
  c(431, 388, 495, 137, 490, 410, 587, 179, 902, 800, 950, 32),
  nrow=4)

chisq.test(mat)

## 
##  Pearson's Chi-squared test
## 
## data:  mat
## X-squared = 213.16, df = 6, p-value < 2.2e-16

x_sq <- chisq.test(mat)$statistic[[1]]
ContCoef <- sqrt(x_sq/(x_sq+sum(mat)))
print(ContCoef)

## [1] 0.1882638

方法2：

# install.packages("DescTools")
mat <- matrix(
  c(431, 388, 495, 137, 490, 410, 587, 179, 902, 800, 950, 32),
  nrow=4)
DescTools::ContCoef(mat)

## [1] 0.1882638

參考：
https://www.rdocumentation.org/packages/DescTools/versions/0.99.19/topics/Association%20measures

  • 雙向有序分組變量的線性趨勢(shì)檢驗(yàn)
    • 線性趨勢(shì)：對(duì)卡方進(jìn)行分解，分解為回歸分量卡方和偏離線性回歸分量卡方
    • 相關(guān)關(guān)系：等級(jí)相關(guān)分析
    • 差別分析：視為單向有序的秩轉(zhuǎn)換非參數(shù)檢驗(yàn)（如不同年齡組的療效差別）
  • 雙向有序?qū)傩韵嗤?
    • 相當(dāng)為配伍資料（不同檢測(cè)方法，同一樣本）
    • 一致性檢驗(yàn)（Kappa檢驗(yàn)）
    • 特殊模型分析
  • 單向有序
    • 分組為有序，指標(biāo)為無序（卡方檢驗(yàn)）
    • 分組為無序，指標(biāo)為有序（秩轉(zhuǎn)換非參數(shù)檢驗(yàn)）

6. 多個(gè)樣本率之間的多重比較

  • 卡方分割法
  • 多個(gè)實(shí)驗(yàn)組兩兩比較（alpha=alpha/(choose(2, k)+1)）
  • 實(shí)驗(yàn)組與同一個(gè)對(duì)照組比較（alpha=alpha/(2*(k-1))）

7. 頻數(shù)分布擬合度的卡方檢驗(yàn)

  • 卡方=sum((理論頻數(shù)-觀察頻數(shù))^2/(理論頻數(shù)))
  • 自由度=行數(shù)-(計(jì)算理論頻數(shù)時(shí)使用的統(tǒng)計(jì)量個(gè)數(shù)如總例數(shù)或均數(shù)等)
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

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

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