R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), 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 a HTML browser interface to help. Type 'q()' to quit R. > # sandm77.r November 2007, 2008 > # > # replicating Swallow & Monahan simulation study > # > # read in function code for 'estims' > source("estims.r") > # > # replication pattern > ni <- c(3,3,5,5,7,7) # this is P4 > a <- length(ni) > N <- sum(ni) > # > # set seed > set.seed(5151917 + 4 ) # for P4 > # > # 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) > 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)) > c(nrow(ests),ncol(ests)) [1] 9 700 > # save results in file for further analysis > write(ests,"estsP47.dat",ncolumns=9) # goofy tranposition! > warnings() Warning messages: 1: NA/Inf replaced by maximum positive value 2: NA/Inf replaced by maximum positive value 3: NA/Inf replaced by maximum positive value 4: NA/Inf replaced by maximum positive value 5: NA/Inf replaced by maximum positive value 6: NA/Inf replaced by maximum positive value 7: NA/Inf replaced by maximum positive value 8: NA/Inf replaced by maximum positive value 9: NA/Inf replaced by maximum positive value 10: NA/Inf replaced by maximum positive value 11: NA/Inf replaced by maximum positive value 12: NA/Inf replaced by maximum positive value 13: NA/Inf replaced by maximum positive value > # > # print a few to see what it looks like > print(t(ests)[1:5,]) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1.0497932 -0.05389674 1.0218755 3.510332e-08 0.9406769 8.197423e-10 1 [2,] 1.4032543 0.05136418 1.3715710 9.735719e-02 1.3177983 4.868822e-02 1 [3,] 1.0052913 -0.15283210 0.9987855 3.424385e-04 0.8182505 4.247799e-08 3 [4,] 1.3037456 -0.03144308 1.2462201 1.466322e-08 1.1971799 5.261036e-10 1 [5,] 0.8962539 0.05755979 0.9259866 3.653538e-02 0.8615218 3.201123e-02 1 [,8] [,9] [1,] 1 0 [2,] 1 0 [3,] 1 0 [4,] 1 0 [5,] 1 0 > print(t(ests)[695:700,]) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1.0466158 2.841449 1.0466521 2.413268 1.0047550 2.948710 1 1 5 [2,] 0.7271426 6.443425 0.7276748 6.126972 0.6983574 7.390812 1 1 5 [3,] 0.4997930 3.574199 0.4995557 2.855804 0.4794910 3.424040 1 1 5 [4,] 1.3365097 7.780763 1.3369926 6.821759 1.2833388 8.263252 1 1 5 [5,] 1.1907002 6.390817 1.1924939 5.589971 1.1441998 6.783055 1 1 5 [6,] 0.9236529 2.382180 0.9275775 2.295918 0.8890882 2.822939 1 1 5 > # some quick diagnostics > table(ests[7,],ests[8,]) # ml vs reml convergence 1 1 671 2 2 3 26 4 1 > table(ests[7,],ests[9,]) # ml convergence vs ratio 0 0.1 0.2 0.5 1 2 5 1 92 98 94 94 99 96 98 2 0 0 0 1 0 0 1 3 8 1 6 5 1 4 1 4 0 1 0 0 0 0 0 > table(ests[8,],ests[9,]) # reml convergence vs ratio 0 0.1 0.2 0.5 1 2 5 1 100 100 100 100 100 100 100 > # done > rm(list=ls()) > q()