### 3.5 ###

# load data matrices, rows=endpts, cols=animals 

 S1 <- read.csv(file="http://lanfiles.williams.edu/~bklingen/SMH/S1.csv",header=FALSE) 

#matrix of ordinal observations under control

 S2 <- read.csv(file="http://lanfiles.williams.edu/~bklingen/SMH/S2.csv",header=FALSE) 

#matrix of ordinal observations under treatment

 source("http://lanfiles.williams.edu/~bklingen/SMH/SMH.R") 

#source necessary functions




 fob.perm<-perm.SMH(S1,S2,c=4,enum=TRUE,dom.index=c(5,9,14,17,22,25)) 

  fob.perm$sig.meanT 
# global test using W_0' as test stat

 fob.perm$adj.P 
# adjusted P-values (via maxT step down method) for the 25 


# endpoints

 fob.perm$adj.P.domain 
# adjusted P-values (via full closed testing) for the domains 

 hist(fob.perm$max.T) 
# histogram of max_h z_h


# If complete enumeration is too computationally demanding, use

 fob.perm<-perm.SMH(S1,S2,c=4,perms=10000,dom.index=c(5,9,14,17,22,25))


# using maximum scores:

 fob.perm.maxscores<-perm.SMH(S1,S2,c=4,enum=TRUE,scores="max",dom.index= c(5,9,14,17,22,25))


 fob.perm.userscores <- perm.SMH(S1,S2,c=4,enum=TRUE,scores=c(1,3,4,5),dom.index=c(5,9,14,17,22,25))

#bootstrap analysis with centering z_h:

 fob.boot <- boot.SMH(S1,S2,c=4,boots=10000,dom.index=c(5,9,14,17,22,25))
 fob.boot$sig.meanT 

#global test using W0` and centering z_h
 fob.boot$adj.P 

#adjusted P-values (via maxT step down method) for the 25 endpoints

fob.boot$adj.P.domain 
#adjusted P-values (via full closed testing) for the 6 domains



#bootstrap analysis with centering s_h and then standardizing:

fob.boot <- boot.SMH(S1,S2,c=4,boots=10000,dom.index=c(5,9,14,17,22,25), centering=2)

fob.boot$sig.meanT 
#global test using W0` and centering s_h

fob.boot$adj.P 
#adjusted P-values (via maxT step down method) for the 25 endpoints

fob.boot$adj.P.domain 
#adjusted P-values (via full closed testing) for the 6 domains