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. > # 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.901050e-07 1.6800472 1.990391e-07 1 [2,] 0.9800917 -0.08707874 1.0099923 1.219897e-01 0.8597961 1.338487e-09 1 [3,] 1.2399172 -0.01093755 1.2127099 6.653325e-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 6633 5 2 2 29 0 0 3 313 0 0 4 18 0 0 > table(estall[7,],estall[9,]) # ml convergence vs ratio 0 0.1 0.2 0.5 1 2 5 1 942 966 934 945 933 947 973 2 5 4 8 3 4 1 4 3 50 29 55 50 61 45 23 4 3 1 3 2 2 7 0 > table(estall[8,],estall[9,]) # reml convergence vs ratio 0 0.1 0.2 0.5 1 2 5 1 998 999 998 998 1000 1000 1000 2 2 0 1 2 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 660 661 665 648 678 664 657 684 656 667 2 3 3 4 3 2 3 5 1 2 3 3 37 32 30 46 18 31 35 15 39 30 4 0 4 1 3 2 2 3 0 3 0 > table(estall[8,],estall[10,]) # reml convergence vs pattern 1 2 3 4 5 6 7 8 9 10 1 698 697 700 699 699 700 700 700 700 700 2 2 3 0 0 0 0 0 0 0 0 4 0 0 0 1 1 0 0 0 0 0 > # done > rm(list=ls()) > q()