定義
Pearson's chi squared test：
Pearson卡方檢驗(yàn)杏死，由著名統(tǒng)計(jì)學(xué)家Karl Pearson提出泵肄，主要是比較兩個(gè)及兩個(gè)以上樣本構(gòu)成比以及兩個(gè)分類(lèi)變量的關(guān)聯(lián)性的分析。其根本思想是比較理論頻數(shù)與實(shí)際頻數(shù)的溫和程度或擬合優(yōu)度的問(wèn)題淑翼。廣泛應(yīng)用于1）分類(lèi)變量的獨(dú)立性（是否具有相關(guān)性）檢驗(yàn)中腐巢；2）是否服從某個(gè)比例的比較檢驗(yàn)中⌒ǎ可應(yīng)用于分類(lèi)變量（categorical data）的獨(dú)立性檢驗(yàn)中系忙，也可用于分類(lèi)變量的比較檢驗(yàn)中。這兩種檢驗(yàn)都需要用到R×C列聯(lián)表（R×C contingency table）惠豺，其中R表示行（Row）银还，C表示列（Column）。最簡(jiǎn)單的情形是行與列都是二分類(lèi)無(wú)序變量洁墙，這樣的數(shù)據(jù)也稱(chēng)為四格表資料蛹疯。
卡方檢驗(yàn)使用條件：1.隨機(jī)樣本數(shù)據(jù)； 2.卡方檢驗(yàn)的理論頻數(shù)不能太小热监。 兩個(gè)獨(dú)立樣本比較可以分以下3種情況： 1.所有的理論數(shù)T≥5并且總樣本量n≥40捺弦，用Pearson卡方進(jìn)行檢驗(yàn)。 2.如果理論數(shù)T＜5但T≥1孝扛，并且n≥40列吼，用連續(xù)性校正的卡方進(jìn)行檢驗(yàn)。 3.如果有理論數(shù)T＜1或n＜40苦始，則用Fisher’s檢驗(yàn)寞钥。 R×C表卡方檢驗(yàn)應(yīng)用條件： 1.R×C表中理論數(shù)小于5的格子不能超過(guò)1/5；2.不能有小于1的理論數(shù)陌选。我的實(shí)驗(yàn)中也不符合R×C表的卡方檢驗(yàn)理郑。可以通過(guò)增加樣本數(shù)咨油、列合并來(lái)實(shí)現(xiàn)您炉。
McNemar‘s test（McNemar's test for correlated proportions，也叫配對(duì)卡方檢驗(yàn)） 和Kappa 一致性檢驗(yàn)：
兩種方法均用于配對(duì)資料率的檢驗(yàn)役电，如檢驗(yàn)兩種治療方式赚爵。二者區(qū)別：1、Kappa檢驗(yàn)旨在評(píng)價(jià)兩種方法是否存在一致性；配對(duì)χ2檢驗(yàn)主要確定兩種方法診斷結(jié)果是否有差別冀膝；2唁奢、Kappa檢驗(yàn)會(huì)利用列聯(lián)表的全部數(shù)據(jù)，而配對(duì)χ2檢驗(yàn)只利用“不一致“數(shù)據(jù)畸写；3驮瞧、Kappa檢驗(yàn)可計(jì)算Kappa值用于評(píng)價(jià)一致性大小，而配對(duì)χ2檢驗(yàn)只能給出兩種方法差別是否具有統(tǒng)計(jì)學(xué)意義的判斷枯芬。Kappa值判斷標(biāo)準(zhǔn)：Kappa≥0.75论笔，說(shuō)明兩種方法診斷結(jié)果一致性較好; 0.4≤Kappa<0.75，說(shuō)明兩種方法診斷結(jié)果一致性一般; Kappa<0.4千所，說(shuō)明兩種方法診斷結(jié)果一致性較差狂魔。
使用示例
Rscript this R <輸入文件> <需要檢驗(yàn)的變量列表隆箩，來(lái)自于輸入文件的列名> <輸出路徑> <輸出文件前綴> <基線/金標(biāo)準(zhǔn)措译，來(lái)自于輸入文件的列名东亦，“需要檢驗(yàn)的變量列表”均與該變量進(jìn)行檢驗(yàn)>
Rscript chisq.test.r data_merge.txt list /lustre/work/user/zhou.yang/project/10_ZL_multi/2_analysis/6_single_sample/10_shannon_diff_mutation_set/6.SVV/group_by_diff_method_compare/statist_test_result Common_mutation Common_mutation
infile文件示例
list文件示例
代碼：
args = commandArgs(T)
if (length(args) !=5){
print("Rscript this R <infile> <list> <outdir> <prefix> <baseline>")
q()
}
library(data.table)
library("irr")
data<-read.table(args[1],sep="\t",header=T,encoding='UTF-8',check.names=F,row.names=1)
list<-read.table(args[2],sep="\t",header=F,encoding='UTF-8',check.names=F)
len=dim(list)[1]
print(dim(data))
header<-paste("Item","X-squared","df","p-value",sep="\t")
write.table(header,file=(paste(args[3],"/",args[4],".chi-square_test_result.xls",sep="")),sep="\t",quote=F,row.names=F,col.names=F,append=F)
header<-paste("Item","X-squared","df","p-value",sep="\t")
write.table(header,file=(paste(args[3],"/",args[4],".mcnemar-square_test_result.xls",sep="")),sep="\t",quote=F,row.names=F,col.names=F,append=F)
header<-paste("Item","method","Kappa","p-value",sep="\t")
write.table(header,file=(paste(args[3],"/",args[4],".kappa-square_test_result.xls",sep="")),sep="\t",quote=F,row.names=F,col.names=F,append=F)
for(i in 1:len){
id=as.character(list[i,])
data1<-as.data.frame(t(data))
traits<-as.data.frame(t(data1[which(rownames(data1)==id),]))
traits$Group=data[,which(colnames(data)==args[5])]
colnames(traits)<-c("Group1","Group2")
print(dim(traits))
#traits$Group1=as.factor(traits$Group1)
#traits$Group2=as.factor(traits$Group2)
traits<-na.omit(traits)
traits <- subset(traits, traits$Group2 != "")
#生成4聯(lián)表
four_table<-xtabs(~traits$Group1+traits$Group2,data=traits)
#獨(dú)立性檢驗(yàn)/差異檢驗(yàn)：卡方檢驗(yàn)
#理論頻數(shù)T＜5（chisq.test(four_table)$expected）時(shí)转砖，進(jìn)行連續(xù)性校正的卡方進(jìn)行檢驗(yàn)：correct=TURE。
k<-"FALSE"
judge<- chisq.test(four_table)$expected >5
if (k %in% judge){chisq_v<-chisq.test(four_table,correct=T)}
else{chisq_v<-chisq.test(four_table,correct=F)}
# print(chisq.test(four_table)$expected)
# print(chisq_v)
chisq_va<-paste(id,chisq_v$statistic[[1]],chisq_v$parameter[[1]],chisq_v$p.value,sep="\t")
write.table(chisq_va,file=(paste(args[3],"/",args[4],".chi-square_test_result.xls",sep="")),sep="\t",quote=F,col.names=F,row.names=F,append=T)
#一致性檢：配對(duì)卡方檢驗(yàn)（McNemar 檢驗(yàn)）
mcnemar_v<-mcnemar.test(four_table,correct=T)
mcnemar_va<-paste(id,mcnemar_v$statistic[[1]],mcnemar_v$parameter[[1]],mcnemar_v$p.value,sep="\t")
write.table(mcnemar_va,file=(paste(args[3],"/",args[4],".mcnemar-square_test_result.xls",sep="")),sep="\t",quote=F,col.names=F,row.names=F,append=T)
#一致性檢：kappa檢驗(yàn)：對(duì)于僅探討陽(yáng)性率是否相同套菜，可采用配對(duì)卡方檢驗(yàn)（又稱(chēng)McNemar檢驗(yàn)）脚牍；如果探討結(jié)果是否具有一致性伐蒋，則采用Kappa檢驗(yàn)火俄。通常犯建，在診斷、檢驗(yàn)領(lǐng)域瓜客，Kappa檢驗(yàn)更為多見(jiàn)适瓦。
#Kappa≥0.75，說(shuō)明兩種方法診斷結(jié)果一致性較好;0.4≤Kappa<0.75谱仪，說(shuō)明兩種方法診斷結(jié)果一致性一般;Kappa<0.4玻熙，說(shuō)明兩種方法診斷結(jié)果一致性較差。
kappa_v<-kappa2(traits,'unweighted')
#Weight參數(shù)是函數(shù)的重點(diǎn)疯攒，如果有多個(gè)檢查項(xiàng)目中有一個(gè)是為0的時(shí)候需要加權(quán)檢驗(yàn)嗦随，其他時(shí)候一般都是非加權(quán)檢驗(yàn)。
kappa_va<-paste(id,kappa_v$method,kappa_v$value,kappa_v$p.value,sep="\t")
write.table(kappa_va,file=(paste(args[3],"/",args[4],".kappa-square_test_result.xls",sep="")),sep="\t",quote=F,col.names=F,row.names=F,append=T)
}
輸出文件
一致性結(jié)果可視化（基于Kappa檢驗(yàn)結(jié)果）
結(jié)果示例：
腳本
args = commandArgs(T)
if (length(args) !=4){
print("Rscript this R <infile> <X> <Y> <outfile>")
q()
}
library("ggplot2")
dat<-read.table(file=args[1],header=T,sep="\t")
data<-na.omit(dat)
data$judge <- ifelse(data$kappa.value >=0.75,"Good",ifelse(data$kappa.value >=0.4,"Moderate","Bad"))
data$judge <- as.factor(data$judge)
data$kappa <- data$kappa.value
data1<-data[data$P.value < 0.05,]
pdf(file=args[4],height=4,width=10)
ggplot(data=data) + geom_point(mapping = aes(x = factor(get(args[2]),levels=as.factor(unique(get(args[2])))), y = factor(get(args[3]),levels=as.factor(unique(get(args[3])))), size=kappa,color=judge))+theme_bw()+theme(panel.grid =element_blank())+ theme(axis.ticks= element_blank()) + scale_color_manual("Consistency (Cohen's Kappa)",values=c(Good = "red2", Moderate = "orange",Bed="Blue"))+scale_size_continuous("Kappa consistency Value")+ geom_point(data= data1,alpha=0.9,pch=21,colour="black",aes(x = factor(get(args[2]),levels=as.factor(unique(get(args[2])))), y = factor(get(args[3]),levels=as.factor(unique(get(args[3])))),size=kappa))+xlab("")+ylab("")+guides(size = guide_legend(order = 1))+theme(axis.text.x = element_text(size=12))+ theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black", size = 9),axis.title.y = element_text(size=10))+ guides(color = guide_legend(override.aes = list(size = 4)))
#ggplot(data=data) + geom_point(mapping = aes(x = factor(get(args[2]),levels=as.factor(unique(get(args[2])))), y = factor(get(args[3]),levels=as.factor(unique(get(args[3])))), size=kappa,color=judge))+theme_bw()+theme(panel.grid =element_blank())+ theme(axis.ticks= element_blank())+theme(panel.border = element_blank()) + scale_color_manual("Consistency",values=c(Good = "red2", Moderate = "orange",Bed="Blue"))+scale_size_continuous("Kappa consistency Value")+ geom_point(data= data1,alpha=0.9,pch=21,colour="black",aes(x = factor(get(args[2]),levels=as.factor(unique(get(args[2])))), y = factor(get(args[3]),levels=as.factor(unique(get(args[3])))),size=kappa))+xlab("")+ylab("")+guides(size = guide_legend(order = 1))+theme(axis.text.x = element_text(size=12))+ theme(axis.text.x = element_text(angle = 45, hjust = 1, color = "black", size = 9),axis.title.y = element_text(size=10))+ guides(color = guide_legend(override.aes = list(size = 4)))
dev.off()
使用方法：
Rscript this R <infile> <X> <Y> <outfile>
Rscript this R <輸入文件> <X軸> <Y軸> <輸出文件>
使用示例：
Rscript kappa_Bubble_plot.r kappa.format shannon Group kappa_Bubble_plot.pdf
輸入文件示例（統(tǒng)計(jì)分析的結(jié)果需要整理一下）：
