# sandm8.r                  November 2007, 2008
 #
 # replicating Swallow & Monahan simulation study
 #
 # read in data from sandm7.r results  *** P4 pattern, 1000 reps
 # put variable 'ests' in same structure as before
   ests <- t( read.table("estsP4.dat",header=F) )      # no header line!
 #
 # print a few to see what it looks like
   print(t(ests)[1:5,])
   print(t(ests)[6995:7000,])
 # some quick diagnostics  --- as before
   table(ests[7,],ests[8,])   # ml vs reml convergence
 # further analysis
   sigmaaf <- factor(ests[9,])  # make sigmaa a factor
   levels(sigmaaf)
 # first -- just means and variances of estimators
   tapply(ests[1,],sigmaaf,mean)
   by(t(ests[1:6,]),sigmaaf,mean)       # means by estimator, sigmaa
   by(t(ests[1:6,]),sigmaaf,var)        # covariance matrix
 #
 # compare estimates of sigma-sq-a using mse
   dseam <- (ests[2,]-ests[9,])^2 - (ests[4,]-ests[9,])^2   # anova vs ml
   meanam <- tapply(dseam,sigmaaf,mean)           # mean anova vs ml
   seam <- sqrt(tapply(dseam,sigmaaf,var)/1000)   # se a vs m
   tam <- meanam / seam                           # t statistic
   print(matrix(c(meanam,seam,tam),7,3))
   dsear <- (ests[2,]-ests[9,])^2 - (ests[6,]-ests[9,])^2   # anova vs reml
   meanar <- tapply(dsear,sigmaaf,mean)           # mean anova vs reml
   sear <- sqrt(tapply(dsear,sigmaaf,var)/1000)   # se a vs r
   tar <- meanar / sear                           # t statistic
   print(matrix(c(meanar,sear,tar),7,3))
   dsemr <- (ests[4,]-ests[9,])^2 - (ests[6,]-ests[9,])^2   # ml vs reml
   meanmr <- tapply(dsemr,sigmaaf,mean)           # mean ml vs reml
   semr <- sqrt(tapply(dsemr,sigmaaf,var)/1000)   # se m vs r
   tmr <- meanmr / semr                           # t statistic
   print(matrix(c(meanmr,semr,tmr),7,3))
 # done
  rm(list=ls())
  q()