code

########## for a=3 and b=8############

#clear all; format compact; format long;

iter=5000
a = 3
b = 8
N = a*b
dfE = (a-1)*(b-1)

blockvec <-c(rep(1,a),rep(2,a),rep(3,a),rep(4,a),
             rep(5,a),rep(6,a),rep(7,a),rep(8,a))
trtvec <-c(rep(1:a,b)) 
mu = 0
tau = matrix(0,nrow=1,ncol=a)
beta= matrix(0,nrow=1,ncol=b)


EY<-matrix(0,nrow=8,ncol=3)

for (j in 1:b){
  for (i in 1:a){
    EY[j,i] = mu + tau[i] + beta[j]
  }
}

#EY

Reject     = matrix(0,nrow=1,ncol=3)
RejectMLE  = matrix(0,nrow=1,ncol=3)
RejectExact= matrix(0,nrow=1,ncol=3)

for (samp in 1:iter){
  block = blockvec
  trt   = trtvec
  #Generate data set
  #Ymat = EY + matrix( rnorm(b*a,mean=0,sd=1), b, a)
  Ymat = EY + randn(b,a)
  Y = matrix(t(Ymat),nrow = 1,ncol=N, byrow = TRUE)
  
  # Calculate summary stats needed for MLE computation
  yidot<-matrix(0,nrow=a,ncol=1)
  for (i in 1:a){
    yidot[i,1] = sum(Ymat[,i])
  }  
  yjdot<-matrix(0,nrow=b,ncol=1)
  for (j in 1:b){
    yjdot[j,1] = sum(Ymat[j,]) 
  } 
  ydotdot = sum(yjdot)
  
  # Run ANOVA for complete data set and save p-value
  #p = anovan(Y,{block,trt},'varnames',{'block','trt'},'display','off')
Y=t(Y)
  p = aov(Y~as.factor(trt)+as.factor(block))
  p=c(summary(p)[[1]][[1,"Pr(>F)"]],summary(p)[[1]][[2,"Pr(>F)"]])
  pvalue=p[1]
  
  # Delete an observation 
  row = ceil(b*rand(1))
  col = ceil(a*rand(1))
  remove = Ymat[row,col]
  Ymat[row,col] = NaN
  Y = matrix(t(Ymat),nrow = 1,ncol=N, byrow = TRUE)
  
  vecID = a*(row-1) + col # location of NaN in block and trt vectors
  block[vecID] = NaN
  trt[vecID] = NaN
  
  # Run exact ANOVA and save p-value
  Y=t(Y)
  p = aov(Y~as.factor(block)+as.factor(trt))
  p=c(summary(p)[[1]][[1,"Pr(>F)"]],summary(p)[[1]][[2,"Pr(>F)"]])
  pvalueExact=p[2]

  
  # Adjust summary stats for missing value and calculate MLE
  ydotdotP = ydotdot    - remove
  yidotP   = yidot[col] - remove
  yjdotP   = yjdot[row] - remove
  
  MLE = (a*yidotP + b*yjdotP - ydotdotP)/dfE
  Y[vecID] = MLE
  block[vecID] = row
  trt[vecID] = col
  
  # Run ANOVA with MLE
  p = aov(Y~as.factor(trt)+as.factor(block))
  #SS=c(summary(p)[[1]][[1,"Sum Sq"]],summary(p)[[1]][[2,"Sum Sq"]],summary(p)[[1]][[3,"Sum Sq"]])
  MStrt = summary(p)[[1]][[1,"Mean Sq"]]
  SSE   = summary(p)[[1]][[3,"Sum Sq"]]
  MSE   = SSE/(dfE-1)
  Ftrt = MStrt/MSE
  pvalueMLE = pf(Ftrt,(a-1),(dfE-1),lower.tail=FALSE)
  
  if( pvalue <= .10) Reject[1,1] = Reject[1,1]+1
  if (pvalue <= .05) Reject[1,2] = Reject[1,2]+1
  if (pvalue <= .01) Reject[1,3] = Reject[1,3]+1
  
  if (pvalueMLE <= .10) RejectMLE[1,1] = RejectMLE[1,1]+1
  if (pvalueMLE <= .05) RejectMLE[1,2] = RejectMLE[1,2]+1
  if (pvalueMLE <= .01) RejectMLE[1,3] = RejectMLE[1,3]+1
  
  if (pvalueExact <= .10) RejectExact[1,1] = RejectExact[1,1]+1
  if (pvalueExact <= .05) RejectExact[1,2] = RejectExact[1,2]+1
  if (pvalueExact <= .01) RejectExact[1,3] = RejectExact[1,3]+1
  
}

Power8 = Reject/iter
PowerMLE8 = RejectMLE/iter
PowerExact8 = RejectExact/iter