#####################################################
#####################################################
InteractionFollowup <- function(data, factorA, factorB, whichLevel, MSerror, DFerror, alpha)
#e.g., InteractionFollowup(DF$bin_tot, DF$go_gng, DF$speed_condition, "go", 81382.29, 24, .05) 
#use detailed=TRUE arg for ezANOVA to get MSerror
{
	rows <- factorA==whichLevel
	data <- data[rows]
	factorB <- factorB[rows]

	#count number of observations per cell
	numObs <- tapply(data[rows],list(factorB), length)
	numGroups <- length(unique(factorB))
	qval = qtukey(1-alpha, numGroups, DFerror)

  	#Get the mean sums of squares for factor B at this level
  	means = tapply(data,list(factorB),mean)
  	SS <- t(scale(means, scale=F))%*%scale(means, scale=F)*min(numObs)
  	MS <- SS/(numGroups-1)
  	#Compute f-test for factor B
  	Fval <- MS/MSerror
  	p <- 1-pf(Fval,numGroups-1, DFerror)
  	print(paste('Fval', Fval))
  	#print(Fval)
  	print(paste('MS', MS))
  	#print(MS)
  	print(paste('SS', SS))
  	#print(SS)
  	print(paste('p', p))
  	#print(p)
  	#Compute p-val
  	
  	#Compute Tukey - Kramer test
  	critValMatrix = matrix(0,nrow=numGroups, ncol=numGroups)
  
  	for (rowCount in 1:numGroups)
  	{
  		for (colCount in 1:numGroups)
  		{
  			critVal = qval * sqrt(MSerror/2 * (1/numObs[rowCount] + 1/numObs[colCount]))
  			critValMatrix[rowCount, colCount] = critVal
  		}
  	}		
  	
  	diffmatrix = t(means%*%t(rep(1,numGroups)))-means%*%t(rep(1,numGroups))
  	sigDiffMatrix = abs(diffmatrix) > critVal
  	print('')
	print('diffmatrix')
      print(diffmatrix)
      print(paste('critVal', critVal))
      #print(critVal)
      print('sigDiffMatrix')
  	print(sigDiffMatrix)
}
#####################################################