R软件学习笔记-6（各种检验）

a<-rnorm(2000,0,1)   ###########产生2000个随机数
plot(a[order(a)],dnorm(a[order(a)]),type="l",main="抽样分布",xlab="统计量",ylab="概率密度")          ######第一个参数是X值，第二个是Y值


###############利用bootstrap模拟 独立样本均值差的抽样分布
par(mfrow=c(2,1),mar=c(4,4,4,4))         #####2*1的布局，每个布局的4个边界均空出4个单位
set.seed(12345)          #####计算机每次产生的随机数并不是固定的，但是通过set.seed()后获得的随机数都是固定的
Pop1<-rnorm(10000,mean=2,sd=2)   ###两总体方差相等
Pop2<-rnorm(10000,mean=10,sd=2)
Diff<-vector()          #一开始赋值为空的向量
Sdx1<-vector()
Sdx2<-vector()
for(i in 1:2000){
 x1<-sample(Pop1,size=100,replace=TRUE)            ##########    replace是否放回抽样，Pop1是要从中抽样的样本总体，size表示抽取多少次
 x2<-sample(Pop2,size=120,replace=TRUE)
 Diff<-c(Diff,(mean(x1)-mean(x2)))        //这是一种递归的做法，通过这样把diff填充完整
 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)        ##在原来的图形上添加点，pch： Plot CHaracter，点的形状，如0为空心正方形，1为空心圆，2为空心三角形等

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)



#############独立样本均值检验示例
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)
t.test(temp~month,data=Tmp,paired=FALSE,var.equal=FALSE)

################levene's方差齐性检验
library("car")
leveneTest(Tmp$temp,Tmp$month, center=mean)


##################利用bootstrap模拟样本均值的抽样分布
set.seed(12345)
Pop<-rnorm(100000,mean=4,sd=2)  #正态总体，均值为4，标准差为2
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)

##############配对样本均值检验示例
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
ReportCard<-na.omit(ReportCard)
t.test(ReportCard$chi,ReportCard$math,paired=TRUE)

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


setwd("C:\\Users\\Administrator\\Desktop\\统计软件\\数据与程序")
################t检验的功效分析
install.packages("pwr")
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")

##############列联表卡方检验的功效分析
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)

####################计算效应量
(prob<-matrix(c(0.42,0.28,0.03,0.07,0.10,0.10),nrow=3,ncol=2,byrow=TRUE))
ES.w2(prob)
pwr.chisq.test(w=ES.w2(prob),df=(3-1)*(2-1),sig.level=0.05,power=0.9)


##################独立样本的曼-惠特尼U检验
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")
wilcox.test(temp~month,data=Tmp)
##################独立样本的K-S检验
x1<-subset(Forest,Forest$month=="jan")
x2<-subset(Forest,Forest$month=="aug")
ks.test(x1$temp,x2$temp)


###############配对样本的Wilcoxon符号秩检验
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
ReportCard<-na.omit(ReportCard)
wilcox.test(ReportCard$chi,ReportCard$math,paired=TRUE)

sum(outer(ReportCard$chi,ReportCard$math,"-")<0)
sum(outer(ReportCard$math,ReportCard$chi,"-")<0)


################置换检验
install.packages("coin")
library("coin")

x1<-c(15,13,11,14,12,10)
x2<-c(1,1,1,2,2,2)
x<-data.frame(x1=x1,x2=x2)
oneway_test(x1~as.factor(x2),data=x,distribution="exact")

#############独立样本均值的置换检验示例
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)
Tmp$month<-as.vector(Tmp$month)
Tmp$month<-as.factor(Tmp$month)
oneway_test(temp~month,data=Tmp,distribution="exact")
oneway_test(temp~month,data=Tmp,distribution="asymptotic")
oneway_test(temp~month,data=Tmp,distribution=approximate(B=1000))


###############spearman等级相关系数置换检验
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
Tmp<-ReportCard[complete.cases(ReportCard),]
cor.test(Tmp[,5],Tmp[,7],alternative="two.side",method="spearman")
set.seed(12345)
spearman_test(math~phy,data=Tmp,distribution=approximate(B=1000))

###############卡方置换检验
Tmp<-ReportCard[complete.cases(ReportCard),]
CrossTable<-table(Tmp[,c(2,12)])  #编制性别和平均分等级的列联表
chisq.test(CrossTable,correct=FALSE)
chisq_test(sex~avScore,data=Tmp,distribution="asymptotic")
set.seed(12345)
chisq_test(sex~avScore,data=Tmp,distribution=approximate(B=1000))

##############配对总体分布差的置换检验
ReportCard<-read.table(file="ReportCard.txt",header=TRUE,sep=" ")
ReportCard<-na.omit(ReportCard)
wilcox.test(ReportCard$chi,ReportCard$math,paired=TRUE)
wilcoxsign_test(chi~math,data=ReportCard,distribution="asymptotic")

###############两样本均值差的自举法检验

DiffMean<-function(DataSet,indices){
 ReSample<-DataSet[indices,]
 diff<-tapply(ReSample[,1],INDEX=as.factor(ReSample[,2]),FUN=mean)
 return(diff[1]-diff[2])
}

install.packages("boot")
library("boot")
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")
Tmp<-cbind(Tmp$temp,Tmp$month)
set.seed(12345)
BootObject<-boot(data=Tmp,statistic=DiffMean,R=20)
BootObject$t0
mean(BootObject$t,na.rm=TRUE)
print(BootObject)
plot(BootObject)
boot.ci(BootObject,conf=0.95,type=c("norm","perc"))