题目
(统计学-基于R-第五版-第九章习题9.1)
一家牛奶公司有4台机器装填牛奶,每桶的容量为4L。下面是从4台机器中抽取的装填量样本数据:
- 检验机器对装填量是否有显著影响(α=0.01)。
- 分析效应量。
- 采用HSD方法比较哪些机器的装填量之间存在显著差异。
- 检验装填量是否满足正态性和方差齐性。
实现代码和结果
1、检验机器对装填量是否有显著影响(α=0.01)。
只考虑一个因子对观测数据产生影响的方差分析为单因子方差分析,本题只考虑机器对牛奶装填量的影响。
首先读入装填量样本数据csv文件,将表格短格式转换为长格式,如下图所示:
使用melt函数。使用aov函数拟合方差分析模型,summary总结输出方差分析模型结果。
设机器对牛奶装填量的影响效应分别为α₁(机器1),α₂(机器2),α₃(机器3),α₄(机器4),提出假设:
H₀:α₁=α₂=α₃=α₄=0(影响不显著)
H₁:α₁、α₂、α₃、α₄至少有一个不等于0(影响显著)。
由于检验值P值(Pr(>F))=0.0135大于α=0.01,接受H₀,所以机器对牛奶装填量影响不显著。
实现代码和输出结果如下:
example7<-read.csv("C:/Users/lenovo/OneDrive/大二/大数据统计学实验/实验表格/装填量样本数据.csv")
#检验机器对装填量是否有显著影响(α=0.01)
library(reshape2)
tab<-melt(example7,variable.name="机器",value.name ="装填量") #短格式转换为长格式
#拟合单因子方差分析模型
model<-aov(装填量~机器,data=tab)
summary(model)
2、分析效应量
效应量用于分析因子(自变量)对观测结果(因变量)的影响程度。
使用EtaSq函数计算效应量。输出结果表示装填量取值的总误差中被机器解释的比例为40.7%,
实现代码及计算结果如下:
#分析效应量
library(DescTools)
model<-aov(装填量~机器,data=tab)
EtaSq(model)
3、采用HSD方法比较哪些机器的装填量之间存在显著差异。
HSD方法用于多重比较,全称是“真实差异显著”。
使用TukeyHSD函数计算P值,P值计算结果表示机器3-机器1之间差异显著,其他差异均不显著。如配对差值置信区间比较图所示,机器3-机器1的置信区间不包含0,表示二者差异显著,其他的都包含0,表示差异不显著。
代码实现和分析结果如下:
#采用HSD方法比较哪些机器的装填量之间存在显著差异
TukeyHSD(model,alpha=0.01)
plot(TukeyHSD(model)) #配对差值置信区间比较图
4、检验装填量是否满足正态性和方差齐性。
检验正态性有两种方法,一个是Q-Q图,一个是Shapiro-Wilk检验。
Q-Q图是状态概率图,可以检验正态性,由图5.2可知,分位点都大致落在直线上及周围,可认为近似服从正态分布。
Shapiro-Wilk中W为检验统计量,最大值为1,统计量越大表明数据越符合正态分布,P值小于指定显著性水平时表示总体不符合正态分布。由图5.3可知,四个机器的W值都接近1,P值也大于0.01,因此可认为符合正态分布。
实现代码和分析结果如下:
#检验装填量是否满足正态性
#用Q-Q图检验正态性
attach(tab)
#总体正态Q-Q图
qqnorm(装填量,xlab="理论分位数",ylab="样本分位数",main="正态Q-Q图")
qqline(装填量,col="red",lwd=2)
par(fig=c(.075,.5,.55,.98),new=TRUE)
hist(装填量,xlab="装填量",ylab="",freq=FALSE,col="lightblue",cex.axis=0.1,cex.lab=0.1,main="")
lines(density(装填量),col="red",lwd=1)
box(col="gray")
#用Shapiro-Wilk检验正态性
shapiro.test(装填量[机器=="机器1"])
shapiro.test(装填量[机器=="机器2"])
shapiro.test(装填量[机器=="机器3"])
shapiro.test(装填量[机器=="机器4"])
检验方差齐性有两种方法,一个是残差图,一个是Levene检验。
残差图的横坐标是模型的拟合值,纵坐标是预测残差,如果拟合值和残差的散点随机分布在一个水平带内且离散程度基本相同,表明满足方差齐性假定。由图5.4可知本题满足方差齐性。
Levene检验可用于正态总体和非正态总体,使用LeveneTest函数计算,由图5.5可知P=0.8189>0.01,满足方差齐性。
实现代码和分析结果如下:
#用残差图检验方差齐性
model<-aov(装填量~机器,data=tab)
plot(model,which=1,main="残差图")
#用Levene检验方差齐性
library(car)
leveneTest(装填量~机器,center=median,data=tab)