# estims.r                  November 2007, 2008
 #
 # replicating Swallow & Monahan simulation study
 #
 # ****************************************************************
 # 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/(vc[1]+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
 # *********************************************************************
 #