找回密码
 立即注册
查看: 4958|回复: 2

生物统计专题 | 用R做方差分析

[复制链接]
发表于 2015-2-10 11:33:09 | 显示全部楼层 |阅读模式
本帖最后由 R语言微信号 于 2015-2-10 11:34 编辑

生物统计专题 | 用R做方差分析2015-02-09bgstnR语言中文网

扫描关注微信号,查看更多内容。

上期内容回顾
上期主要介绍了对两独立样本、多个独立样本以及随机区组设计的数据的非参数检验,下面一起来回顾一下:
两独立样本的秩和检验:检验两个样本之间是否存在差异,直接调用wilcox()函数,对x,y赋值即可进行;
多个独立样本比较的Kruskal-Wallis H检验:检验多个样本之间是否存在差异,作用相当于参数检验中的单因素方差分析,直接调用kruskal.test()函数即可
随机区组设计的多样本比较Friedman M检验:直接调用friedman.test()函数,通过用friedman.test(y~A|B,data)的格式进行调用进行相应的检验。
温馨提示:如果您对此感到有些陌生,那就赶快翻开公众号历史消息温习一下,点击文章末尾“阅读原文”可见~温故而知新哦~
本期内容导读
前几期的内容涉及到的是两样本均值及中位数的检验,然而当涉及多于两样本均值检验时,就需要用到方差分析ANOVA(analysis of variance),比如研究研究不同种类的肥料对作物产量的影响时。本期将会对该方法进行详细介绍。
方差分析的思想及原理
一句话概括:方差分析(ANOVA)通过分析数据的误差来源来检验多组数据均值是否相等,其实质是研究分类型自变量对数值型因变量的影响。
当然在使用方差分析之前,需要确定数据是否符合方差分析的应用条件:1. 各样本相互独立且均来自正态总体;2. 相互比较的各样本的总体方差相等,即具有方差齐性。
调用函数
aov(formula, data)
参数解释
formula是一个公式,比如var~var1+var2+…,var表示因变量,var1,var2,…表示自变量,data是数据在R中的名称。
单因素方差分析示例
某职业病防治所对30名矿工分别测定各期(Phase0,Phase0_1,Phase1)血清铜蓝蛋白含量(μmol/l),资料如下

Q:各期血清铜蓝蛋白含量测定结果有无差别?
A: 详细的操作步骤如下
⊙ 第1步:重铸数据
实际操作中读者会发现数据格式与函数要求格式不同,因此要想直接调用函数需要对数据进行改动,为了达到这个目的,就需要用到reshape包中的melt()函数(需要读者输入install.packages(“reshape”)来安装)。

最终mydata2的数据格式便是调用函数所需的格式了。
⊙ 第2步:变量重命名
对于某些处女座的读者会认为数据框中的变量名称不符合自己的要求,那么可以使用reshape包中的rename()函数进行修改。

当然也可以通过fix()函数在数据框中修改(这里不再演示)。
⊙ 第3步:方差分析

但这里并没有说明组内和组间差异是否显著,为了得到其差异,还需要对结果使用summary()函数才能得到其F检验值。
图中可以看到,P=0.002<0.01,说明组内和组间差异显著,说明了各期铜蓝蛋白测试结果有显著性差异。
⊙ 第4步:TukeyHSD检验进行多重比较
方差分析只能告诉我们多组数据整体之间是否存在差异,但没告诉具体哪些组之间存在差异。本例中我们知道了各期铜蓝蛋白测试结果有显著性差异,具体差异仍然不明,此时就需要用到TukeyHSD()函数,它提供了对各组均值差异的成对检验。

这个结果也可以看到,1期与0期有显著差异,但1期与0到1期差异不显著,0期与0到1期差异也不显著,单从这个测试上看,说明了铜蓝蛋白测试对于0期与1期的测试有显著性差异。

拓展延伸:
⊙ multcomp包中的cholesterol数据集数据集包含了50个患者接受降低胆固醇药物五种疗法治疗后的数据,请确定哪种疗法效果最好。
⊙ 在更改数据格式时,我们用到了reshape包中的melt()函数,但我们也可以通过循环和一些基本函数来达到相同的效果,有兴趣的读者可以尝试一下。
注:R中很多函数对数据格式的要求是类似本例mydata2的格式,此时需要对不符合格式的数据进行变换,建议读者自己编写一个较为通用的处理函数保存在R中,以备不时之需。
关注我们
—官方网站—
R语言中文网  www.r-china.net
—官方QQ群—
R语言中文论坛-2(1000人群):427060123
R语言中文论坛(2000人群,已满):74076289
Biostatistician(500):186701945
—官方微博—
新浪微博:@R语言中文网官网
—官方微信—
微信名:R语言中文网  微信号:rchinanet


点击“阅读原文”温习上期多样本非参数检验内容!干货多多,喜欢就点个鼓励一下呗~

[url=]阅读原文[/url]




微信扫一扫
关注该公众号





回复

使用道具 举报

发表于 2015-2-13 10:51:21 | 显示全部楼层
C:\Users\Administrator\Desktop\1.png
C:\Users\Administrator\Desktop\2.png

课后练习答案,不知道是否操作正确
回复

使用道具 举报

发表于 2015-2-13 10:52:30 | 显示全部楼层
> myresult<-aov(response~trt,cholesterol)
> myresult
Call:
   aov(formula = response ~ trt, data = cholesterol)

Terms:
                      trt Residuals
Sum of Squares  1351.3690  468.7504
Deg. of Freedom         4        45

Residual standard error: 3.227488
Estimated effects may be unbalanced
> summary(myresult)
            Df Sum Sq Mean Sq F value   Pr(>F)   
trt          4 1351.4   337.8   32.43 9.82e-13 ***
Residuals   45  468.8    10.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> TukeyHSD(myresult)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = response ~ trt, data = cholesterol)

$trt
                  diff        lwr       upr     p adj
2times-1time   3.44300 -0.6582817  7.544282 0.1380949
4times-1time   6.59281  2.4915283 10.694092 0.0003542
drugD-1time    9.57920  5.4779183 13.680482 0.0000003
drugE-1time   15.16555 11.0642683 19.266832 0.0000000
4times-2times  3.14981 -0.9514717  7.251092 0.2050382
drugD-2times   6.13620  2.0349183 10.237482 0.0009611
drugE-2times  11.72255  7.6212683 15.823832 0.0000000
drugD-4times   2.98639 -1.1148917  7.087672 0.2512446
drugE-4times   8.57274  4.4714583 12.674022 0.0000037
drugE-drugD    5.58635  1.4850683  9.687632 0.0030633
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Archiver|手机版|小黑屋|R语言中文网

GMT+8, 2024-11-25 18:53 , Processed in 0.033225 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表