### 1.5.1 ###

 setwd("C:/R springer/R CODE 1")
 
  set.seed(1)
  n1<-3
  n2<-4
  n<-n1+n2
  y1<-rnorm(n1)
  y2<-rnorm(n2,mean=1)
  y<-c(y1,y2)
  label<-rep(c(1,2),c(n1,n2))
  y1 ; y2 ; label


### Exact distribution

  C<-choose(n,n1)

  library(combinat)
  index<-seq(1,n)
  index

  pi1<-matrix(unlist(combn(n,n1)),nrow=C,byrow=TRUE)
  pi2<-t(apply(pi1,1,function(x){index[-x]}))
  pi<-cbind(pi1,pi2)
  pi

 
 y.perm<-matrix(y[pi],nrow=C,byrow=F)
 round(y.perm,digits=5)


  T<-apply(y.perm,1,function(x){mean(x[c(1:n1)]) -mean(x[-c(1:n1)])})
  T<-array(T,dim=c(C,1)) 
  T

 sigma0<-sqrt(sum(y^2)/n-mean(y)^2)

 sigma0^2*n*(1/n1+1/n2)/(n-1)
 mean(T)
 sum(T^2)/C-mean(T)^2

## two-sided alternative p-value

  T.star<-T^2
  p<-(sum(T.star >= T.star[1,]))/C
  p

## one-sided alternative p-value

 p<-sum(T<=T[1,])/C
 p



### Monte Carlo distribution

   set.seed(11)
   B<-1000
   T<-array(0,dim=c((B+1),1))
   T[1,]<-mean(y[1:n1])-mean(y[(n1+1):n])
   for(bb in 2:(B+1)){
   y.perm<-sample(y)
   T[bb,]<-mean(y.perm[1:n1])-mean(y.perm[(n1+1):n])
   }
   T.star<-T^2
   p<-sum(T.star[-1,]>=T.star[1,])/B
   p



################# 1.5.2 ##############################



  library(combinat)
  library(mvtnorm)
  n<-4
  q<-2
  rho <- 0.5
  I<-diag(q)
  J<-array(1,dim=c(q,q))
  S<-rho*J+(1-rho)*I
  S

  mu = c(0,1)
  set.seed(101)
  y<-rmvnorm(n,mean=mu,sigma=S)
  y

### Monte Carlo distribution

# C <- 500 ; 
# sgn = matrix(rbinom(C*n,1,.5),ncol=n)
# sgn = rbind(rep(1,n),sgn)

# T.star<-array(0,dim=c((C+1),2))
# T.star[,1]<-sgn%*%y[,1]
# T.star[,2]<-sgn%*%y[,2]

# T.star[1,]			## observed values of the partial test statistics

# partial.p<-apply(T.star,2,function(x){mean(x[-1]>=x[1])})	# partial p-values
# partial.p				


### exact distribution

  C<-2^n
  z<-rep(2,n)
  sgn<- hcube(z)%/%2
  sgn<-apply(sgn,2,function(x){ifelse(x==0,1,-1)})
  sgn


  T.star<-array(0,dim=c(C,2))
  T.star[,1]<-sgn%*%y[,1]
  T.star[,2]<-sgn%*%y[,2]
  T.star


  apply(T.star,2,mean)
  apply(T.star,2,var)*(C-1)/C
  apply(y,2,function(x){sum(x^2)})



#### nonparametric combination (direct combination)


 psi<-apply(T.star,1,sum)
 T.ind<-expand.grid(T.star[,1],T.star[,2])

## bivariate null distribution representation of the partial test statistics [T1*,T2*]
## white/black dots: dependent/independent permutations within variables



 plot(T.ind,pch=16,xlim=c(-3,3),xlab=expression(T[1]^{"*"}),ylab=expression(T[2]^{"*"}))
 points(T.star,col="white",pch=16)
 points(T.star,col="black")
 text(T.star[1,1],5.3,expression(paste("[",T[1],",",T[2],"]")))

 psi[1]				## observed value of the Direct global test statistic

 p.glob = mean(psi>=psi[1])
 p.glob				## global p-value


## observed rejection region (i.e. rejection region at a significance level equal to p.glob)

 tx<-seq(-3,3,length.out=100)
 lines(tx,psi[1]-tx,type="l",col="grey")
 for(i in 1:100){
 lines(c(tx[i],tx[i]),c(psi[1]-tx[i],10),type="l",col="grey")
 }

  partial.p<-(C+1-apply(T.star,2,rank))/C
  partial.p



### bivariate null distribution representation of the partial p-values [p1*,p2*]

 P.ind<-expand.grid(partial.p[,1],partial.p[,2])
 plot(P.ind,pch=16,xlab=expression(p[1]^{"*"}),ylab=expression(p[2]^{"*"}),xlim=c(0,1),ylim=c(0,1))
 points(partial.p,col="white",pch=16)
 points(partial.p,col="black")
 text(partial.p[1,1],.1,expression(paste("[",p[1],",",p[2],"]")))

### Fisher's combining function:

psi = apply(partial.p,1,function(x){-2*log(prod(x))})


 psi[1]				## observed value of Fisher's global test statistic

p.glob = mean(psi>=psi[1])
p.glob				## global p-value


## observed rejection region (i.e. rejection region at a significance level equal to p.glob): area below the grey curve

 tx<-seq(0,1,length.out=100)
 ty<-exp(-.5*(psi[1]+2*log(tx)))	
 lines(tx,ty,type="l",col="grey")

 lines(tx,ty,type="l",col="grey")
 for(i in 1:100){
 lines(c(tx[i],tx[i]),c(0,ty[i]),type="l",col="grey")
 }



##################################
### Liptak's combining function:

psi = apply(partial.p,1,function(x){sum(qnorm(1-x))})

 psi[1]				## observed value of Liptak's global test statistic

p.glob = mean(psi>=psi[1])
p.glob				## global p-value


## observed rejection region (i.e. rejection region at a significance level equal to p.glob): area below the black curve

 tx<-seq(1/100,1,length.out=100)
 ty<-1-pnorm(psi[1]-qnorm(1-tx))
 lines(tx,ty,type="l")



##################################
### Tippett's combining function:

psi = apply(partial.p,1,min)

 psi[1]				## observed value of Tippett's global test statistic

p.glob = mean(psi<=psi[1])
p.glob				## global p-value



## observed rejection region (i.e. rejection region at a significance level equal to p.glob): area below the dotted curve

lines(c(psi[1],psi[1]),c(psi[1],1),lty="dotted")
lines(c(psi[1],1),c(psi[1],psi[1]),lty="dotted")




### 1.6 ###

 require(matlab)
 source("ptest2s.r")
 load("westdata.Rdata")
 ls()


 source("ptest2s.r")
 set.seed(0)
 B <- 5000
 T <- ptest2s(Y,x,5000, "Student")
 dim(T)
 T <- T^2

 source("t2p.r")
 P <- t2p(T)
 P[1,]

## Direct combining function

 source("clostest.r")
 adjP <- clostest(T,combi="sum")
 adjP

## MaxT combining function

 adjP <- clostest(T,combi="max")
 adjP

# or

 source("stepdown.r")
 adjP <- stepdown(T)
 adjP



## Fisher's combining function

 m2lP <- -2*log(P)
 adjP <- clostest(m2lP,combi="sum")
 adjP

## Tippett's combining function

 mP <- -P		# (or: mP <- 1/P)
 adjP <- stepdown(mP)
 adjP