R的两均值比较检验（参数检验）

1.由于抽样的随机性，样本均值在不同总体上的差距很可能是由抽样误差引起的，而这种差距不被认为具有统计上的显著性。

2.反之，若分析发现样本均值在不同总体上差距较大，但不是由抽样误差引起的，则数值型变量在不同总体上的分布参数存在显著差异。

检验两个样本上的均值差是否统计显著的方法：参数检验&非参检验，步骤：

h0&h1
构造检验统计。该检验统计量在原假设成立条件下，服从某个已知的理论分布，这称为抽样分布。
依据样本数据计算在原假设成立的条件下，检验统计量的观测值与概率P值。检验统计量反映了观测值与原假设之间的差距，p反映了在原假设成立条件下检验统计量取当前观测值或更极端的可能性。
指定显著新水平α，原假设成立却拒绝的概率
power:1-β,p(H0|H1)

1.两独立样本的均值检验

1.1.概述

适用数据：

观测样本来自总体中的两个独立样本，抽样个过程中互不干扰

检验目标：

量样本均值是否具有统计上的显著性。不具有显著性：均值差是由抽样误差导致的。

理论依据：

1.2抽样自举：

###############利用bootstrap模拟独立样本均值差的抽样分布
par(mfrow=c(2,1),mar=c(4,4,4,4))
set.seed(12345)

#总体方差相等
Pop1<-rnorm(10000,mean=2,sd=2)   
Pop2<-rnorm(10000,mean=10,sd=2)
Diff<-vector()
Sdx1<-vector()
Sdx2<-vector()
#重复M次
for(i in 1:2000){
 x1<-sample(Pop1,size=100,replace=TRUE)#随机选出100个
 x2<-sample(Pop2,size=120,replace=TRUE)
 Diff<-c(Diff,(mean(x1)-mean(x2)))
 Sdx1<-c(Sdx1,sd(x1))
 Sdx2<-c(Sdx2,sd(x2))
}
plot(density(Diff),xlab="mean(x1)-mean(x2)",ylab="Density",main="均值差的抽样分布(等方差)",cex.main=0.7,cex.lab=0.7) 
points(mean(Diff),sd(Diff),pch=1,col=1)
S1<-mean(Sdx1)
S2<-mean(Sdx2)
Sp<-((100-1)*S1^2+(120-1)*S2^2)/(100+120-2)
#理论上的均值与方差：红三角
points((2-10),sqrt(Sp/100+Sp/120),pch=2,col=2)

###两方差不等
set.seed(12345)
Pop1<-rnorm(10000,mean=2,sd=2)    
Pop2<-rnorm(10000,mean=10,sd=4)
Diff<-vector()
Sdx1<-vector()
Sdx2<-vector()
for(i in 1:2000){
 x1<-sample(Pop1,size=100,replace=TRUE)
 x2<-sample(Pop2,size=120,replace=TRUE)
 Diff<-c(Diff,(mean(x1)-mean(x2)))
 Sdx1<-c(Sdx1,sd(x1))
 Sdx2<-c(Sdx2,sd(x2))
 }
plot(density(Diff),xlab="mean(x1)-mean(x2)",ylab="Density",main="均值差的抽样分布(不等方差)",cex.main=0.7,cex.lab=0.7) 
points(mean(Diff),sd(Diff),pch=1,col=1)
S1<-mean(Sdx1)
S2<-mean(Sdx2)
points((2-10),sqrt(S1^2/100+S2^2/120),pch=2,col=2)

1.3检验单原假设与检验统计量

实现两独立样本均值检验单R程序：R.test

