使用R语言进行卡方检验(chi-square test)

Python013

使用R语言进行卡方检验(chi-square test),第1张

文|程瑞林(山东大学第二医院足踝外科)

来源|(微信公众号)云中瑞麟(ID:ruilinfly)

瑞麟导读:

对于计量资料,临床医学研究中常用的统计分析方法是t检验;而对于计数资料,卡方检验是一个常用的统计分析方法。

最近看到一篇文章,里面分析了骨巨细胞瘤患者术后复发的比例,其中计数资料使用卡方检验(又称χ 2 检验),下面针对卡方检验的使用方法及其R语言实现方法进行简单介绍。

卡方检验是一种用途很广的 计数资料 的假设检验方法,由卡尔·皮尔逊提出。它属于非参数检验的范畴,主要是比较两个及两个以上样本率( 构成比)以及两个分类变量的关联性分析。其根本思想就是在于比较理论频数和实际频数的吻合程度或拟合优度问题。

它在分类资料统计推断中的应用,包括:两个率或两个构成比比较的卡方检验;多个率或多个构成比比较的卡方检验以及分类资料的相关分析等。

可以分为成组比较(不配对资料)和个别比较(配对,或同一对象两种处理的比较)两类。

通常卡方检验的应用主要为:

1、 卡方拟合优度检验

2、卡方独立性检验

我们想知道喝牛奶对感冒发病率有没有影响,以下为数据统计的四格表:

通过简单的统计我们得出喝牛奶组和不喝牛奶组的感冒率为30.94%和25.00%,两者的差别可能是抽样误差导致,也有可能是牛奶对感冒率真的有影响。

为了确定真实原因,我们先假设喝牛奶对感冒发病率是没有影响的,即喝牛奶喝感冒时独立无关的,所以我们可以得出感冒的发病率实际是(43+28)/(43+28+96+84)= 28.29%

所以,理论的四格表应该如下表所示:

即下表:

如果喝牛奶喝感冒真的是独立无关的,那么四格表里的理论值和实际值差别应该会很小。

那如何来描述这种差别呢,我们定义卡方值为

其中,A为实际值,T为理论值。

x2用于衡量实际值与理论值的差异程度(也就是卡方检验的核心思想),包含了以下两个信息:

根据卡方检验公式我们可以得出例1的卡方值为:

卡方 = (43 - 39.3231)平方 / 39.3231 + (28 - 31.6848)平方 / 31.6848 + (96 - 99.6769)平方 / 99.6769 + (84 - 80.3152)平方 / 80.3152 = 1.077

卡方值(理论值与实际值差异大小)的意义是什么呢?为此我们再引入一个概念:

上一步我们得到了卡方的值,但是如何通过卡方的值来判断喝牛奶和感冒是否真的是独立无关的?也就是说,怎么知道无关性假设是否可靠?

答案是,通过查询卡方分布的临界值表。

第一行表示显著性水平α

第一列表示自由度

这里需要用到一个 自由度 的概念,自由度等于V = (行数 - 1) * (列数 - 1),对四格表,自由度V = 1。

对V = 1,喝牛奶和感冒(95%概率)不相关的卡方分布的临界值(最大)是:3.84。即如果卡方大于3.84,则认为喝牛奶和感冒(有95%的概率)相关。

【瑞麟描述】临界值3.84的意义表示:如果卡方值>3.84,则纵列因素与横行因素不相关的的概念<0.05(即显著性水平),也即纵列因素与横行因素相关的概念>0.95。

显然1.077<3.84,没有达到卡方分布的临界值,所以喝牛奶和感冒独立不相关的假设没有被推翻。

【瑞麟】 简单说,如果我们计算出的卡方值(表示实际值与理论值的差异,越大表示实际值与理论值越不符,即越有可能纵列因素会影响横行数值)大于临界值(列因素不影响横行值的范围:0~临界值),我们就排斥原假设(H0,即纵列因素不影响横行的因素的变化),接受备择假设(H1:纵列因素对横行的因素变化有影响);反之,卡方值小于临界值,即在(纵列与横行互不影响这一假设)理论范围内,无法推翻原假设,即无统计差异。

