## randomization test
## effect0 is the hypothesized treatment effect; its default is 0.
## b is the number of simulation runs; its default is 4000)

randt <- function(tr, y, effect0 = 0, b=4000)
  {
    N <- length(y)
    s1 <- (1:N)[tr == 1]
    n1 <- length(s1)
    diffobs <- mean(y[-s1])-mean(y[s1])
    y0[s1] <- y[s1]
    y0[-s1] <- y[-s1] - effect0
    
    tre <- numeric(b)
    stderror <- sqrt(var(y[s1])/length(y[s1])+var(y[-s1])/length(y[-s1]))
    print(c(diffobs,stderror))
    tobs <- diffobs/stderror


    
    for(k in 1:b)
      {
        s1 <- sample(N,n1)
        diff <- mean(y0[-s1]+effect0)-mean(y0[s1])
        stderror <-
      sqrt(var(y[s1])/length(y[s1])+var(y[-s1])/length(y[-s1]))
        tre[k] <- diff/stderror
      }

    hist(tre, freq=F)
    points(tobs,0,cex=2,pch=21,bg="blue")
    abline(v=tobs,col="blue")
    curve(dt(x,df=6),-4,4,100,add=T)
    
    pvalueup <- length(tre[tre >= tobs])/length(tre)
    pvaluedown <- length(tre[tre <= tobs])/length(tre)
    pvalue <- min(pvalueup, pvaluedown)
    print(c(effect0, pvalue))

    
  }