#############独立样本均值检验示例
Forest<-read.table(file="ForestData.txt",header=TRUE,sep="	")
#设置因子
Forest$month<-factor(Forest$month,levels=c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec"))
Tmp<-subset(Forest,Forest$month=="jan" | Forest$month=="aug")
#方差不同
t.test(temp~month,data=Tmp,paired=FALSE,var.equal=TRUE)

Two Sample t-test

data: temp by month
t = -4.8063, df = 184, p-value = 3.184e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-23.106033 -9.657011
sample estimates:
mean in group jan mean in group aug
5.25000 21.63152

#方差相同

t.test(temp~month,data=Tmp,paired=FALSE,var.equal=FALSE)

Welch Two Sample t-test

data: temp by month
t = -45.771, df = 177.49, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-17.08782 -15.67522
sample estimates:
mean in group jan mean in group aug
5.25000 21.63152

　　T值都大，P<α=0.05拒绝原假设，两个月平均气温相差很大。

1.4方差是否相等

采用上述那个分析结论，取决于量总体方差是否相等。通常采用F-test，也可采用更为稳健但不依赖总体分布具体形式的Levene's方差同质检验

检验1、8月份温度总体方差是否相等

library("car")
leveneTest(Tmp$temp,Tmp$month, center=mean)

Levene's Test for Homogeneity of Variance (center = mean)
Df F value Pr(>F)
group 1 2.6773 0.1035
184

　　p>α=0.05,不能拒绝原假设，方差齐性，选择第一种检验结论

2.两配对样本的均值检验

2.1概述

2.1.1使用数据：

观测样本来自总体的两个配对样本，表象为两个样本样本量一样，且两样本的观测具有一一对应关系，可视为观测样本的“前后”或多“侧面”的数据。

例如，某班学生各科测试中，语文和数学的测试成绩。

2.1.2检验目标：

即检验数学成绩总体分布于语文成绩总体分布（均值）是否存在显著差异。

2.1.3理论依据：

##################利用bootstrap模拟样本均值的抽样分布

set.seed(12345)
Pop<-rnorm(100000,mean=4,sd=2)  #正态总体，均值为4，标准差为2
#从总体样本中随机抽取2000个大小为1000的样本，然后测试样本的均值分布
MeanX<-vector()
for(i in 1:2000){
 x<-sample(Pop,size=1000,replace=TRUE)
 MeanX<-c(MeanX,mean(x))
}
plot(density(MeanX),xlab="mean(x)",ylab="Density",main="样本均值的抽样分布",cex.main=0.7,cex.lab=0.7)
points(mean(MeanX),sd(MeanX),pch=1,col=1)
points(4,sqrt(2^2/1000),pch=2,col=2)

2.1.4检验的原假设和检验统计量：

R函数t.test

检验语文成绩与数学成绩是否具有显著差异：

##############配对样本均值检验示例
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")

ReportCard<-na.omit(ReportCard)
t.test(ReportCard$chi,ReportCard$math,paired=TRUE)

Paired t-test

data: ReportCard$chi and ReportCard$math
t = 11.712, df = 57, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
18.48871 26.11474
sample estimates:
mean of the differences
22.30172

存在显著差异

t.test也可以用来进行单样本均值检验（结果一样）

###############单样本的均值检验示例
Diff<-ReportCard$chi-ReportCard$math
t.test(Diff,mu=0)

3.样本均值检验的功效分析

3.1概述

功效-power：Type2 error(原假设为错却接受了)，若原假设为错并且拒绝，则作者这一正确决策的概率为1-β,统计功效。

首先显著性水平：type1 error是影响统计功效的重要因素之一。alpha小，1-beta也小
其次，若事实上两样本均值差异大，Power也高; 若增大样本量会导致均值差的样本分布的方差减小，是的检验统计量t的观测值增大，更容易拒绝两总体均值无差异这个错误的原假设，Power大，所以样本量大小也影响。
由于样本均值是一个绝对量，会受数据计量单位和数量级的影响。所以找到一个可反映两总体分布重叠的相对指标更有意义--effect size(ES)
单侧检验相对于双边检验更容易拒绝原假设。检验方向也影响。

效应量太小，意味着处理即使达到了显著水平，也缺乏实用价值。

常见的几种ES：

a) 两个平均数间的标准差异；