我们想知道不吃晚饭对体重下降有没有影响,并获得以下数据:

H0:r1=r2,不吃晚饭对体重下降没有影响,即吃不吃晚饭的体重下降率相等;

H1:r1≠r2,不吃晚饭对体重下降有显著影响,即吃不吃晚饭的体重下降率不相等。α=0.05

【瑞麟:H0为纵列因素对横行因素无影响;H1为有影响】

3.计算卡方值

根据图1所示公式,计算出卡方值为5.498

在查表之前应知本题自由度。按卡方检验的自由度v=(行数-1)×(列数-1),则该题的自由度v=(2-1)(2-1)=1,查卡方界值表,找到3.84,而本题卡方=5.498即卡方>3.84,P<0.05,差异有显著统计学意义,按显著性水平α=0.05水准,拒绝H0,可以认为两组的体重下降率有明显差别。

通过实例计算,对卡方的基本公式有如下理解:若各理论数与相应实际数相差越小,卡方值越小;如两者相同,则卡方值必为零。

x2值表是数理统计根据正态分布的定义计算出来的。 是一种近似,在自由度大于1、理论数皆大于5时,这种近似很好;当自由度为1时,尤其当1<T<5,而n>40时,应用以下校正公式:

如果观察资料的T<1或n<40时,四格表资料用上述校正法也不行,可参考预防医学专业用的医学统计学教材中的精确检验法【瑞麟:Fisher检验?】直接计算概率以作判断。

1.一般认为行×列表中不宜有1/5以上格子的理论数小于5,或有小于1的理论数。当理论数太小可采取下列方法处理:①增加样本含量以增大理论数;②删去上述理论数太小的行和列;③将太小理论数所在行或列与性质相近的邻行邻列中的实际数合并,使重新计算的理论数增大。由于后两法可能会损失信息,损害样本的随机性,不同的合并方式有可能影响推断结论,故不宜作常规方法。另外,不能把不同性质的实际数合并,如研究血型时,不能把不同的血型资料合并。

2.如检验结果拒绝检验假设,只能认为各总体率或总体构成比之间总的来说有差别,但不能说明它们彼此之间都有差别,或某两者间有差别。

R语言自带卡方检测的方法,只要调用方法chisq.test(),会自行输出X-squared卡方值, df自由度, p-value概率。

判断5种品牌啤酒的爱好者有无显著差异:

P值越大,支持原假设的证据就越强,给定显著性水平α(取0.05), 当P值小于α时,就拒绝原假设。

H0:两种药物疗效相同

H1:有效率不等

为何会提示算法可能不准确呢?计算理论值:

文献1中的数据列表为

文章提及计数资料使用χ 2 检验,而数据列表中多处数据小于5,显然应该视理论值大小选择连续性修正的卡方检验或Fisher检验更合适一些。

参考文献:

1.同志超,等。四肢骨巨细胞瘤的外科治疗分析。中华解剖与临床杂志,2018,23(3)

2.snowdroptulip, 统计学——卡方检验和卡方分布 , CDSN博客,2017

3.lijinxiu123, 卡方检验及R语言实现 ,CDSN博客,2017-3-27

4.Knowlege_上下求索, 卡方检验x2检验(chi-square test) ,CSDN博客,2016-7-7

5.x2yline, 统计学第七章 卡方检验【R语言实现】 ,,2017.10.11

6.嘉儿jy 《卡方检验中非连续性校正与连续性校正的区别!》 百度知道,2016-1-19

7.薛毅、陈立萍 编著《统计建模与R软件》,清华大学出版社,2006

8.qazonly123 《求助,下面几种状况SPSS交叉表分别该使用哪一种卡方分析,是Pearson卡方,还是Fisher‘s,还是连续性校正》 ,百度知道,2016-5-11

