【R<-练习】双因子方差分析

对于两因素的方差分析, 基本思想和方法与单因素的方差分析相似, 前提 条件仍然是要满足独立、正态、方差齐性. 所不同的是在双因素方差分析中, 有时会出现交互作用, 即二因素的不同水平交叉搭配对指标产生影响. 我们先 讨论无交互作用的双因素方差分析.

1 无交互作用的方差分析

在R软件中, 方差分析函数aov( )既适合于单因素方差分析, 也同样适用 于双因素方差分析, 其中方差模型公式为x~A+B, 加号表示两个因素具有 可加的. 下面用一个例子来说明.

例 题1
原来检验果汁中含铅量有三种方法A1、A2、A3, 现研究出另 一种快速检验法A4, 能否用A4代替前三种方法, 需要通过实验考察. 观察的对 象是果汁, 不同的果汁当做不同的水平: B1为苹果, B2为葡萄汁, B3为西红柿 汁, B4为苹果饮料汁, B5桔子汁, B6菠萝柠檬汁. 现进行双因素交错搭配试验, 即用四种方法同时检验每一种果汁, 其检验结果如表8.4所示. 问因素A(检验方 法)和B(果汁品种) 对果汁的含铅量是否有显著影响?

统计表

解答:

# 将数据导入到命名为juice的数据框中
 juice<-data.frame( X = c(0.05, 0.46, 0.12, 0.16, 0.84, 1.30, 0.08, 0.38, 0.4, 0.10, 0.92, 1.57, 0.11, 0.43, 0.05, 0.10, 0.94, 1.10, 0.11, 0.44, 0.08, 0.03, 0.93, 1.15), A = gl(4, 6), B = gl(6, 1, 24) )
> juice
      X A B
1  0.05 1 1
2  0.46 1 2
3  0.12 1 3
4  0.16 1 4
5  0.84 1 5
6  1.30 1 6
7  0.08 2 1
8  0.38 2 2
9  0.40 2 3
10 0.10 2 4
11 0.92 2 5
12 1.57 2 6
13 0.11 3 1
14 0.43 3 2
15 0.05 3 3
16 0.10 3 4
17 0.94 3 5
18 1.10 3 6
19 0.11 4 1
20 0.44 4 2
21 0.08 4 3
22 0.03 4 4
23 0.93 4 5
24 1.15 4 6

注: 这里函数gl( )用来给出因子水平
其调用格式为 gl(n, k, length=n*k, labels=1:n, ordered=FALSA)
说明: n是水平数, k是每一水平上的重复次数, length是总观测值数, ordered指明各水平是否先排序.

# 接下来作双因子方差分析
> juice.aov <- aov(X~A+B, data = juice)
> summary(juice.aov)
            Df Sum Sq Mean Sq F value Pr(>F)    
A            3  0.057  0.0190   1.629  0.225    
B            5  4.902  0.9804  83.976  2e-10 ***
Residuals   15  0.175  0.0117                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

结论: p值说明果汁品种(因素B)对含铅量有显著影响, 而没有充分理由说明检 验方法(因素A)对含铅量有显著影响.

# 最后用函数bartlett.test( )分别对因素A和因素B作方差的齐性检验:

> bartlett.test(X~A,data = juice)

    Bartlett test of homogeneity of variances

data:  X by A
Bartlett's K-squared = 0.26802, df = 3, p-value = 0.9659

> bartlett.test(X~B,data = juice)

    Bartlett test of homogeneity of variances

data:  X by B
Bartlett's K-squared = 17.422, df = 5, p-value = 0.003766

结论: 对因素A, p值(0.966)远大于0.05, 接受原假设, 认为因素A的各水平下 的数据是等方差的; 对因素B, p值(0.003766)小于0.05, 拒绝原假设, 即认为因 素B不满足方差齐性要求

2 有交互作用的方差分析

