# sandm6.r                  November 2007, 2008
 #
 # replicating Swallow & Monahan simulation study
 #
  ni <- c(3,3,5,5,7,7)      # this is P4
  a <- length(ni)
  N <- sum(ni)
 #
 # ****************************************************************
 # function for computing estimators -- expects a+1 values as input
  estims <- function(ybsse) { # compute all estimators
     sse <- ybsse[a+1]        # last component is SSE
     yb <- ybsse[1:a]         # sample means
     ybg <- sum(yb*ni)/N      # grand mean or ybar..
 # first compute ANOVA estimators
     ssa <- sum(ni*(yb-ybg)*(yb-ybg))
     sgh1 <- sse/(N-a)
     sgh2 <- (ssa-(a-1)*sse/(N-a))/(N-sum(ni*ni)/N)
     eanova <-c(sgh1,sgh2)    # used as starters for others       
     modanova <- c(sgh1,abs(sgh2))
 #
 # function for ML estimators -- computed iteratively
   # define ML log-likelihood function (-2*log-likelihood)
   # reparameterize so always positive
      u2ll <- function(vc) { 
        vc1 = exp(vc[1])
        vc2 = exp(vc[2])
        muhat <- sum(ni*(vc1/(vc1+ni*vc2))*yb)/
                 sum(ni*(vc1/(vc1+ni*vc2)))
        u2ll <- (N-a)*log(vc1) + sum(log(vc1+ni*vc2)) +
                  sse/vc1 + 
      sum(ni*(yb-muhat)*(yb-muhat)/(vc1+ni*vc2)) }
   # maximize ML log-likelihood
      ml <- nlm(u2ll, log(modanova)) # all defaults -- no printing
    eml <- exp(as.vector(ml$estimate))
 #
 # now reml is not much different  
   # define ML log-likelihood function (-2*log-likelihood)
      u2reml <- function(vc) { 
        vc1 = exp(vc[1])
        vc2 = exp(vc[2])
        muhat <- sum(ni*(vc1/(vc1+ni*vc2))*yb)/
                 sum(ni*(vc1/(vc1+ni*vc2)))
        u2reml <- (N-a)*log(vc1) + sum(log(vc1+ni*vc2)) +
          log( sum(ni*(vc1/(vc1+ni*vc2))) ) + sse/vc1 + 
      sum(ni*(yb-muhat)*(yb-muhat)/(vc1+ni*vc2)) }
   # maximize REML 
     rml <- nlm(u2reml, log(modanova) ) # all defaults -- no printing
    reml <- exp(as.vector(rml$estimate))
 #
 # done with estimation
 # put three pairs of estimators together with warnings
    estims <- c(eanova,eml,reml,ml$code,rml$code)  }  # end of estimation
 # *********************************************************************
 #
 #
 # generate data
   nrep <- 10000
   sigmaa <- 2      # set true value
   sigmae <- 1      # of course
 #
 # set seed
   set.seed(5151917 + 4)
 # first mean vectors
   ybars <- matrix(rnorm(a*nrep),a,nrep)  # matrix of standard normals
   ybars <- ybars*sqrt(sigmaa+sigmae/ni)  #***weird***
 # each col of ybars is MVN with cov mx (sigmaa*I+diag(sigmae/ni) )
 #
 # now sse/1 has chi-square distribution
   sses <- rchisq(nrep,(N-a)) # with (N-a) df
 #
 # put ybar and sse together as single input vector
 #   but columns of matrix  
   data <- rbind(ybars,unname(sses))    # drop labels
 # 
 # apply function to cols of matrix
   ests <- apply(data,2,estims)
   c(nrow(ests),ncol(ests))
 # save all results in file for further analysis
   write(ests,"ests20P4.dat",ncolumns=8)       # goofy!
   warnings()
 #w print(t(ests))
 #
 # some simple analysis
   apply(ests,1,mean)
   cov(t(ests))
   table(ests[7,],ests[8,])
 # done
  rm(list=ls())
  q()