#  R code to simulate the exact 
#  randomization distribution of
#  F-statistics for randomized block 
#  designs. This code assumes each 
#  treatment is used exactly once in 
#  each block.  The code is posted as
#
#        randomhw.R 
#
#  Enter the responses or subject effects
#  with one row for each block.

 
  re <- matrix(c(1.4, 2.8, 1.8, 1.9,
                 2.6, 1.4, 1.6, 1.5,
                 3.5, 1.9, 2.4, 5.6,
                 3.4, 2.2, 2.8, 2.4,
                 4.9, 3.4, 3.6, 3.5,
                 5.9, 6.1, 4.4, 4.7,
                 2.4, 1.3, 1.6, 1.5,
                 4.3, 2.4, 3.7, 3.1,
                 1.5, 2.5, 1.2, 1.4,
                 4.6, 2.2, 2.9, 3.1),  
                  ncol=4, byrow=T)
  re
 

#  Compute number of blocks and number
#  of units per block

  nb <- dim(re)[1]
  ns <- dim(re)[2]
  nt <- nb*ns

# Contruct the model matrix and 
# projection operator

  b <-  kronecker(diag(rep(1,nb)),rep(1,ns))
  xb <-  cbind(matrix(1,nrow=nt,ncol=1),b) 
  xb <- xb[ , -(nb+1)]
  xt <- kronecker(rep(1,nb),diag(rep(1,ns)))
  x <- cbind(xb, xt[ , -ns])
  pb <- xb%*%solve(t(xb)%*%xb)%*%t(xb)
  px <- x%*%solve(t(x)%*%x)%*%t(x)

# Sample nsim of the possible random
# assignments of subjects to treatments
# with blocks 

  nsim <- 10000
  fdat <- matrix(0,nrow=nsim,ncol=1)  

  for(ii in 1:nsim) {
	 y <- matrix(0, nrow=nt, ncol=1)
	 nc <- 0
	 for(i in 1:nb) {
		rv <- runif(ns,0,1)
		rvr <- rank(rv)
	 for(j in 1:ns) {
		nc<-nc+1
		y[nc,1]<-re[i,rvr[j]] }}
  yhat <- px%*%y
  r <- y-yhat
  mserror <- sum(r*r)/(nt-nb-ns+1)
  fdat[ii,1] <- (t(y)%*%(px-pb)%*%y)/(ns-1)/mserror
  }  

#  Compute desired quantiles

  p <- c(.5, .9, .95, .99)
  fstats <-quantile(fdat, probs=p)
  fdist <- qf(p, ns-1, nt-nb-ns+1)

#  Print percentiles

  cbind(p, fstats, fdist)

# Compute p-value using the permuation
# distribution as a reference

  y<-matrix(t(re),ncol=1)
  yhat <- px%*%y
  r <- y-yhat
  mserror <- sum(r*r)/(nt-nb-ns+1)
  fvalue <- (t(y)%*%(px-pb)%*%y)/(ns-1)/mserror
  pvalue <- 1-length(fdat[fdat<matrix(fvalue,nrow=nsim,ncol=1)])/nsim
   
  fvalue
  pvalue

# Compute p-value using the central F
# distribution as a reference

  pvalueF <- 1-pf(fvalue, ns-1, nt-nb-ns+1)

  pvalueF


# Make a histogram 

  hist(fdat,nclass=50,probability=T, cex=1.5, 
           xlim=c(0,8), ylim=c(0,1.0),
           xlab="F-values",
           main="F-values for 5000 Permutations",
           col=0)
 
 fdat.dns1 <- density(fdat, n=200)
 fdat.dns2 <- list( x = fdat.dns1$x, 
	           y = df(fdat.dns1$x,ns-1,nt-nb-ns+1))

  lines(fdat.dns2,lty=1,lwd=4)