b) 分组自变量与个体因变量分数间的相关--相关效应大小。

c) 方差分析中处理效应的效应大小

Cohen( 1988) 将ES 定义为“总体中存在某种现象的程度”; 具体到NHST 体系中，ES 即“虚无假设H0错误的程度”［9］。这种错误程度可形象理解为虚无假设H0和备择假设H1所代表的两抽样分布分离程度或面积重叠程度。如图1 所示，ES 越大，H0偏离H1而犯错误的程度越明显，两分布的分离程度越高，重叠面积越小，反之亦然.

3.2理论基础

3.3 R程序

pwr.t.test / per.t2n.test

library("pwr")
pwr.t2n.test(n1=2,n2=184,d=4.8,sig.level=0.05,alternative="two.sided")
pwr.t.test(n=58,sig.level=0.05,power=0.8,type="paired",alternative="two.sided")

ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
Tmp<-ReportCard[complete.cases(ReportCard),]
cor.test(Tmp[,5],Tmp[,7],alternative="two.side",method="pearson")
library("pwr")
pwr.r.test(r=0.75,sig.level=0.05,n=58,alternative="two.sided")

Paired t test power calculation

n = 58
d = 0.3742143
sig.level = 0.05
power = 0.8
alternative = two.sided

NOTE: n is number of *pairs*

4.其他功效检验

4.1相关系数检验

pwr.r.test

##############相关系数检验的功效分析
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
Tmp<-ReportCard[complete.cases(ReportCard),]#complete.cases 和 na.omit去除有空值的行
cor.test(Tmp[,5],Tmp[,7],alternative="two.side",method="pearson")
library("pwr")
pwr.r.test(r=0.75,sig.level=0.05,n=58,alternative="two.sided")

Pearson's product-moment correlation

data: Tmp[, 5] and Tmp[, 7]
t = 8.5775, df = 56, p-value = 8.753e-12
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6149204 0.8469769
sample estimates:
cor
0.7535317

approximate correlation power calculation (arctangh transformation)

n = 58
r = 0.75
sig.level = 0.05
power = 0.9999999
alternative = two.sided

说明在显著性水平0.05，样本量58，样本相关系数0.75的条件下，做出拒绝相关的正确决策的概率为0.99

2.列联表卡放检验的功效分析

##############列联表卡方检验的功效分析
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
Tmp<-ReportCard[complete.cases(ReportCard),]
(CrossTable<-table(Tmp[,c(2,12)]))
(ResChisq<-chisq.test(CrossTable,correct=FALSE))
library("pwr")
pwr.chisq.test(sig.level=0.05,N=58,power=0.9,df=3)

avScore
sex B C D E
F 2 13 10 3
M 2 11 12 5
> (ResChisq<-chisq.test(CrossTable,correct=FALSE))

Pearson's Chi-squared test

data: CrossTable
X-squared = 0.78045, df = 3, p-value = 0.8541

Chi squared power calculation

w = 0.4943029
N = 58
df = 3
sig.level = 0.05
power = 0.9

NOTE: N is the number of observations

延伸：

计算效应量：ES.w2()

1.计算未来顾客年龄与消费行为的关系（phi效应量）应出现多大变化才能打破原来的格局

####################计算效应量
(prob<-matrix(c(0.42,0.28,0.03,0.07,0.10,0.10),nrow=3,ncol=2,byrow=TRUE))
ES.w2(prob)

> (prob<-matrix(c(0.42,0.28,0.03,0.07,0.10,0.10),nrow=3,ncol=2,byrow=TRUE))
[,1] [,2]
[1,] 0.42 0.28
[2,] 0.03 0.07
[3,] 0.10 0.10
> ES.w2(prob)
[1] 0.1853198

2.计算样本量

pwr.chisq.test(w=ES.w2(prob),df=(3-1)*(2-1),sig.level=0.05,power=0.9)

> pwr.chisq.test(w=ES.w2(prob),df=(3-1)*(2-1),sig.level=0.05,power=0.9)