#########################################################
#########################################################
##########                                      #########
##########	R Scripts for MATH 204 Winter 2013  #########
##########                                      #########
#########################################################
#########################################################

#DNA data Chisq Test
DNA.data<-c(38514,24631,25685,38249)
DNA.fitted<-rep(1/4,4)*sum(DNA.data)

(Xsq<-sum((DNA.data-DNA.fitted)^2/DNA.fitted))
(pval<-1-pchisq(Xsq,df=3))

#DNA data Chisq Test: alternative null
DNA.data<-c(38514,24631,25685,38249)
DNA.fitted<-c(0.3,0.2,0.2,0.3)*sum(DNA.data)

(Xsq<-sum((DNA.data-DNA.fitted)^2/DNA.fitted))
(pval<-1-pchisq(Xsq,df=3))


#########################################################
#Hair & Eye colour

xtabs(count~eye+hair,data=HairEyes)

chisq.test(xtabs(count~eye+hair,data=HairEyes))

#Direct calculation
HairEyes.table<-xtabs(count~eye+hair,data=HairEyes)
HE.rowtots<-apply(HairEyes.table,1,sum)
HE.coltots<-apply(HairEyes.table,2,sum)

#Form the fitted values from the row and column totals
HE.fitted<-outer(HE.rowtots,HE.coltots,"*")/sum(HairEyes.table)

Xsq<-sum((HairEyes.table-HE.fitted)^2/HE.fitted)

(pval<-1-pchisq(Xsq,df=9))

#########################################################
#########################################################
#Non-parametric tests

#Mann-Whitney-Wilcoxon test

Sample1<-c(0.31,0.48,1.02,3.11)
Sample2<-c(0.16,0.20,1.97)
Y<-c(Sample1,Sample2)
Samp<-rep(c(1,2),c(4,3))

rY<-rank(Y)
(R1<-sum(rY[Samp==1]))
(R2<-sum(rY[Samp==2]))

wilcox.test(Y~Samp)
wilcox.test(Y~Samp,alternative="less")
wilcox.test(Y~Samp,alternative="greater")

require(coin)
wilcox_test(Y~as.factor(Samp)) #Asymptotic approximation
wilcox_test(Y~as.factor(Samp),distribution="exact") #Exact null
wilcox_test(Y~as.factor(Samp),distribution="approximate") #Approximate null using permutation Monte Carlo

######################################################################
#Sign test

eta<-9
binom.test(sum(Sign1$h2o>eta),nrow(Sign1),p=0.5,alternative="greater")
binom.test(sum(Sign1$h2o>eta),nrow(Sign1),p=0.5,alternative="less")
binom.test(sum(Sign1$h2o>eta),nrow(Sign1),p=0.5,alternative="two.sided")

#Test statistic
S<-sum(Sign1$h2o>9)
#First p-value
1-pbinom(S-1,nrow(Sign1),0.5)
#Second p-value
pbinom(S,nrow(Sign1),0.5)
#Third p-value
2*(1-pbinom(max(S,nrow(Sign1)-S)-1,nrow(Sign1),0.5))

#####################################################################
eta<-9
binom.test(sum(Sign2$h2o>eta),nrow(Sign2),p=0.5,alternative="greater")
binom.test(sum(Sign2$h2o>eta),nrow(Sign2),p=0.5,alternative="less")
binom.test(sum(Sign2$h2o>eta),nrow(Sign2),p=0.5,alternative="two.sided")

#Test statistic
S<-sum(Sign2$h2o>9)
#First p-value
1-pbinom(S-1,nrow(Sign2),0.5)
#Second p-value
pbinom(S,nrow(Sign2),0.5)
#Third p-value
2*(1-pbinom(max(S,nrow(Sign2)-S)-1,nrow(Sign2),0.5))


#####################################################################

#Two sample Wilcoxon:Birthweight

attach(Birthweight)
(R1<-sum(rank(BW)[gp==1]))
(R2<-sum(rank(BW)[gp==2]))
n1<-sum(gp==1)
n2<-sum(gp==2)
detach(Birthweight)

W<-R2
U<-W-0.5*n2*(n2+1)
Z<-(U-0.5*n1*n2)/sqrt(n1*n2*(n1+n2+1)/12)

wilcox.test(BW~Group,data=Birthweight)
wilcox.test(BW~Group,data=Birthweight,alternative="greater")
wilcox.test(BW~Group,data=Birthweight,alternative="less")

