# This R code simulates the exact # randomization distribution of an # F-statistic for a randomized block # design. This code assumes each # treatment is used exactly once in # each block. The code is posted as # # randomblock.ssc # # Enter the responses or subject effects # with one row for each block. re <- matrix(c(1.4, 2.8, 1.8, 2.6, 1.4, 1.6, 3.5, 1.9, 2.4, 3.4, 2.2, 2.8, 4.9, 3.4, 3.6, 5.9, 6.1, 4.4, 2.4, 1.3, 1.6, 4.3, 2.4, 3.7, 1.5, 2.5, 1.2, 4.6, 2.2, 2.9), ncol=3, 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 <- 5000 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