R软件中仍用函数aov( )进行有交互作用的方差分析, 但其中的方差模型格式 为x A+B+A : B. 下面用一个例子来全面展示有交互作用方差分析过程.

例 2
有一个关于检验毒品强弱的试验, 给48只老鼠注射I、II、III三 种毒药(因素A), 同时有A、B、C、D 4种治疗方案(因素B), 这样的试验在每一 种因素组合下都重复四次测试老鼠的存活时间, 数据如表8.5所示. 试分析毒药 和治疗方案以及它们的交互作用对老鼠存活时间有无显著影响.

统计表2

首先以数据框形式输入数据, 并用函数plot( )作图.

# 将统计表的数据输入到命名为rats的数据框中
> rats<-data.frame( Time=c(0.31, 0.45, 0.46, 0.43, 0.82, 1.10, 0.88, 0.72, 0.43, 0.45, 0.63, 0.76, 0.45, 0.71, 0.66, 0.62, 0.38, 0.29, 0.40, 0.23, 0.92, 0.61, 0.49, 1.24, 0.44, 0.35, 0.31, 0.40, 0.56, 1.02, 0.71, 0.38, 0.22, 0.21, 0.18, 0.23, 0.30, 0.37, 0.38, 0.29, 0.23, 0.25, 0.24, 0.22, 0.30, 0.36, 0.31, 0.33),
+ toxicant=gl(3,16,48,labels = c("I","II","III")),
+ cure=gl(4,4,48,labels = c("A","B","C","D")))
> op <- par(mfrow=c(1,2)) #更改图形显示的界面排版
> plot(Time~toxicant + cure, data = rats)
Hit <Return> to see next plot: 

显示出的统计图:


图2

图2显示两因素的各水平均有较大差异存在.

下面再用函数interaction.plot( )作出交互效应图, 以考查因素之间交互作用是否存在, R程序为

> with(rats,interaction.plot(toxicant, cure, Time, trace.label = "cure"))
> with(rats,interaction.plot(cure, toxicant, Time, trace.label = "toxicant"))

交互效应图如下所示:

图3

两图中的曲线并没有明显的相交情况出现, 因此我们初步认为两个因素没有交互作用

尽管如此, 由于实验误差的存在, 我们用方差分析函数aov( )对此进行确 认, 其中方差模型格式为x~A*�B, 或A+B+A : B, 表示不仅考虑因素A、B各 自的效应, 还考虑两者的交互效应. 若仅考虑A与B的交互效应则方差模型格式为A : B.

> rats.aov <- aov(Time~toxicant+cure+toxicant:cure,data = rats)
> summary(rats.aov)
              Df Sum Sq Mean Sq F value   Pr(>F)    
toxicant       2 1.0356  0.5178  23.225 3.33e-07 ***
cure           3 0.9146  0.3049  13.674 4.13e-06 ***
toxicant:cure  6 0.2478  0.0413   1.853    0.116    
Residuals     36 0.8026  0.0223                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

根据p值知, 因素Toxicant和Cure对Time的影响是高度显著的, 而交互作用对Time的影响却是不显著的.

再 进 一 步 使 用 前 面 的Bartlett检 验 因 素Toxicant和Cure下的数据是否满足方差齐性的要求, R程序如下.

> bartlett.test(Time~toxicant,data=rats)

    Bartlett test of homogeneity of variances

data:  Time by toxicant
Bartlett's K-squared = 25.806, df = 2, p-value = 2.49e-06

> bartlett.test(Time~cure,data = rats)

    Bartlett test of homogeneity of variances

data:  Time by cure
Bartlett's K-squared = 13.055, df = 3, p-value = 0.004519

其中各p值均小于0.05表明在0.05显著性水平下两因素下的方差 不满足齐性的要求, 这与图2是一致的.

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 202,905评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,140评论 2 379
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,791评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,483评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,476评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,516评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,905评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,560评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,778评论 1 296
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,557评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,635评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,338评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,925评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,898评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,142评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,818评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,347评论 2 342