wilcox_test(BW~Group,data=Birthweight,dist="exact",ties.method="mid-ranks")
wilcox_test(BW~Group,data=Birthweight,dist="exact",alternative="greater")
wilcox_test(BW~Group,data=Birthweight,dist="exact",alternative="less")

#####################################################################

#Two sample Wilcoxon:Treadmill

attach(Treadmill)
(R1<-sum(rank(Time)[Group==1]))
(R2<-sum(rank(Time)[Group==2]))
n1<-sum(Group==1)
n2<-sum(Group==2)
detach(Treadmill)

W<-R2
U<-W-0.5*n2*(n2+1)
Z<-(U-0.5*n1*n2)/sqrt(n1*n2*(n1+n2+1)/12)

wilcox.test(Time~Group,data=Treadmill)
wilcox.test(Time~Group,data=Treadmill,alternative="greater")
wilcox.test(Time~Group,data=Treadmill,alternative="less")

wilcox_test(Time~as.factor(Group),data=Treadmill,dist="exact",
		ties.method="mid-ranks")
wilcox_test(Time~as.factor(Group),data=Treadmill,dist="exact",
		alternative="greater")
wilcox_test(Time~as.factor(Group),data=Treadmill,dist="exact",
		alternative="less")


#####################################################################
#Paired test
#Haemodialysis
wilcox.test(Haemo$after,Haemo$before,paired=T)
wilcox.test(Haemo$after,Haemo$before,paired=T,alternative="greater")
wilcox.test(Haemo$after,Haemo$before,paired=T,alternative="less")

#Direct calculation
xi<-Haemo$after-Haemo$before
xi<-xi[xi!=0]
si<-abs(xi)
n<-length(si)

Tplus<-sum(rank(si)[xi > 0])
Tminus<-sum(rank(si)[xi < 0])

(Z<-(Tplus-0.25*n*(n+1))/sqrt(n*(n + 1)*(2*n + 1)/24))

pnorm(Z)

#PEFRAsthma
wilcox.test(PEFRAsthma$after,PEFRAsthma$before,paired=T)
wilcox.test(PEFRAsthma$after,PEFRAsthma$before,paired=T,alternative="greater")
wilcox.test(PEFRAsthma$after,PEFRAsthma$before,paired=T,alternative="less")

#Direct calculation
xi<-PEFRAsthma$after-PEFRAsthma$before
xi<-xi[xi!=0]
si<-abs(xi)
n<-length(si)

Tplus<-sum(rank(si)[xi > 0])
Tminus<-sum(rank(si)[xi < 0])

(Z<-(Tplus-0.25*n*(n+1))/sqrt(n*(n + 1)*(2*n + 1)/24))

pnorm(Z)



#####################################################################
#Mucociliary data
#Mean ranks
xtabs(rank(y)~group,data=Muco)/xtabs(~group,data=Muco)
kruskal.test(y~group,data=Muco)

#MemoryTask
#Mean ranks
xtabs(rank(Memory)~Task,data=Memory)/xtabs(~Task,data=Memory)

kruskal.test(Memory~Task,data=Memory)

#####################################################################
#Skin potential data (Hypnosis)

friedman.test(SkinResponse~EmoCode | Subject,data=Hypnosis)

#Soil
friedman.test(sulphur~solvent | soil,data=Soil)

#####################################################################
#Rank correlation:DZT data

plot(DZT$Lat,DZT$DZT,pch=20,xlab="Latitude (degrees)",ylab="DZT Rate (per 1000, standardized)",xlim=range(37,65),ylim=range(5,12.5))
DZT.ylab<-DZT$DZT+0.25
DZT.ylab[c(11,13,18)]<-DZT.ylab[c(11,13,18)]-0.45
text(DZT$Lat,DZT.ylab,DZT$Country,cex=0.6)

cor(DZT$Lat,DZT$DZTRate,method="spearman")
cor(DZT$Lat,DZT$DZTRate,method="pearson")

#####################################################################
#####################################################################
#Permutation/small sample distribution of the Wilcoxon statistic

n1<-7
n2<-3
n<-n1+n2

Rvec<-c(1:n)
require(combinat)
Cmat<-combn(n,n2)
Scores<-apply(Cmat,2,sum)
Prob<-table(Scores)/length(Scores)
CProb<-cumsum(Prob)