201808282046更新

R语言常用函数整理本篇是基础篇,即R语言自带的函数。 vector:向量 numeric:数值型向量 logical:逻辑型向量 character;字符型向量 list:列表 data.frame:数据框 c:连接为向量或列表 length:求长度 subset:求子集 seq,from:to,sequence:等差序列 rep:重复 NA:缺失值 NULL:空对象 sort,order,unique,rev:排序 unlist:展平列表 attr,attributes:对象属性 mode,class,typeof:对象存储模式与类型 names:对象的名字属性 字符型向量 nchar:字符数 substr:取子串 format,formatC:把对象用格式转换为字符串 paste()、paste0()不仅可以连接多个字符串,还可以将对象自动转换为字符串再相连,另外还能处理向量。 strsplit:连接或拆分 charmatch,pmatch:字符串匹配 grep,sub,gsub:模式匹配与替换 complex,Re,Im,Mod,Arg,Conj:复数函数 factor:因子 codes:因子的编码 levels:因子的各水平的名字 nlevels:因子的水平个数 cut:把数值型对象分区间转换为因子 table:交叉频数表 split:按因子分组 aggregate:计算各数据子集的概括统计量 tapply:对“不规则”数组应用函数 dev.new() 新建画板 plot()绘制点线图,条形图,散点图. barplot( ) 绘制条形图 dotchart( ) 绘制点图 pie( )绘制饼图. pair( )绘制散点图阵 boxplot( )绘制箱线图 hist( )绘制直方图 scatterplot3D( )绘制3D散点图. par()可以添加很多参数来修改图形 title( ) 添加标题 axis( ) 调整刻度 rug( ) 添加轴密度 grid( ) 添加网格线 abline( ) 添加直线 lines( ) 添加曲线 text( ) 添加标签 legend() 添加图例 +, -, *, /, ^, %%, %/%:四则运算 ceiling,floor,round,signif 1、round() #四舍五入 例:x <- c(3.1416, 15.377, 269.7) round(x, 0) #保留整数位 round(x, 2) #保留两位小数 round(x, -1) #保留到十位 2、signif() #取有效数字(跟学过的有效数字不是一个意思) 例:略 3、trunc() #取整 floor() #向下取整 ceiling() #向上取整 例:xx <- c(3.60, 12.47, -3.60, -12.47) trunc(xx) floor(xx) ceiling(xx) max,min,pmax,pmin:最大最小值 range:最大值和最小值 sum,prod:向量元素和,积 cumsum,cumprod,cummax,cummin:累加、累乘 sort:排序 approx和approx fun:插值 diff:差分 sign:符号函数 abs,sqrt:绝对值,平方根 log, exp, log10, log2:对数与指数函数 sin,cos,tan,asin,acos,atan,atan2:三角函数 sinh,cosh,tanh,asinh,acosh,atanh:双曲函数 beta,lbeta,gamma,lgamma,digamma,trigamma,tetragamma,pentagamma,choose ,lchoose:与贝塔函数、伽玛函数、组合数有关的特殊函数 fft,mvfft,convolve:富利叶变换及卷积 polyroot:多项式求根 poly:正交多项式 spline,splinefun:样条差值 besselI,besselK,besselJ,besselY,gammaCody:Bessel函数 deriv:简单表达式的符号微分或算法微分 array:建立数组 matrix:生成矩阵 data.matrix:把数据框转换为数值型矩阵 lower.tri:矩阵的下三角部分 mat.or.vec:生成矩阵或向量 t:矩阵转置 cbind:把列合并为矩阵 rbind:把行合并为矩阵 diag:矩阵对角元素向量或生成对角矩阵 aperm:数组转置 nrow, ncol:计算数组的行数和列数 dim:对象的维向量 dimnames:对象的维名 rownames,colnames:行名或列名 %*%:矩阵乘法 crossprod:矩阵交叉乘积(内积) outer:数组外积 kronecker:数组的Kronecker积 apply:对数组的某些维应用函数 tapply:对“不规则”数组应用函数 sweep:计算数组的概括统计量 aggregate:计算数据子集的概括统计量 scale:矩阵标准化 matplot:对矩阵各列绘图 cor:相关阵或协差阵 Contrast:对照矩阵 row:矩阵的行下标集 col:求列下标集 solve:解线性方程组或求逆 eigen:矩阵的特征值分解 svd:矩阵的奇异值分解 backsolve:解上三角或下三角方程组 chol:Choleski分解 qr:矩阵的QR分解 chol2inv:由Choleski分解求逆 ><,>,<=,>=,==,!=:比较运算符 !,&,&&,|,||,xor(): 逻辑运算符 logical: 生成逻辑向量 all, any:逻辑向量都为真或存在真 ifelse():二者择一 match, %in%:查找 unique:找出互不相同的元素 which:找到真值下标集合 duplicated:找到重复元素 optimize,uniroot,polyroot:一维优化与求根 if,else, ifelse, switch: 分支 for,while,repeat,break,next: 循环 apply,lapply,sapply,tapply,sweep:替代循环的函数。 function:函数定义 source:调用文件 ’ call:函数调用 . C,.Fortran:调用C或者Fortran子程序的动态链接库。 Recall:递归调用 browser,debug,trace,traceback:程序调试 options:指定系统参数 missing:判断虚参是否有对应实参 nargs:参数个数 stop:终止函数执行 on.exit:指定退出时执行 eval,expression:表达式计算 system.time:表达式计算计时 invisible:使变量不显示 menu:选择菜单(字符列表菜单) 其它与函数有关的还有: delay, delete.response, deparse, do.call, dput, environment , formals, format.info, interactive, is.finite, is.function, is.language, is.recursive , match.arg, match.call, match.fun, model.extract, name, parse 函数能将字符串转换为表达式expression deparse 将表达式expression转换为字符串 eval 函数能对表达式求解 substitute, sys.parent , warning, machine cat,print:显示对象 sink:输出转向到指定文件 dump,save,dput,write:输出对象 scan,read.table,readlines, load,dget:读入 ls,objects:显示对象列表 rm, remove:删除对象 q,quit:退出系统 .First,.Last:初始运行函数与退出运行函数。 options:系统选项 ?,help,help.start,apropos:帮助功能 data:列出数据集 head()查看数据的头几行 tail()查看数据的最后几行 每一种分布有四个函数: d―density(密度函数),p―分布函数,q―分位数函数,r―随机数函数。 比如,正态分布的这四个函数为dnorm,pnorm,qnorm,rnorm。下面我们列出各分布后缀,前面加前缀d、p、q或r就构成函数名: norm:正态, t:t分布, f:F分布, chisq:卡方(包括非中心) unif:均匀, exp:指数, weibull:威布尔, gamma:伽玛, beta:贝塔 lnorm:对数正态, logis:逻辑分布, cauchy:柯西, binom:二项分布, geom:几何分布, hyper:超几何, nbinom:负二项, pois:泊松 signrank:符号秩, wilcox:秩和, tukey:学生化极差 sum, mean, var, sd, min, max, range, median, IQR(四分位间距)等为统计量, sort,order,rank与排序有关, 其它还有ave,fivenum,mad,quantile,stem等。 R中已实现的有chisq.test,prop.test,t.test。 cor,cov.wt,var:协方差阵及相关阵计算 biplot,biplot.princomp:多元数据biplot图 cancor:典则相关 princomp:主成分分析 hclust:谱系聚类 kmeans:k-均值聚类 cmdscale:经典多维标度 其它有dist,mahalanobis,cov.rob。 ts:时间序列对象 diff:计算差分 time:时间序列的采样时间 window:时间窗 lm,glm,aov:线性模型、广义线性模型、方差分析 quo()等价于quote() enquo()等价于substitute()