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


>  # sandm10.r                  November 2007, 2008
>  #
>  # replicating Swallow & Monahan simulation study
>  #  -- this time, the whole thing
>  #
>  # read in function code for 'estims'
>   source("estims.r")
>  #
>  # replication patterns
>   nilist <- list( c(3,5,7), c(1,5,9), c(3,3,5,5,7,7), c(1,1,5,5,9,9),
+          c(1,1,1,1,13,13), c(3,3,3,5,5,5,7,7,7), c(1,1,1,5,5,5,9,9,9),
+          c(1,1,1,1,1,1,1,19,19), c(2,10,18), c(3,15,27) )
>  #
>  # initialize estall
>    estall <- matrix(0,10,0)
>  # set seed
>    set.seed(5151917)          # usual seed
>  #
>  #                            main loop start
>   for (pattern in (1:10)) {
+  # 
+   ni <- nilist[[pattern]]
+   a <- length(ni)
+   N <- sum(ni)
+   print(ni)
+  #
+  # generate data
+    nrep <- 100      # only a few replicates
+    sigmaa <- rep(c(0,1,2,5,10,20,50)/10,rep(nrep,7)) # set true values
+    sigmae <- 1      # of course
+  #
+  # first mean vectors
+    ybars <- matrix(rnorm(a*nrep*7),a,nrep*7)      # 7 levels of sigmaa
+    ybars <- ybars*sqrt(rep(sigmaa,a)+sigmae/ni)
+  #
+  # now sse/1 has chi-square distribution
+    sses <- rchisq(nrep*7,(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)
+    ests <- rbind(ests,unname(sigmaa),rep(pattern,ncol(ests)))
+    print(pattern,nrow(ests),ncol(ests))
+    estall <- cbind(estall,ests)
+                           }   # end of pattern loop       
[1] 3 5 7
[1] 1
[1] 1 5 9
[1] 2
[1] 3 3 5 5 7 7
[1] 3
[1] 1 1 5 5 9 9
[1] 4
[1]  1  1  1  1 13 13
[1] 5
[1] 3 3 3 5 5 5 7 7 7
[1] 6
[1] 1 1 1 5 5 5 9 9 9
[1] 7
[1]  1  1  1  1  1  1  1 19 19
[1] 8
[1]  2 10 18
[1] 9
[1]  3 15 27
[1] 10
>  #
>  # save results in file for further analysis
>    write(estall,"estsall.dat",ncolumns=10)       # goofy tranposition!
>  #
>  # print a few to see what it looks like
>    print(t(estall)[1:5,])
          [,1]        [,2]      [,3]         [,4]      [,5]         [,6] [,7]
[1,] 2.0957595 -0.43731542 1.7005747 1.900822e-07 1.6800472 1.990139e-07    1
[2,] 0.9800917 -0.08707874 1.0099923 1.219897e-01 0.8597961 1.332397e-09    1
[3,] 1.2399172 -0.01093755 1.2127099 6.653326e-02 1.1155268 5.327370e-02    1
[4,] 0.9948197  0.14798138 1.0422754 1.591407e-01 0.9424845 1.415562e-01    1
[5,] 0.6883350  0.18922916 0.7040313 1.220734e-01 0.6272762 1.914745e-01    1
     [,8] [,9] [,10]
[1,]    1    0     1
[2,]    1    0     1
[3,]    1    0     1
[4,]    1    0     1
[5,]    1    0     1
>  # some quick diagnostics
>    table(estall[7,],estall[8,])   # ml vs reml convergence
   
       1    2    4
  1 6632    3    2
  2   23    0    0
  3  327    0    0
  4   13    0    0
>    table(estall[7,],estall[9,])   # ml convergence vs ratio
   
      0 0.1 0.2 0.5   1   2   5
  1 941 964 934 944 934 946 974
  2   5   3   5   3   4   0   3
  3  51  33  58  52  61  50  22
  4   3   0   3   1   1   4   1
>    table(estall[8,],estall[9,])   # reml convergence vs ratio
   
       0  0.1  0.2  0.5    1    2    5
  1 1000  998  998  999 1000 1000 1000
  2    0    1    1    1    0    0    0
  4    0    1    1    0    0    0    0
>    table(estall[7,],estall[10,])  # ml convergence vs pattern
   
      1   2   3   4   5   6   7   8   9  10
  1 658 659 667 649 679 664 657 683 656 665
  2   4   2   2   2   1   3   3   1   2   3
  3  38  36  30  48  20  31  38  16  38  32
  4   0   3   1   1   0   2   2   0   4   0
>    table(estall[8,],estall[10,])  # reml convergence vs pattern
   
      1   2   3   4   5   6   7   8   9  10
  1 698 699 700 699 699 700 700 700 700 700
  2   2   1   0   0   0   0   0   0   0   0
  4   0   0   0   1   1   0   0   0   0   0
>  # done
>   rm(list=ls())
>   q()