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")

　　

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)

Chi squared power calculation

w = 0.1853198
N = 368.4528
df = 2
sig.level = 0.05
power = 0.9

NOTE: N is the number of observations

需要對368個顧客做調查才有90%的目標實現既定目標。