R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>  # 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))
[1]     8 10000
>  # save all results in file for further analysis
>    write(ests,"ests20P4.dat",ncolumns=8)       # goofy!
>    warnings()
NULL
>  #w print(t(ests))
>  #
>  # some simple analysis
>    apply(ests,1,mean)
[1] 1.0020020 1.9727709 1.0017125 1.6121521 0.9615894 1.9890997 1.0000000
[8] 1.0000000
>    cov(t(ests))
            [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
[1,]  0.08207907 -0.01793168  0.08184332 -0.01905400  0.07861220 -0.01862517
[2,] -0.01793168  1.97403583 -0.01806848  1.56243486 -0.01719051  1.86982681
[3,]  0.08184332 -0.01806848  0.08189211 -0.01840254  0.07856116 -0.01758549
[4,] -0.01905400  1.56243486 -0.01840254  1.31777194 -0.01780332  1.58081095
[5,]  0.07861220 -0.01719051  0.07856116 -0.01780332  0.07540811 -0.01715138
[6,] -0.01862517  1.86982681 -0.01758549  1.58081095 -0.01715138  1.89704907
[7,]  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000
[8,]  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000
     [,7] [,8]
[1,]    0    0
[2,]    0    0
[3,]    0    0
[4,]    0    0
[5,]    0    0
[6,]    0    0
[7,]    0    0
[8,]    0    0
>    table(ests[7,],ests[8,])
   
        1
  1 10000
>  # done
>   rm(list=ls())
>   q()