R_meandiff_2.Rmd

---
title: "R的假设检验之二"
author: "李峰"
date: “2019/05”
output:
  html_document: default
---




```{r}
library(sm)
```



---


####  1. 单因素方差分析


#####　1.1 F分布密度函数


```{r}
###############F分布
set.seed(12345)
x<-rnorm(1000,0,1)
Ord<-order(x,decreasing=FALSE)
x<-x[Ord]
y<-dnorm(x,0,1)
plot(x,y,xlim=c(-1,5),ylim=c(0,2),type="l",ylab="密度",main="标准正态分布与不同自由度下的F分布密度函数",lwd=1.5)
#######不同自由度的F分布
df1<-c(10,15,30,100)
df2<-c(10,20,25,110)
for(i in 1:4){
 x<-rf(1000,df1[i],df2[i])
 Ord<-order(x,decreasing=FALSE)
 x<-x[Ord]
 y<-df(x,df1[i],df2[i])
 lines(x,y,lty=i+1)
}
legend("topright",title="自由度",c("标准正态分布",paste(df1,df2,sep="-")),lty=1:5)
```



#####　1.2 单因素方差分析示例


```{r}
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
aov(MPG~ModelYear,data=CarData)
OneWay<-aov(MPG~ModelYear,data=CarData)
anova(OneWay)
summary(OneWay)
```


* 关于aov中R的表达式，下图来自**R in Action**
![](http://i2.tiimg.com/611786/02527d8825a33882.jpg)



#####　1.3 单因素方差分析的可视化


```{r}
# install.packages("gplots")
library("gplots")
plotmeans(MPG~ModelYear,data=CarData,p=0.95,use.t=TRUE,xlab="年代车型",ylab="平均MPG",main="不同年代车型MPG总体均值变化折线图(95%置信区间)")

```


**p=0.95,use.t=TRUE**表示95%的置信区间，用t统计量估计总体均值的置信区间。



#####　1.4 单因素方差分析的正态性和方差齐性检验


```{r}
par(mfrow=c(3,5),mar=c(4,4,4,4))
for(i in unique(CarData$ModelYear)){
 T<-subset(CarData,CarData$ModelYear==i)
 qqnorm(T$MPG,main=paste(i,"年车型mpg Q-Q图"),cex=0.7,cex.main=0.9)
 qqline(T$MPG,distribution = qnorm)
}
```



```{r}
library(car)
 qqPlot(lm(MPG ~ ModelYear, data = CarData), simulate = TRUE, 
    main = "车型mpg Q-Q图", labels = FALSE)
```



```{r}
library("lattice")
qqmath(~MPG|ModelYear,data=CarData)
```

* KS检验

```{r}
ks.test(rnorm(100),"pnorm")
ks.test(rnorm(100),"punif")
for(i in unique(CarData$ModelYear)){
 T<-subset(CarData,CarData$ModelYear==i)
 R<-ks.test(T$MPG,"pnorm")
 print(R)
}
```

* 方差齐性检验

```{r}
leveneTest(CarData$MPG,CarData$ModelYear, center=mean)
```


```{r}
bartlett.test(MPG ~ ModelYear, data = CarData)
```



#####　1.5 单因素方差分析的事后检验



```{r}
OneWay<-aov(MPG~ModelYear,data=CarData)
OneWay$coefficients

```




```{r}
TukeyHSD(OneWay,ordered=FALSE,conf.level=0.95)
Result<-TukeyHSD(OneWay,ordered=TRUE,conf.level=0.95)

```


```{r}
LineCol<-vector()
LineCol[Result[[1]][,4]<0.05]<-2
LineCol[Result[[1]][,4]>=0.05]<-1
par(las=2) # 旋转标签
par(mar=c(5,8,4,2)) # 增大左边界位置
plot(Result,cex.axis=0.5,col=LineCol)
```


#####　1.6 单因素方差分析的功效分析


方差分析里效应量的定义：

\begin{equation}\label{eq:1}
  f = \sqrt{\frac{\sum_{i=1}^{k}{\frac{n_i}{n}(\bar{x_i}-\bar{x})^2}}{MSE}},
\end{equation}


```{r}
library("pwr")
pwr.anova.test(k=13,f=0.25,sig.level=0.05,power=0.8)
pwr.anova.test(k=13,f=0.25,sig.level=0.05,n=30)
```


```{r}
ES<-seq(from=0.1,to=0.8,by=0.01)
SampleSize<-matrix(nrow=length(ES),ncol=8)
for(i in 3:10){
 for(j in 1:length(ES)){
  result<-pwr.anova.test(k=i,f=ES[j],sig.level=0.05,power=0.8)
  SampleSize[j,i-2]<-ceiling(result$n)
  }
 }
plot(SampleSize[,1],ES,type="l",ylab="效应量",xlab="样本量（每个水平）",main="单因素方差分析(Alpha=0.05,Power=0.8)")
for(i in 2:8){
 lines(SampleSize[,i],ES,type="l",col=i)
}
legend("topright",title="水平数",paste("k",3:10,sep="="),lty=1,col=1:8)

```