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. > # sandm7.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 <- 1000 # 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) # matrix of std normals > ybars <- ybars*sqrt(rep(sigmaa,a)+sigmae/ni) #***weird***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*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)) # add row of sigmaa level > c(nrow(ests),ncol(ests)) [1] 9 7000 > # save results in file for further analysis > write(ests,"estsP4.dat",ncolumns=9) # goofy tranposition! > warnings() Warning messages: 1: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 2: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 3: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 4: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 5: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 6: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 7: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 8: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 9: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 10: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 11: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 12: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 47: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 48: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 49: In nlm(u2ll, log(modanova)) : NA/Inf replaced by maximum positive value 50: In nlm(u2ll, log(modanova)) : 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.1764664 -0.11677249 1.1077960 1.890532e-07 1.0420153 1.562887e-07 1 [2,] 0.6690832 0.15195628 0.6923927 2.662505e-01 0.6341990 1.588799e-01 1 [3,] 0.8028353 -0.14687911 0.8822460 3.491727e-02 0.6562857 3.651575e-08 3 [4,] 0.8546202 0.01440330 0.9190613 1.421776e-02 0.8175841 2.484703e-02 3 [5,] 1.1639303 -0.02872148 1.1432367 6.912314e-08 1.1017089 1.100769e-08 1 [,8] [,9] [1,] 1 0 [2,] 1 0 [3,] 1 0 [4,] 1 0 [5,] 1 0 > print(t(ests)[6995:7000,]) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1.0458838 2.044095 1.0508156 1.957027 1.0071047 2.406670 1 1 5 [2,] 0.8865204 3.816670 0.8869188 3.128988 0.8513084 3.804852 1 1 5 [3,] 0.4497424 5.377049 0.4497681 4.451575 0.4317212 5.352820 1 1 5 [4,] 1.4594521 8.167333 1.4585298 6.401599 1.4005417 7.747390 1 1 5 [5,] 0.6874366 9.528379 0.6874527 8.350344 0.6599455 10.056693 1 1 5 [6,] 0.9402320 4.442200 0.9404206 3.716604 0.9027318 4.511552 1 1 5 > # some quick diagnostics > table(ests[7,],ests[8,]) # ml vs reml convergence 1 2 1 6645 2 2 19 0 3 309 0 4 25 0 > table(ests[7,],ests[9,]) # ml convergence vs ratio 0 0.1 0.2 0.5 1 2 5 1 943 936 949 934 959 959 967 2 4 3 2 2 4 2 2 3 50 57 48 60 33 35 26 4 3 4 1 4 4 4 5 > table(ests[8,],ests[9,]) # reml convergence vs ratio 0 0.1 0.2 0.5 1 2 5 1 999 1000 1000 1000 1000 999 1000 2 1 0 0 0 0 1 0 > # done > rm(list=ls()) > q()