 上传我的文档
 下载
 收藏
该文档贡献者很忙，什么也没留下。
 下载此文档
正在努力加载中...
由ANOVA结果得到差异显著
下载积分：400
内容提示：由ANOVA结果得到差异显著
文档格式：PDF|
浏览次数：48|
上传日期： 21:59:46|
文档星级：
全文阅读已结束，如果下载本文需要使用
 400 积分
下载此文档
阅读此文档的用户还读了
由ANOVA结果得到差异显著
官方公共微信For enrollment inquiries at ACE Sonoma, ACE Marin, ACE Contra Costa, or our new Transition Program, please contact our In-take Coordinator Stacy Kalember at 707-527-7032 or skalember@anovaeducation.org. You may also submit an inquiry via our .&
Please note: Anova's Admin office has recently moved to 220 Concourse Blvd., Santa Rosa, CA 95403
&(ages 18-22)
Loading events...
Please help us give our ACE Sonoma students a playground!&
Looking to make a positive difference?&Current career opportunities:
&- ACE Campuses in Santa Rosa and Concord
&-&ACE Campuses in Santa Rosa and Concord
&- ACE Campuses in Santa Rosa and Concord, as well as San Francisco and the East Bay Area.&
&&Loading videos from YouTube...
Registration is open!
Designed to build understanding about our social world using traditional camp activities, Anova's Super Social Summer Camp teaches students to foster social thinking, build friendships, and recognize different perspectives. Activities include&Drama, Art, Sports, Music, and much more!&
DATES:July 17th & July 21stJuly 24th & July 28thTIME: 8:30 am & 2:30 pmAGES: 5 - 13 years old
Applications:&&&or&ACE Sonoma front&office.&
&&Loading videos from YouTube...
Our&&is a super way for you to help us sustain and grow your child's ACE campus! Sign-up today&to make ongoing monthly donations to support your child's&school.&ACE Campus pages:&
ANOVA ADMINISTRATIVE OFFICE
220 Concourse Blvd., Santa Rosa, CA 95403 - 707.527.7032
Anova is a non-profit 501(c)(3) organization IRS Tax ID 94-3370998
School CalendarFaculty & Staff DirectorySchool News《R语言实战》第三部分第九章-方差分析复习笔记 - 知乎专栏
{"debug":false,"apiRoot":"","paySDK":"/api/js","wechatConfigAPI":"/api/wechat/jssdkconfig","name":"production","instance":"column","tokens":{"X-XSRF-TOKEN":null,"X-UDID":null,"Authorization":"oauth c3cef7c66aa9e6a1e3160e20"}}
{"database":{"Post":{"":{"title":"《R语言实战》第三部分第九章-方差分析复习笔记","author":"realphy","content":"当包含的变量时解释变量时，我们关注的重点通常会从预测转向组间差异的分析，这种分析法称为方差分析（）。换一句话说：方差分析就是通过检验各总体的均值是否相等来判断分类型自变量对数值型因变量是否有显著影响。1 术语以焦虑症治疗为例，现有两种治疗方案：认知行为疗法（CBT）和眼动脱敏再加工法（EMDR）。我们招募10位焦虑症患者作为志愿者，随机分配一半的人接受为期五周的CBT，另外一半接受为期五周的EMDR，设计方案下表所示。在治疗结束时，要求每位患者都填写状态特质焦虑问卷（STAI），也就是一份焦虑度测量的自我评测报告。在这个实验设计中，治疗方案是两水平（CBT、EMDR）的组间因子。之所以称其为组间因子，是因为每位患者都仅被分配到一个组别中，没有患者同时接受CBT和EMDR。表中字母s代表受试者（患者）。STAI是因变量，治疗方案是自变量。由于在每种治疗方案下观测数相等，因此这种设计也称为均衡设计（balanced design）；若观测数不同，则称作非均衡设计（unbalanceddesign）。单因素方差分析（）：指对单因素试验结果进行分析，检验因素对试验结果有无显著性影响的方法，是两个比较的引伸，它是用来检验多个平均数之间的差异，从而确定因素对试验结果有无显著性影响的一种。单因素组内方差分析及重复测量方差分析：在上述的实例中，假如只对CBT的效果感兴趣，则需将10个患者都放在CBT组中，然后在治疗五周和六个月后分别评价疗效，设计方案如下表所示。此时，时间（time）是两水平（五周、六个月）的组内因子。因为每位患者在所有水平下都进行了测量，所以这种统计设计称单因素组内方差分析；又由于每个受试者都不止一次被测量，也称作重复测量方差分析。此时，时间（time）是两水平（五周、六个月）的组内因子。因为每位患者在所有水平下都进行了测量，所以这种统计设计称单因素组内方差分析；又由于每个受试者都不止一次被测量，也称作重复测量方差分析。协方差分析（）：假设招募患者时使用抑郁症的自我评测报告，比如白氏抑郁症量表（BDI），记录了他们的抑郁水平，那么你可以在评测疗法类型的影响前，对任何抑郁水平的组间差异进行统计性调整。本案例中，BDI为协变量，该设计为协方差分析（ANCOVA）。协方差分析是关于如何调节协变量对因变量的影响效应，从而更加有效地分析实验处理效应的一种统计技术，也是对实验进行统计控制的一种综合方差分析和回归分析的方法。协方差是用来度量两个变量之间“协同变异”大小的总体参数，即二个变量相互影响大小的参数，协方差的绝对值越大，二个变量相互影响越大。多元方差分析（）：当因变量不止一个时，以上述实验为例，为增强研究的有效性，可以对焦虑症进行其他的测量（比如家庭评分、医师评分，以及焦虑症对日常行为的影响评价），此时称作多元方差分析。多元协方差分析（）：在多元方差分析的基础上，如果引入协变量，则称作多元协方差分析。以下两个图表，引自能够比较好的说明各种术语，向原作者表示致意。2 ANOVA模型拟合ANOVA模型拟合采用aov()函数，语法如下：aov(formula, data = dataframe)\n其中formula为表达式，表达式在之前的章节中遇到过，这里总结复习一下，R表达式中的符号意义如下：基于上述符号，对于不同种方差分析R中常用研究设计表达式如下：基于上述符号，对于不同种方差分析R中常用研究设计表达式如下：其中，小写字母表示定量变量，大写字母表示组别因子，Subject是对被试者独有的标识变量。其中，小写字母表示定量变量，大写字母表示组别因子，Subject是对被试者独有的标识变量。此外，当因子不止一个时，以上表达式中因子出现的顺序也非常重要。在具体实际操作当中，样本大小越不平衡，效应项的顺序对结果的影响越大。，越基础性的效应越需要放在表达式前面。首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项。在R语言当中，car包中的Anova()函数（在具体实际操作当中，样本大小越不平衡，效应项的顺序对结果的影响越大。，越基础性的效应越需要放在表达式前面。首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项。在R语言当中，car包中的Anova()函数（不要与标准anova()函数混淆）提供了使用类型Ⅱ或类型Ⅲ方法的选项，而aov()函数使用的是类型I方法。3 单因素方差分析单因素方差分析主要比较分类因子定义的不同组别的因变量均值。以multcomp包中的cholesterol数据集为例，研究哪种药物疗法降低胆固醇最多。&install.packages(\"multcomp\")\n& library(multcomp)\n载入需要的程辑包：mvtnorm\n载入需要的程辑包：survival\n载入需要的程辑包：TH.data\n载入需要的程辑包：MASS\n\n载入程辑包：‘TH.data’\n\nThe following object is masked from ‘package:MASS’:\n\n
geyser\n\n& attach(cholesterol)\n& head(cholesterol)\n
trt response\n1 1time
2.7139\n& table(trt)\ntrt\n 1time 2times 4times
10 \n& aggregate(response,by=list(trt),FUN = mean)\n
9.2times 12.37478\n4
drugD 15.36117\n5
drugE 20.94752\n& aggregate(response,by=list(trt),FUN = sd)\n
drugE 3.345003\n& fit_91 &- aov(response ~ trt)\n& summary(fit_91)\n
Df Sum Sq Mean Sq F value
32.43 9.82e-13 ***\nResiduals
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n& library(gplots)\nError in library(gplots) : 不存在叫‘gplots’这个名字的程辑包\n& install.packages(\"gplots\")\n& library(gplots)\n\n载入程辑包：‘gplots’\n\nThe following object is masked _by_ ‘.GlobalEnv’:\n\n
residplot\n\nThe following object is masked from ‘package:stats’:\n\n
lowess\n\n& plotmeans(response ~ trt, xlab = \"Treatment\", ylab = \"Responese\", main = \"Mean Plot\\nwith 95% CI\")\n从输出结果看ANOVA对治疗方式（trt）的F检验非常显著（p&0.0001），说明五种疗法的效果不同。顺便提一句，方差分析主要依靠F分布（从输出结果看ANOVA对治疗方式（trt）的F检验非常显著（p&0.0001），说明五种疗法的效果不同。顺便提一句，方差分析主要依靠F分布（）为概率分布的依据，利用平方和与自由度（）所计算的组间与组内均方（Mean of square）估计出F值，若有显著差异则考量进行或称多重比较。gplots包中的plotmeans()函数可以用来绘制带有置信区间的组均值图形，默认置信区间为95%。ANOVA仅仅表明五种药物疗法效果不同，但是并没有告知哪种疗法与其他疗法不同。多重比较TukeyHSD()函数提供了对各组均值差异的成对检验。& TukeyHSD(fit_91)\n
Tukey multiple comparisons of means\n
95% family-wise confidence level\n\nFit: aov(formula = response ~ trt)\n\n$trt\n
p adj\n2times-1time
3.4282 0.times-1time
6.5092 0.0003542\ndrugD-1time
9.5482 0.0000003\ndrugE-1time
15.832 0.times-2times
3.1092 0.2050382\ndrugD-2times
6.1482 0.0009611\ndrugE-2times
11.7832 0.0000000\ndrugD-4times
2.9672 0.2512446\ndrugE-4times
8.5022 0.0000037\ndrugE-drugD
5.5632 0.0030633\n\n& par(las=2)\n& par(mar=c(5,8,4,2))\n& plot(TukeyHSD(fit_91))\n从结果可以看出，1time和2times的均值差异不显著，而1time和4times间的差异非常显著，从图形的角度，置信置信区间包含0的疗法说明差异不显著。从结果可以看出，1time和2times的均值差异不显著，而1time和4times间的差异非常显著，从图形的角度，置信置信区间包含0的疗法说明差异不显著。multcomp包中的glht()函数提供了多重均值更为全面的方法。& par(mar=c(5,4,6,2))\n& tuk &- glht(fit_91, linfct = mcp(trt=\"Tukey\"))\n& par(las = 1)\n& plot(cld(tuk, level = .05),col = \"red\")\n注意，记得par(las = 1)函数将轴标签旋转回来。同样的，我们对于结果的信心依赖于做统计检验时数据满足假设条件的程度。单因素方差分析中，我们假设因变量服从正态分布，各组方差相等。使用Q-Q图来检验正态性假设：& library(car)\n& qqPlot(lm(response ~ trt, data=cholesterol),\nsimulate=TRUE, main=\"Q-Q Plot\", labels=FALSE)\nqqPlot()要求用lm()拟合，从结果上看，数据落在95%的置信区间内，满足正态性假设。此外，R还提供了一些可用来做方差齐性检验的函数，比如bartlett.test()函数（）、fligner.test()函数（Fligner-Killeen检验）以及HH包中的hov()函数（ ）。考虑到，方差齐性分析对离群点非常敏感，可利用car包中的outlierTest(0函数检测离群点。& outlierTest(fit_91)\n\nNo Studentized residuals with Bonferonni p & 0.05\nLargest |rstudent|:\n
rstudent unadjusted p-value Bonferonni p\n19 2.251149
NA\n& \n4 单因素协方差分析单因素协方差分析（ANCOVA）扩展了单因素方差分析（ANOVA），包含一个或多个定量的协变量。以multcomp包中的litter数据集为例，四组怀孕小鼠接受不同剂量（0、5、50或500）的药物处理，怀孕时间为协变量，产下幼崽的体重均值为因变量。& data(litter,package = \"multcomp\")\n& attach(litter)\n& table(dose)\ndose\n
50 500 \n 20
17 \n& aggregate(weight, by = list(dose), FUN = mean)\n
0 32.30850\n2
5 29.30842\n3
50 29.86611\n4
500 29.64647\n& fit_93 &- aov(weight ~ gesttime + dose)\n& summary(fit_93)\n
Df Sum Sq Mean Sq F value
\ngesttime
8.049 0.00597 **\ndose
2.739 0.04988 * \nResiduals
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n& \n我们不难发现，协方差分析同样使用aov()函数，但是formula中的变量顺序是将协变量放到实际实际水平变量前面。另外，使用aggregate()函数求均值时，由于均值收到协变量的影响，需要用effects包中的effects()函数计算去除协变量效应后的组均值。& library(effects)\n\n载入程辑包：‘effects’\n\nThe following object is masked from ‘package:car’:\n\n
Prestige\n\n& effect(\"dose\",fit_93)\n\n dose effect\ndose\n
500 \n32.72 30.60 \n& \n结果表明：(a)怀孕时间与幼崽出生体重相关；(b)控制怀孕时间，药物剂量与出生体重相关。控制怀孕时间，确实发现每种药物剂量下幼崽出生体重均值不同。多重比较使用multcomp包中的glht()函数，mcp参数由rbind(c(3,-1,-1,-1))给出，意思为用量为0的和其它三组用量不为0的作比较，也可以使得c(1,-1,0,0)得到用量0和用量5的做比较。代码及结果如下：& library(multcomp)\n& comtrast &- rbind(\"no drrug vs. drug\" = c(3, -1, -1, -1))\n& summary(glht(fit_93,linfct = mcp(dose = comtrast)))\n\n\t Simultaneous Tests for General Linear Hypotheses\n\nMultiple Comparisons of Means: User-defined Contrasts\n\n\nFit: aov(formula = weight ~ gesttime + dose)\n\nLinear Hypotheses:\n
Estimate Std. Error t value Pr(&|t|)
\nno drrug vs. drug == 0
0.012 *\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n(Adjusted p values reported -- single-step method)\n最后用HH保重的ancova()函数绘制因变量、协变量和因子之间的关系图。代码和结果如下：&library(HH)\n载入需要的程辑包：lattice\n载入需要的程辑包：grid\n载入需要的程辑包：latticeExtra\n载入需要的程辑包：RColorBrewer\n载入需要的程辑包：gridExtra\n\n载入程辑包：‘HH’\n\nThe following object is masked _by_ ‘.GlobalEnv’:\n\n
residplot\n\nThe following objects are masked from ‘package:car’:\n\n
logit, vif\n\nThe following object is masked from ‘package:gplots’:\n\n
residplot\n\n& ancova(weight ~ gesttime + dose, data = litter)\nAnalysis of Variance Table\n\nResponse: weight\n
Sum Sq Mean Sq F value
\ngesttime
134.30 134.304
8.971 **\ndose
2.883 * \nResiduals 69 .685
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n& \n结果表明，在不同的剂量下，重量随着怀孕时间增加而增加。，还可以看到0剂量组截距项最大，5剂量组截距项最小。由于你上面的设置，直线会保持平行，若用ancova(weight ~ gesttime*dose)，生成的图形将允许斜率和截距项依据组别而发生变化，这对可视化那些违背回归斜率同质性的实例非常有用。ANCOVA的假设验证除了正态性和同方差性，还需要进行回归斜率检验。本例中，假定四个处理组通过怀孕时间来预测出生体重的回归斜率都相同。ANCOVA模型包含怀孕时间×剂量的交互项时，可对回归斜率的同质性进行检验。交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。ANCOVA的假设验证除了正态性和同方差性，还需要进行回归斜率检验。本例中，假定四个处理组通过怀孕时间来预测出生体重的回归斜率都相同。ANCOVA模型包含怀孕时间×剂量的交互项时，可对回归斜率的同质性进行检验。交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。& fit_932 &- aov(weight ~ gesttime*dose, data = litter)\n& summary(fit_932)\n
Df Sum Sq Mean Sq F value
\ngesttime
8.289 0.00537 **\ndose
2.821 0.04556 * \ngesttime:dose
1.684 0.17889
\nResiduals
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n& \n5 双因素方差分析与单因素协方差分析类似，双因素方差分析的区别就在于他们有交互效应，这个与协方差分析中后面的回归斜率检验有些许类似，区别在于交互效应是否显著。以基础安装中的ToothGrowth数据集为例，随机分配60只豚鼠，分别采用两种喂食方法（橙汁或维生素C），各喂食方法中抗坏血酸含量有三种水平（0.5mg、1mg或2mg），每种处理方式组合都被分配10只豚鼠。其中牙齿长度为因变量。& attach(ToothGrowth)\nThe following object is masked from litter:\n\n
dose\n\n& table(supp, dose)\n
dose\nsupp 0.5
10 10 10\n
10 10 10\n& aggregate(len, by = list(supp,dose), FUN = mean)\n
Group.1 Group.2
0.5 13.23\n2
1.0 22.70\n4
1.0 16.77\n5
2.0 26.06\n6
2.0 26.14\n& aggregate(len, by = list(supp,dose), FUN = sd)\n
Group.1 Group.2
2.0 4.797731\n& dose &- factor(dose)\n& fit_96 &- aov(len ~ supp*dose)\n& summary(fit_96)\n
Df Sum Sq Mean Sq F value
15.572 0.000231 ***\ndose
& 2e-16 ***\nsupp:dose
4.107 0.021860 *
\nResiduals
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n结果表明这是一个均衡设计，aggregate()函数获得各单元的均值和标准差。summary()函数的结果表明主效应和交互效应都非常显著。双因素方差分析可视化有以下三种方法：& interaction.plot(dose, supp, len, type=\"b\",\n
col=c(\"red\",\"blue\"), pch=c(16, 18),\n
main = \"Interaction between Dose and Supplement Type\")\n& library(gplots)\n& plotmeans(len ~ interaction(supp, dose, sep = \" \"),\n+
connect = list(c(1,3,5),c(2,4,6)),\n+
col = c(\"red\", \"darkgreen\"),\n+
main = \"Interaction Plot with 95% CIs\",\n+
xlab = \"Treatment and Dose Combination\")\n& \n& library(HH)\n& interaction2wt(len~supp*dose)\n三种绘图方法中，我更推荐HH包中的interaction2wt()函数，因为它能展示任意复杂度设计（双因素方差分析、三因素方差分析等）的主效应（箱线图）和交互效应。6 重复测量方差分析所谓重复测量方差分析，即测试对象会被重复测量其结果，打个比方就很容易明白了，比如，两个人在同一年中各月的身高变化，这里就包含了一个因变量（身高），一个组间因子（人，有两个，所以是两个水平），一个组内因子（时间，12个月，就是12个水平），此时，对这两个人来说（测试对象），他们每个人要测试12次，即组内因子的水平数，所以叫重复测量方差分析。 书中实例是研究生命系统的生理和生化过程如何响应环境因素的变异，对美国北方和南方牧草类植物Echinochloa crus-galli的寒冷容忍度研究结果进行了分析，比较了在某浓度二氧化碳的环境中，对寒带植物与非寒带植物的光合作用率。#将conc转换成因子变量\n& CO2$conc &- factor(CO2$conc)\n#去除寒带植物的数据\n& w1b1 &- subset(CO2, Treatment == 'chilled')\n#aov拟合，formula带Error\n& fit_97 &- aov(uptake ~ conc*Type + Error(Plant/(conc)),w1b1)\n#显示拟合结果\n& summary(fit_97)\n\nError: Plant\n
Df Sum Sq Mean Sq F value
60.41 0.00148 **\nResiduals
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nError: Plant:conc\n
Df Sum Sq Mean Sq F value
52.52 1.26e-12 ***\nconc:Type
15.30 3.75e-07 ***\nResiduals 24
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n#设置坐标轴旋转\n& par(las=2)\n#图面布局\n& par(mar=c(10,4,4,2))\n#绘制交互效应\n& with(w1b1,interaction.plot(conc,Type,uptake,\n+
type=\"b\",col = c(\"red\",\"blue\"),pch = c(16,18),\n+
main= \"Interaction Plot for Plant Type and Concentration\"))\n#绘制箱型交互图\n& boxplot(uptake ~ Type*conc, data = w1b1,col=(c(\"gold\",\"green\")),\n+
main = \"Chilled Quebec and Mississippi Plants\",\n+
ylab = \"Carbon dioxide uptake rate (umon/m^2 sec)\")\n从以上任意一幅图都可以看出，魁北克省的植物比密西西比州的植物二氧化碳吸收率高，而且随着CO2浓度的升高，差异越来越明显。本例使用了传统的重复测量方差分析。该方法假设任意组内因子的协方差矩阵为球形，并且任意组内因子两水平间的方差之差都相等。但在现实中这种假设不可能满足，于是衍生了一系列备选方法：使用lme4包中的lmer()函数拟合线性混合模型（Bates，2005） 使用car包中的Anova()函数调整传统检验统计量以弥补球形假设的不满足（例如Geisser-Greenhouse校正） 使用nlme包中的gls()函数拟合给定方差-协方差结构的广义最小二乘模型（UCLA，2009） 使用多元方差分析对重复测量数据进行建模（Hand，1987）。 7 多元方差分析多元方差分析（MANOVA）用于应对因变量（结果变量）不止一个时的情形。以MASS包中的UScereal数据集为例，研究美国谷物中的卡路里、脂肪和糖含量是否会因为储存架位置的不同而发生变化；其中1代表底层货架，2代表中层货架，3代表顶层货架。卡路里、脂肪和糖含量是因变量，货架是三水平（1、2、3）的自变量。& library(MASS)\n& attach(UScereal)\n& shelf &- factor(shelf)\n& y &- cbind(calories, fat, sugars)\n& aggregate(y, by = list(shelf), FUN = mean)\n
Group.1 calories
1 119.493\n2
2 129.670\n3
3 180.821\n& cov(y)\n
sugars\ncalories
60..380317\nfat
3.995474\nsugars
34.050018\n& fit_98 &- manova(y ~ shelf)\n& summary(fit_98)\n
Df Pillai approx F num Df den Df
122 0.0001015 ***\nResiduals 62
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n& \n需要注意的是，shelf需要因子化，manova()传入的formula参数～号左边的y必须用cbind()函数进行合成。从结果上看，F值显著，三个组的营养成分测量值不同。在多元检验是显著的前提下，可以使用summary.aov()函数对每一个变量再做单因素方差分析。单因素单因素多元方差分析有两个前提假设，一个是多元正态性，一个是方差-协方差矩阵同质性。第一个假设即指因变量组合成的向量服从一个多元正态分布。多元正态性可以由QQ图检验，理论补充及代码如下：#若有一个p*1的多元正态随机向量x，均值为μ，协方差矩阵为Σ，那么x与μ的马氏距离的平方服从自由度为p的卡方分布。Q-Q图展示卡方分布\n& center &- colMeans(y)\n& n &- nrow(y)\n& p &- ncol(t)\n& cov &- cov(y)\n& p &- ncol(y)\n& d &- mahalanobis(y,center,cov)\n& coord &- qqplot(qchisq(ppoints(n),df = p),\n+
d, main = \"Q-Q Plot Assessing Multivariate Normality\",\n+
ylab = \"Mahalanobis D2\")\n& abline(a=0, b=1)\n& identify(coord$x, coord$y, labels=row.names(UScereal))\n警告: 已经找到了最近的点\n警告: 已经找到了最近的点\n警告: 已经找到了最近的点\n警告: 没有点在0.25英尺内\n警告: 已经找到了最近的点\n警告: 已经找到了最近的点\n[1] 62 63 64 65\n& \n方差-协方差矩阵同质性即指各组的协方差矩阵相同，通常可用Box’s M检验来评估该假设，不过遗憾的是目前R中并没有Box‘s M函数。mvoutlier包中的ap.plot()函数可以用来检验多元离群点。代码和图形如下：#需要提前按照mvoutlier包\n& library(mvoutlier)\n载入需要的程辑包：sgeostat\nsROC 0.1-2 loaded\n& outliers &- aq.plot(y)\nProjection to the first and second robust principal components.\nProportion of total variation (explained variance): 0.9789888\n& outliers\n$outliers\n [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
TRUE\n[12]
TRUE FALSE FALSE FALSE FALSE FALSE
TRUE FALSE FALSE FALSE FALSE\n[23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
TRUE FALSE\n[34] FALSE FALSE FALSE FALSE FALSE FALSE
TRUE FALSE FALSE FALSE
TRUE\n[45] FALSE FALSE FALSE FALSE FALSE
TRUE FALSE FALSE FALSE FALSE FALSE\n[56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE\n\n& \n在多元正态性或方差-协方差同质性或担心离群点时，可以使用稳健多元方差分析或非参数多元方差分析，稳健多元单因素方差分析可以通过rrcov包中Wilks.test(y.dataframe,x,method=”mcp”)，另外，vegan包中的adonis()函数则提供了非参数MANOVA的等同形式，代码和结果如下：在多元正态性或方差-协方差同质性或担心离群点时，可以使用稳健多元方差分析或非参数多元方差分析，稳健多元单因素方差分析可以通过rrcov包中Wilks.test(y.dataframe,x,method=”mcp”)，另外，vegan包中的adonis()函数则提供了非参数MANOVA的等同形式，代码和结果如下：& library(rrcov)\n载入需要的程辑包：robustbase\n\n载入程辑包：‘robustbase’\n\nThe following object is masked from ‘package:survival’:\n\n
heart\n\nScalable Robust Estimators with High Breakdown Point (version 1.4-3)\n\n& Wilks.test(y,shelf,method = \"mcd\")\n\n#Windows运行超级慢\n& Wilks.test(y, shelf, method = \"mcd\")\n\n\tRobust One-way MANOVA (Bartlett Chi2)\n\ndata:
x\nWilks' Lambda = 0.51073, Chi2-Value = 22.8480, DF = 4.8047,\np-value = 0.0003013\nsample estimates:\n
sugars\n1 119.143\n2 128.533\n3 160.646\n& library(vegan)\n载入需要的程辑包：permute\nThis is vegan 2.4-2\n& adonis(y~shelf)\n\nCall:\nadonis(formula = y ~ shelf) \n\nPermutation: free\nNumber of permutations: 999\n\nTerms added sequentially (first to last)\n\n
Df SumsOfSqs
MeanSqs F.Model
0.002 **\nResiduals 62
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\nWilks.test()函数在Windows系统下运行超级慢，E31230 16G内存，大约1分28秒，CPU使用率大约15%，回头我可以在Ubuntu系统下测试下做个对比。从结果来看，这两个方差分析的结果都是显著的。再一次验证了货架的不同层对谷物的三种物质的含量是有影响的。 8 用回归来做ANOVA前面提到ANOVA和回归都是广义线性模型的特例，本章所有的设计都可以用lm()函数来进行分析，不过由于回归要求预测变量是数值型，当lm()函数碰到因子时，它会用一系列与因子水平相对应的数值型对照变量来代替因子。如果因子有k个水平，将会创建k–1个对照变量。R提供了五种创建对照变量的内置方法。以第三节的单因素ANOVA问题为例，用lm()函数比较五种降低胆固醇药物疗法（trt）的影响。以第三节的单因素ANOVA问题为例，用lm()函数比较五种降低胆固醇药物疗法（trt）的影响。& library(multcomp)\n& levels(cholesterol$trt)\n[1] \"1time\"
\"2times\" \"4times\" \"drugD\"
\"drugE\" \n& fit.aov &- aov(response ~ trt, data = cholesterol)\n& summary(fit.aov)\n
Df Sum Sq Mean Sq F value
32.43 9.82e-13 ***\nResiduals
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n& \n& fit.lm &- lm(response ~ trt, data = cholesterol)\n& summary(fit.lm)\n\nCall:\nlm(formula = response ~ trt, data = cholesterol)\n\nResiduals:\n
Max \n-6.2 -0.1
6.6008 \n\nCoefficients:\n
Estimate Std. Error t value Pr(&|t|)
\n(Intercept)
5.665 9.78e-07 ***\ntrt2times
\ntrt4times
4.568 3.82e-05 ***\ntrtdrugD
6.637 3.53e-08 ***\ntrtdrugE
10.507 1.08e-13 ***\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nResidual standard error: 3.227 on 45 degrees of freedom\nMultiple R-squared:
0.7425,\tAdjusted R-squared:
0.7196 \nF-statistic: 32.43 on 4 and 45 DF,
p-value: 9.819e-13\n不难发现，方差分析所得的F检验的P值与回归得到的P值是一样的，具体原因可以参考这篇文献，。通常，解释变量有连续的，只能用回归分析。在回归分析时的变量对照，可以通过设置contrasts选项来修改lm()中默认的对照方法。fit.lm &- lm(response ~ trt, data=cholesterol, contrasts=\"contr.helmert\")\n上面的代码将使用Helmert对照。也可以通过options()函数修改R会话当中的默认对照方法，举个例子：options(contrasts = c(\"contr.SAS\", \"contr.helmert\"))\n9 小结我们回顾了基本实验和准实验设计的分析方法，包括ANOVA/ANCOVA/MANOVA。然后通过组内和组间设计的示例介绍了基本方法的使用，如单因素ANOVA、单因素ANCOVA、双因素ANOVA、重复测量ANOVA和单因素MANOVA及各模型的假设检验。下一章，我们将介绍功效分析，功效分析可以帮助我们在给定置信度的情况下，判断达到要求效果所需的样本大小。","updated":"T05:51:53.000Z","canComment":false,"commentPermission":"anyone","commentCount":0,"collapsedCount":0,"likeCount":4,"state":"published","isLiked":false,"slug":"","isTitleImageFullScreen":false,"rating":"none","titleImage":"/v2-e3f3897dbdb0ae83df5b16_r.png","links":{"comments":"/api/posts//comments"},"reviewers":[],"topics":[{"url":"/topic/","id":"","name":"R（编程语言）"},{"url":"/topic/","id":"","name":"方差分析"}],"adminClosedComment":false,"titleImageSize":{"width":478,"height":601},"href":"/api/posts/","excerptTitle":"","column":{"slug":"datacruiser","name":"数海拾荒"},"tipjarState":"inactivated","annotationAction":[],"sourceUrl":"","pageCommentsCount":0,"snapshotUrl":"","publishedTime":"T13:51:53+08:00","url":"/p/","lastestLikers":[{"bio":null,"isFollowing":false,"hash":"a3dedcf1f7cc","uid":480400,"isOrg":false,"slug":"mu-shui-69-58","isFollowed":false,"description":"","name":"暮水","profileUrl":"/people/mu-shui-69-58","avatar":{"id":"da8e974dc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},{"bio":"ggplot(R,aes(zhao,he))+geom_爱(qiu)","isFollowing":false,"hash":"98b363cdaaa0b6fcf924b2","uid":529200,"isOrg":false,"slug":"sa-ma-er-han-han-wang","isFollowed":false,"description":"环境科学本科，生态学硕士，R语言爱好者，ggplot2数据可视化狂热爱好者","name":"撒马尔罕汗王","profileUrl":"/people/sa-ma-er-han-han-wang","avatar":{"id":"v2-47e620e3adc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},{"bio":null,"isFollowing":false,"hash":"dca55811b4fab822a98c370","uid":819800,"isOrg":false,"slug":"twq-80","isFollowed":false,"description":"","name":"tworf","profileUrl":"/people/twq-80","avatar":{"id":"da8e974dc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},{"bio":"量化，模拟，投资，生物信息","isFollowing":false,"hash":"784c49bba61f1b30e19ab","uid":101800,"isOrg":false,"slug":"simingnie","isFollowed":false,"description":"","name":"SimingNie","profileUrl":"/people/simingnie","avatar":{"id":"da8e974dc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false}],"summary":"当包含的变量时解释变量时，我们关注的重点通常会从预测转向组间差异的分析，这种分析法称为方差分析（）。换一句话说：方差分析就是通过检验各总体的均值是否相等来判断分类型自变量对数值型因变量是否有显著影响。1 术语以焦虑症治疗为例，现有两种…","reviewingCommentsCount":0,"meta":{"previous":{"isTitleImageFullScreen":false,"rating":"none","titleImage":"/v2-5f310e9e1d3129edafbc76a8f7841c6f_r.jpg","links":{"comments":"/api/posts//comments"},"topics":[{"url":"/topic/","id":"","name":"R（编程语言）"}],"adminClosedComment":false,"href":"/api/posts/","excerptTitle":"","author":{"bio":"数据分析/机器学习 修炼中","isFollowing":false,"hash":"20c04cea624a1ad","uid":28,"isOrg":false,"slug":"realphy","isFollowed":false,"description":"人生没有起跑线，终生学习践行者","name":"Datacruiser","profileUrl":"/people/realphy","avatar":{"id":"v2-6f18fa1ab887d9f7a722596","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},"column":{"slug":"datacruiser","name":"数海拾荒"},"content":"之前在研究生毕业论文当中用SPSS采用多元线性及非线性模型用于预测配煤动力学特性，本节主要结合R语言提供的函数和工具重现一下。1 几个概念：单煤：就是从一个矿区开采出来的煤炭，一般性质相近；配煤：为了取得特定的特性，由不同种单煤按照不同的比例混合而成，由于煤炭性质的复杂性，通常配煤的动力学特性无法由各种单煤加权平均求得；煤的工业分析（）：煤的工业分析是指包括煤的水分（M ）、灰分（A ）、挥发分（V ）和固定碳（Fc )四个分析项目指标的测定的总称。煤的工业分析是了解煤质特性的主要指标，也是评价煤质的基本依据。通常煤的水分、灰分、挥发分是直接测出的，而固定碳是用差减法计算出来的。广义上讲，煤的工业分析还包括煤的全硫分和发热量的测定， 又叫煤的全工业分析。本文的实验数据介于狭义与广义之间，测量了发热量，但是全硫在元素分析中进行。煤的元素分析（）：煤的元素分析是对煤中的元素含量进行检测和分析（一般用质量百分数表示），包括常规的C、H、O、N、S、Al、Si、Fe、Ca等元素含量，还可检测煤中的痕量元素包括Ti、Na、K等。本文的实验数据主要测量了最主要的C、H、O、N、S五种。煤的着火温度：煤的着火点是随着热力条件变化而变化的，并不是一个物理常数，只是在一定条件下得出的相对特征值。通常把由缓慢氧化状态转变为高速燃烧状态的瞬间过程称为着火，此时的温度即称为着火温度（Ignition Temperature）。煤的活化能（Activation Energy）：即让某一种煤开始发生（通常指燃烧）所需要克服的能量障碍。2 数据描述：有16种单煤及48种这些单煤按不同比例混合而成配煤的工业分析以及元素分析数据，并通过热重分析、热重微分分析及差热分析得出的曲线结合Arrhenius方程及Doyle单曲线法求得各个不同煤样的着火温度（）和活化能。这样一共64组数据。考虑到从实验的角度求取着火温度及活化能需要特殊的仪器并耗费较长的时间，希望能基于上述64组数据找到一个用于预测工业分析、元素分析已知煤样的动力学特性，包括着火温度、活化能等，特别是在混配煤动力学特性预测方面进行应用。3 数据导入：& library(xlsx)\n载入需要的程辑包：rJava\n载入需要的程辑包：xlsxjars\n\n载入程辑包：‘xlsx’\n\nThe following objects are masked from ‘package:openxlsx’:\n\n
createWorkbook, loadWorkbook, read.xlsx, saveWorkbook,\n
write.xlsx\n\n& filepath &- \"G:/DataCruiser/workspace/Coal Analysis/Coal Analysis data.xlsx\"\n& CoalData &- read.xlsx(filepath, 1, encoding = \"UTF-8\")\n& class(CoalData)\n[1] \"data.frame\"\n& head(CoalData)\n
煤样编号 煤样名称 Ignition_Temp Activation_Ene Mad.
Vad. FCad.\n1
.85 23.68 26.98 45.49\n2
.00 14.62 10.44 73.94\n3
.50 14.42 26.06 57.02\n4
.47 14.19 30.02 46.32\n5
.21 18.26 29.39 50.14\n6
.62 25.47 23.32 49.59\n
Qnet.ad.J.g.
Cad. Had. Nad. St.add.
.93 3.42 1.12
.90 3.14 1.15
.55 3.92 0.81
.44 2.92 0.83
0.44 13.71\n5
.13 3.74 1.18
0.71 10.77\n6
.26 3.32 1.09
6.72\n4 数据重命名：着火温度：ignitionTemperature活化能：activationEnergy水分：moisture灰分：ash挥发分：volatileMatter固定碳：fixedCarbon发热量：LHV碳：carbon氢：hydrogen氮：nitrogen硫：sulfur& names(CoalData) &- c(\"sampleNO\",\"sampleName\",\"ignitionTemperature\",\"activationEnergy\",\"moistrue\",\n+
\"ash\",\"volatileMatter\",\"fixedCarbon\",\"LHV\",\"carbon\",\"hydrogen\",\"nitrogen\",\"sulfur\",\"oxygen\")\n& head(CoalData)\n
sampleNO sampleName ignitionTemperature activationEnergy moistrue
3.85 23.68\n2
1.00 14.62\n3
2.50 14.42\n4
9.47 14.19\n5
2.21 18.26\n6
1.62 25.47\n
volatileMatter fixedCarbon
LHV carbon hydrogen nitrogen sulfur\n1
6.72\n4 数据分组：按不同因变量（ignitionTemperature和activationEnergy）进行数据分组。& IT_Estimate &- as.data.frame(CoalData[,c(\"ignitionTemperature\",\"moistrue\",\n+
\"ash\",\"volatileMatter\",\"fixedCarbon\",\"LHV\",\"carbon\",\"hydrogen\",\"nitrogen\",\"sulfur\",\"oxygen\")])\n& head(IT_Estimate)\n
ignitionTemperature moistrue
ash volatileMatter fixedCarbon
3.85 23.68
1.00 14.62
2.50 14.42
9.47 14.19
2.21 18.26
1.62 25.47
49.59 25702.6\n
carbon hydrogen nitrogen sulfur oxygen\n1
6.72\n& AE_Estimate &- as.data.frame(CoalData[,c(\"activationEnergy\",\"moistrue\",\n+
\"ash\",\"volatileMatter\",\"fixedCarbon\",\"LHV\",\"carbon\",\"hydrogen\",\"nitrogen\",\"sulfur\",\"oxygen\")])\n& head(AE_Estimate)\n
activationEnergy moistrue
ash volatileMatter fixedCarbon
3.85 23.68
1.00 14.62
2.50 14.42
9.47 14.19
2.21 18.26
1.62 25.47
49.59 25702.6\n
carbon hydrogen nitrogen sulfur oxygen\n1
6.72\n5 数据相关性检验：这里我们采用Pearson相关性分析，需要调用psych包中的corr.test()函数，结果如下：& corr.test(IT_Estimate, use = \"complete\")\nCall:corr.test(x = IT_Estimate, use = \"complete\")\nCorrelation matrix \n
ignitionTemperature moistrue
ash volatileMatter\nignitionTemperature
-0.65\nmoistrue
1.00 -0.42
-0.20\nvolatileMatter
0.31 -0.20
1.00\nfixedCarbon
-0.12 -0.59
-0.65\nLHV
-0.14 -0.79
-0.14\ncarbon
-0.06 -0.82
-0.22\nhydrogen
-0.33 -0.25
0.60\nnitrogen
-0.49 -0.23
0.00\nsulfur
0.30\noxygen
0.61 -0.29
fixedCarbon
LHV carbon hydrogen nitrogen sulfur\nignitionTemperature
-0.07\nmoistrue
-0.12 -0.14
-0.11\nash
-0.59 -0.79
0.29\nvolatileMatter
-0.65 -0.14
0.30\nfixedCarbon
-0.44\nLHV
-0.31\ncarbon
-0.39\nhydrogen
0.26\nnitrogen
-0.06\nsulfur
-0.44 -0.31
1.00\noxygen
-0.40 -0.16
oxygen\nignitionTemperature
-0.57\nmoistrue
-0.29\nvolatileMatter
0.67\nfixedCarbon
-0.40\nLHV
-0.16\ncarbon
-0.27\nhydrogen
0.04\nnitrogen
-0.40\nsulfur
-0.04\noxygen
1.00\nSample Size \n[1] 64\nProbability values (Entries above the diagonal are adjusted for multiple tests.) \n
ignitionTemperature moistrue
ash volatileMatter\nignitionTemperature
0.00\nmoistrue
1.00\nvolatileMatter
0.00\nfixedCarbon
0.27\ncarbon
0.08\nhydrogen
0.00\nnitrogen
0.99\nsulfur
0.02\noxygen
fixedCarbon
LHV carbon hydrogen nitrogen sulfur\nignitionTemperature
1.00\nmoistrue
0.46\nvolatileMatter
0.40\nfixedCarbon
0.33\ncarbon
0.05\nhydrogen
0.80\nnitrogen
1.00\nsulfur
0.00\noxygen
oxygen\nignitionTemperature
0.00\nmoistrue
0.51\nvolatileMatter
0.00\nfixedCarbon
1.00\ncarbon
0.68\nhydrogen
1.00\nnitrogen
0.03\nsulfur
1.00\noxygen
0.00\n\n To see confidence intervals of the correlations, print with the short=FALSE option\n&
corr.test(AE_Estimate, use = \"complete\")\nCall:corr.test(x = AE_Estimate, use = \"complete\")\nCorrelation matrix \n
activationEnergy moistrue
ash volatileMatter\nactivationEnergy
-0.40 -0.03
-0.56\nmoistrue
1.00 -0.42
-0.20\nvolatileMatter
0.31 -0.20
1.00\nfixedCarbon
-0.12 -0.59
-0.65\nLHV
-0.14 -0.79
-0.14\ncarbon
-0.06 -0.82
-0.22\nhydrogen
-0.33 -0.25
0.60\nnitrogen
-0.49 -0.23
0.00\nsulfur
0.30\noxygen
0.61 -0.29
fixedCarbon
LHV carbon hydrogen nitrogen sulfur oxygen\nactivationEnergy
-0.46\nmoistrue
-0.12 -0.14
-0.59 -0.79
-0.29\nvolatileMatter
-0.65 -0.14
0.67\nfixedCarbon
-0.40\nLHV
-0.16\ncarbon
-0.27\nhydrogen
0.04\nnitrogen
-0.40\nsulfur
-0.44 -0.31
-0.04\noxygen
-0.40 -0.16
1.00\nSample Size \n[1] 64\nProbability values (Entries above the diagonal are adjusted for multiple tests.) \n
activationEnergy moistrue
ash volatileMatter fixedCarbon\nactivationEnergy
0.00\nmoistrue
0.00\nvolatileMatter
0.00\nfixedCarbon
0.00\ncarbon
0.00\nhydrogen
0.22\nnitrogen
0.01\nsulfur
0.00\noxygen
LHV carbon hydrogen nitrogen sulfur oxygen\nactivationEnergy 0.04
0.01\nmoistrue
0.51\nvolatileMatter
0.00\nfixedCarbon
1.00\ncarbon
0.68\nhydrogen
1.00\nnitrogen
0.03\nsulfur
1.00\noxygen
0.00\n\n To see confidence intervals of the correlations, print with the short=FALSE option\n从相关性检验中可以得出无论是着火温度还是活化能，与水分、挥发分、固定碳、发热量以及氧这五个变量相关较大，因此，以下考虑以此五个自变量采用多元线性回归模型分别对着火温度和活化能进行模拟，首先剔除不相关数据，仅包含我们感兴趣的变量。& AE_Estimate &- as.data.frame(AE_Estimate[,c(\"activationEnergy\",\"moistrue\",
\"volatileMatter\",\"fixedCarbon\",\"LHV\",\"oxygen\")])\n& IT_Estimate &- as.data.frame(IT_Estimate[,c(\"ignitionTemperature\",\"moistrue\",
\"volatileMatter\",\"fixedCarbon\",\"LHV\",\"oxygen\")])\n& \n用car包中的scatterplotMatrix()函数生成散点图矩阵。& library(car)\n\n载入程辑包：‘car’\n\nThe following object is masked from ‘package:psych’:\n\n
logit\n\n& scatterplotMatrix(IT_Estimate, spread = FALSE, smoother.args = list(lty = 2), main = \"Scatter Plot Matrix ignitionTemperature\")& scatterplotMatrix(AE_Estimate, spread = FALSE, smoother.args = list(lty = 2), main = \"Scatter Plot Matrix of activationEnergy\")\n从上图我们可以看出，着火温度是一个单峰曲线，随着水分、挥发分、氧含量的升高而降低，随着固定碳以及发热量升高而升高。活化能也有类似的规律。5 回归分析下面使用lm()函数分别拟合各自的多元回归模型，代码如下：& fit_AE &- lm(activationEnergy ~ moistrue + volatileMatter + fixedCarbon + LHV + oxygen, data = AE_Estimate)\n& summary(fit_AE)\n\nCall:\nlm(formula = activationEnergy ~ moistrue + volatileMatter + fixedCarbon + \n
LHV + oxygen, data = AE_Estimate)\n\nResiduals:\n
\n\nCoefficients:\n
Estimate Std. Error t value Pr(&|t|)
\n(Intercept)
3.535 0.000808 ***\nmoistrue
-1.090 0.280255
\nvolatileMatter
-2.493 0.015536 *
\nfixedCarbon
-1.099 0.276134
1.875 0.065797 .
0.408 0.684715
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nResidual standard error: 33440 on 58 degrees of freedom\nMultiple R-squared:
0.465,\tAdjusted R-squared:
0.4188 \nF-statistic: 10.08 on 5 and 58 DF,
p-value: 5.577e-07\n\n& fit_IT &- lm(ignitionTemperature ~ moistrue + volatileMatter + fixedCarbon + LHV + oxygen, data = IT_Estimate)\n& summary(fit_IT)\n\nCall:\nlm(formula = ignitionTemperature ~ moistrue + volatileMatter + \n
fixedCarbon + LHV + oxygen, data = IT_Estimate)\n\nResiduals:\n
Max \n-57.360
57.860 \n\nCoefficients:\n
Estimate Std. Error t value Pr(&|t|)
\n(Intercept)
& 2e-16 ***\nmoistrue
0.00580 ** \nvolatileMatter
5.9e-05 ***\nfixedCarbon
0.00271 ** \noxygen
\n---\nSignif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n\nResidual standard error: 18.61 on 58 degrees of freedom\nMultiple R-squared:
0.6926,\tAdjusted R-squared:
0.6661 \nF-statistic: 26.13 on 5 and 58 DF,
p-value: 1.03e-13\n采用标准方法对模型进行诊断：& confint(fit_AE)\n
97.5 %\n(Intercept)
9. 66\nmoistrue
\nvolatileMatter -1.
-\nfixedCarbon
27.62129\noxygen
\n& confint(fit_IT)\n
97.5 %\n(Intercept)
429.0.\nmoistrue
-1.\nvolatileMatter
-3.\nfixedCarbon
0.\noxygen
5.\n& par(mfrow=c(2,2))\n& plot(fit_AE)\n& plot(fit_AE,main = \"活化能诊断图\")\n& plot(fit_IT,main = \"着火温度诊断图\")\n我们可以看出总体上数据符合正态性、线性的假设，在活化能的模型当中可能需要适当的加入一些二次项，另外，两个模型都有部分离群点，可能是实验的误差导致。","state":"published","sourceUrl":"","pageCommentsCount":0,"canComment":false,"snapshotUrl":"","slug":,"publishedTime":"T21:52:56+08:00","url":"/p/","title":"配煤动力学特性回归模型的R语言实现","summary":"之前在研究生毕业论文当中用SPSS采用多元线性及非线性模型用于预测配煤动力学特性，本节主要结合R语言提供的函数和工具重现一下。1 几个概念：单煤：就是从一个矿区开采出来的煤炭，一般性质相近；配煤：为了取得特定的特性，由不同种单煤按照不同的比例混…","reviewingCommentsCount":0,"meta":{"previous":null,"next":null},"commentPermission":"anyone","commentsCount":0,"likesCount":0},"next":{"isTitleImageFullScreen":false,"rating":"none","titleImage":"/v2-e3f3897dbdb0ae83df5b16_r.png","links":{"comments":"/api/posts//comments"},"topics":[{"url":"/topic/","id":"","name":"R（编程语言）"}],"adminClosedComment":false,"href":"/api/posts/","excerptTitle":"","author":{"bio":"数据分析/机器学习 修炼中","isFollowing":false,"hash":"20c04cea624a1ad","uid":28,"isOrg":false,"slug":"realphy","isFollowed":false,"description":"人生没有起跑线，终生学习践行者","name":"Datacruiser","profileUrl":"/people/realphy","avatar":{"id":"v2-6f18fa1ab887d9f7a722596","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},"column":{"slug":"datacruiser","name":"数海拾荒"},"content":"功效分析可以帮助在给定置信度的情况下，判断检测到给定效应值时所需的样本量。在给定置信度水平情况下，计算在某样本量内能检测到给定效应值的概率。在开始学习R语言中功效分析工具包之前首先回顾下统计学中假设检验的概念。1 假设检验回顾在统计假设检验中，首先要对总体分布参数设定一个假设（零假设，H0），然后从总体分布中抽样，通过样本计算所得的统计量来对总体参数进行推断。假定零假设为真，如果计算获得观测样本的统计量的概率非常小，便可以拒绝原假设，接受它的对立面（称作备择假设或者研究假设，H1）。在使用取自总体的样本数据来对总体做推断时，通常会有以下四种情况：在研究过程当中，我们一般关注以下四个量：样本大小、显著性水平、功效和效应值。样本大小指的是实验设计中每种条件/组中观测的数目。显著性水平（也称为alpha）由Ⅰ型错误的概率来定义。也可以把它看作发现效应不发生的概率功效通过1减去Ⅱ型错误的概率来定义。我们可以把它看作真实效应发生的概率。效应值指的是在备择或研究假设下效应的量。效应值的表达式依赖于假设检验中使用的统计方法。其中，样本大小和显著性水平可以直接控制，而功效和效应值的影响时间接的。四个量紧密相关，给定其中任意三个量，便可以推算第四个量。这也是下面学习如何使用R语言中pwr包实现功效分析的基础。2 pwr包做功效分析t检验格式：pwr.t.test(n=, d=, sig.level=, power=, type=, alternative=)\n其中：1, n为样本大小2, d为效应值，即标准化的均值之差
3, sig.level表示显著性水平（默认为0.05）4, power为功效水平5, type指检验类型：双样本t检验（“two.sample”）,单样本t检验（\"one.sample\"）或相依样本t检验（\"paired\"）。默认为双样本t检验6, \"alternative\"指统计检验是双侧检验（\"two.sided\"）还是单侧检验（\"less\"或\"greater\"）。默认为双侧检验实例：求样本大小：& library(pwr)\n& pwr.t.test(d=.8,sig.level = .05,power = .9,type = \"two.sample\",alternative = \"two.sided\")\n\n
Two-sample t test power calculation \n\n
n = 33.82554\n
sig.level = 0.05\n
power = 0.9\n
alternative = two.sided\n\nNOTE: n is number in *each* group\n\n& \n求功效水平：& pwr.t.test(n=20,d=.5,sig.level = .01)\n\n
Two-sample t test power calculation \n\n
sig.level = 0.01\n
power = 0.1439551\n
alternative = two.sided\n\nNOTE: n is number in *each* group\n方差分析格式：pwr.anova.test(k=, n=, f=, sig.level=, power=)\n其中：1, k为组的个数2, n为各组中样本大小2, f为效应值，对于单因素方差分析，3, sig.level表示显著性水平（默认为0.05）4, power为功效水平实例：求样本大小：& library(pwr)\n& pwr.anova.test(k = 5, f = .25, sig.level = .05,power = .8)\n\n
Balanced one-way analysis of variance power calculation \n\n
n = 39.1534\n
f = 0.25\n
sig.level = 0.05\n
power = 0.8\n\nNOTE: n is number in each group\n相关性格式：pwr.r.test(n=, r=, sig.level=, power=, alternative=)\n其中：1, n为观测数目2, r为效应值（通过线性相关系数衡量）3, sig.level是显著性水平4, power为功效水平5, \"alternative\"指统计检验是双侧检验（\"two.sided\"）还是单侧检验（\"less\"或\"greater\"）。默认为双侧检验实例：求样本大小：& library(pwr)\n& pwr.r.test(r = .25, sig.level = 0.05, power = .90,alternative = \"greater\")\n\n
approximate correlation power calculation (arctangh transformation) \n\n
n = 133.2803\n
r = 0.25\n
sig.level = 0.05\n
power = 0.9\n
alternative = greater\n\n& \n线性模型格式：pwr.f2.test(u=, v=, f2=, sig.level=, power=)\n其中：1, u为分子自由度，v为分母自由度2,为效应值当评价目的的不同，有两种不同的表达式，如果评价一组预测变量对结果的影响程度时，则：，其中为多重相关性的总体平方值如果要评价的是预测变量对结果的影响超过第二组变量—即协变量—多少时，则：，其中为集合A中变量对总体方差的解释率，为集合A和B中变量对总体方差的解释率3, sig.level是显著性水平4, power为功效水平5, \"alternative\"指统计检验是双侧检验（\"two.sided\"）还是单侧检验（\"less\"或\"greater\"）。默认为双侧检验实例：研究老板的领导风格对员工满意度的影响，是否超过薪水和工作小费对员工满意度的影响。已知领导风格能够解释35%的方差，薪水和消费能够解释约30%的员工满意度的方差。假定显著性水平为0.05，在90%的置信度下。求样本大小：& pwr.f2.test(u=3,f2 = 0.0769,sig.level = 0.05,power = 0.9)\n\n
Multiple regression power calculation \n\n
v = 184.2426\n
f2 = 0.0769\n
sig.level = 0.05\n
power = 0.9\n\n&\n在多元回归中，分母的自由度等于N-k-1（），其中N是总观测数，k是预测变量数。本例中，N-7-1 = 185，则样本大小为193。比例检验格式（各组中n相同）：pwr.2p.test(h=, n=, sig.level=, power=)\n其中：1，h为效应值h的定义如下：，在R语言当中可以用ES.h(p1,p2)函数来计算。2，n为各组相同的样本量格式（各组中n不相同）：pwr.2p2n.test(h=, n1=, n2=, sig.level=, power=)\n实例：药物差异的药物受试者样本量问题。求样本大小：& pwr.2p.test(h=ES.h(.65,.6), sig.level = .05, power = .9,alternative = \"greater\")\n\n
Difference of proportion power calculation for binomial distribution (arcsine transformation) \n\n
h = 0.1033347\n
sig.level = 0.05\n
power = 0.9\n
alternative = greater\n\nNOTE: same sample sizes\n\n& \n卡方检验卡方检验常用于检验两个类型变量之间的关系，典型的零假设是变量之间独立，备择假设是不独立。格式：pwr.chisq.test(w=, N=, df=, sig.level=, power=)\n其中：1，N是总样本大小2，df是自由度3，w是效应值，定义如下：，其中时第i单元格中的概率，时第i单元格中的概率，并且从1到m进行求和，m指的是列联表中单元格的数目。R语言当中ES.w2(P)函数可以计算双因素列联表中的备则假设的效应值。实例：研究对象的双因素列联表如下：则，效应值计算如下：& prob &- matrix(c(.42, .28, .03, .07, .10, .10), byrow = TRUE, nrow=3)\n& ES.w2 (prob)\n[1] 0.1853198\n另外，双因素列联表的自由度为(r-1)(c-1)。& pwr.chisq.test(w=.1853, df = 2, sig.level = .05, power = .9)\n\n
Chi squared power calculation \n\n
w = 0.1853\n
N = 368.5317\n
sig.level = 0.05\n
power = 0.9\n\nNOTE: N is the number of observations\n效应值的选择功效分析当中，难在预期效应值的选择。在行为科学领域，Cohen曾经提出一个基准。当然，Cohen的基准值只是根据许多社科类研究提出的一般性建议，特定问题的研究还是需要具体问题具体分析。另外，我们也可以改变研究参数，并记录其对诸如样本大小和功效等方面的影响。现以五分组的单因素方差分析为例，检测一系列的效应值所需要的样本大小。library(pwr)\nes &- seq(.1, .5, .01)\nnes &- length(es)\nsamsize &- NULL\nfor(i in 1:nes){\n
result &- pwr.anova.test(k = 5, f = es[i], sig.level = .05, power = .9)\n
samsize[i] &- ceiling(result$n)\n}\n\nplot(samsize, es, type = \"l\", lwd = 2, col = \"red\",\n
ylab = \"Effect Size\",\n
xlab = \"Sample Size (per cell)\",\n
main = \"One Way ANOVA with Power = .90 and Alpha = .05\")\n2 绘制功效分析图形其实跟上一节当中绘制的图形类似，用for循环来计算一系列效应值和功效水平下所需的样本量，并以相关性的功效分析为例。#生成相关系数和功效值系列\nlibrary(pwr)\nr &- seq(.1, .5, .01)\nnr &- length(r)\n\np &- seq(.4, .9, .1)\nnp &- length(p)\n\n#获取样本大小\nsamsize &- array(numeric(nr*np), dim=c(nr,np))\nfor (i in 1:np){\n
for (j in 1:nr){\n
result &- pwr.r.test(n = NULL, r = r[j],\n
sig.level = .05, power = p[i],\n
alternative = \"two.sided\")\n
samsize[j,i] &- ceiling(result$n)\n
}\n}\n\n#创建图形\nxrange &- range(r)\nyrange &- round(range(samsize))\ncolors &- rainbow(length(p))\nplot(xrange, yrange, type=\"n\",\n
xlab=\"Correlation Coefficient (r)\",\n
ylab=\"Sample Size (n)\" )\n\n#添加功效曲线\nfor (i in 1:np){\n
lines(r, samsize[,i], type=\"l\", lwd=2, col=colors[i])\n}\n\n#添加网格线\nabline(v=0, h=seq(0,yrange[2],50), lty=2, col=\"grey89\")\nabline(h=0, v=seq(xrange[1],xrange[2],.02), lty=2, col=\"gray89\")\n\n#添加注释\ntitle(\"Sample Size Estimation for Correlation Studies\\n\nSig=0.05 (Two-tailed)\")\nlegend(\"topright\", title=\"Power\", as.character(p),\n
fill=colors)\n4 其他软件包除了pwr包，在研究的规划阶段，R还提供了不少其他有用的软件包可供选择，限于篇幅，原文并没有做过多的叙述和举例说明。其中，后面五个包主要聚焦于基因研究中的功效分析，识别基因与可观测特征的关联性的研究称为全基因组关联研究（GWAS）。此外，MBESS包也包含了可供各种形式功效分析所用的函数，有兴趣的朋友可以自行help了解。5 小结在前面三章，主要探索了各种各样的统计假设检验函数，本章更加关注研究的准备阶段，功效分析不仅可以帮助你判断在给定置信度和效应值的前提下所需的样本量，也能说明在给定样本量时检测到要求效应值的概率。对于限定误报效应显著性的可能性（Ⅰ型错误）和正确检测真实效应（功效）的可能性的平衡，也有了初步的了解。本章主要内容是接受pwr包中各个函数的具体用法，这些函数可以对常见的统计方法进行功效分析，特别是计算样本量。最后也给出了一些专业化的功效分析方法。我们必须记住的是，典型的功效分析是一个交互性的过程。研究者会通过改变样本量、效应值、预期显著性水平和预期功效水平等参数，来观测它们对于其他参数的影响。PS：最近连环出差，而且是日常排得满满当当的差，这篇笔记其实上上周就开始提笔了，直到现在才完成。另外，虽然书已经看到第13章了，但是复习笔记才更新到第十章，也是一个回顾总结的过程，今天因为 同学昨晚在群里的一个问题又回过头去复习了一遍各个检验之间的区别，甚好甚好。","state":"published","sourceUrl":"","pageCommentsCount":0,"canComment":false,"snapshotUrl":"","slug":,"publishedTime":"T21:18:41+08:00","url":"/p/","title":"《R语言实战》第三部分第十章-功效分析学习笔记","summary":"功效分析可以帮助在给定置信度的情况下，判断检测到给定效应值时所需的样本量。在给定置信度水平情况下，计算在某样本量内能检测到给定效应值的概率。在开始学习R语言中功效分析工具包之前首先回顾下统计学中假设检验的概念。1 假设检验回顾在统计假设检验…","reviewingCommentsCount":0,"meta":{"previous":null,"next":null},"commentPermission":"anyone","commentsCount":0,"likesCount":0}},"annotationDetail":null,"commentsCount":0,"likesCount":4,"FULLINFO":true}},"User":{"realphy":{"isFollowed":false,"name":"Datacruiser","headline":"人生没有起跑线，终生学习践行者","avatarUrl":"/v2-6f18fa1ab887d9f7a722596_s.jpg","isFollowing":false,"type":"people","slug":"realphy","bio":"数据分析/机器学习 修炼中","hash":"20c04cea624a1ad","uid":28,"isOrg":false,"description":"人生没有起跑线，终生学习践行者","profileUrl":"/people/realphy","avatar":{"id":"v2-6f18fa1ab887d9f7a722596","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false,"badge":{"identity":null,"bestAnswerer":null}}},"Comment":{},"favlists":{}},"me":{},"global":{},"columns":{"datacruiser":{"following":false,"canManage":false,"href":"/api/columns/datacruiser","name":"数海拾荒","creator":{"slug":"realphy"},"url":"/datacruiser","slug":"datacruiser","avatar":{"id":"v2-87c2b69a215bae3b8858fca","template":"/{id}_{size}.jpg"}}},"columnPosts":{},"postComments":{},"postReviewComments":{"comments":[],"newComments":[],"hasMore":true},"favlistsByUser":{},"favlistRelations":{},"promotions":{},"switches":{"couldAddVideo":false},"draft":{"titleImage":"","titleImageSize":{},"isTitleImageFullScreen":false,"canTitleImageFullScreen":false,"title":"","titleImageUploading":false,"error":"","content":"","draftLoading":false,"globalLoading":false,"pendingVideo":{"resource":null,"error":null}},"drafts":{"draftsList":[]},"config":{"userNotBindPhoneTipString":{}},"recommendPosts":{"articleRecommendations":[],"columnRecommendations":[]},"env":{"isAppView":false,"appViewConfig":{"content_padding_top":128,"content_padding_bottom":56,"content_padding_left":16,"content_padding_right":16,"title_font_size":22,"body_font_size":16,"is_dark_theme":false,"can_auto_load_image":true,"app_info":"OS=iOS"},"isApp":false},"sys":{}}}