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. > # 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,]) V1 V2 V3 V4 V5 V6 V7 V8 [1,] 1.1764660 -0.11677250 1.1077960 1.890532e-07 1.0420150 1.562887e-07 1 1 [2,] 0.6690832 0.15195630 0.6923927 2.662505e-01 0.6341990 1.588799e-01 1 1 [3,] 0.8028353 -0.14687910 0.8822460 3.491727e-02 0.6562857 3.651575e-08 3 1 [4,] 0.8546202 0.01440330 0.9190613 1.421776e-02 0.8175840 2.484703e-02 3 1 [5,] 1.1639300 -0.02872148 1.1432370 6.912314e-08 1.1017090 1.100769e-08 1 1 V9 [1,] 0 [2,] 0 [3,] 0 [4,] 0 [5,] 0 > print(t(ests)[6995:7000,]) V1 V2 V3 V4 V5 V6 V7 V8 V9 [1,] 1.0458840 2.044095 1.0508160 1.957027 1.0071050 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.4594520 8.167333 1.4585300 6.401599 1.4005420 7.747390 1 1 5 [5,] 0.6874366 9.528379 0.6874527 8.350344 0.6599455 10.056690 1 1 5 [6,] 0.9402320 4.442200 0.9404206 3.716604 0.9027318 4.511552 1 1 5 > # some quick diagnostics --- as before > table(ests[7,],ests[8,]) # ml vs reml convergence 1 2 1 6652 2 2 21 0 3 314 0 4 11 0 > # further analysis > sigmaaf <- factor(ests[9,]) # make sigmaa a factor > levels(sigmaaf) [1] "0" "0.1" "0.2" "0.5" "1" "2" "5" > # first -- just means and variances of estimators > tapply(ests[1,],sigmaaf,mean) 0 0.1 0.2 0.5 1 2 5 0.9980696 1.0063280 0.9990151 0.9743010 1.0063603 0.9804888 0.9818066 > by(t(ests[1:6,]),sigmaaf,mean) # means by estimator, sigmaa INDICES: 0 V1 V2 V3 V4 V5 V6 0.998069642 -0.003925977 0.971602399 0.055983966 0.918786989 0.054178490 ------------------------------------------------------------ INDICES: 0.1 V1 V2 V3 V4 V5 V6 1.00632804 0.09290774 0.99330665 0.11938695 0.94412497 0.12994841 ------------------------------------------------------------ INDICES: 0.2 V1 V2 V3 V4 V5 V6 0.9990151 0.1930901 0.9948929 0.1955347 0.9468991 0.2175725 ------------------------------------------------------------ INDICES: 0.5 V1 V2 V3 V4 V5 V6 0.9743010 0.4916847 0.9745506 0.4323460 0.9319278 0.5060546 ------------------------------------------------------------ INDICES: 1 V1 V2 V3 V4 V5 V6 1.0063603 0.9870452 1.0071189 0.8325217 0.9646306 0.9943912 ------------------------------------------------------------ INDICES: 2 V1 V2 V3 V4 V5 V6 0.9804888 2.0155536 0.9806297 1.6899055 0.9407852 2.0341334 ------------------------------------------------------------ INDICES: 5 V1 V2 V3 V4 V5 V6 0.9818066 5.0553580 0.9819274 4.2013060 0.9426338 5.0665903 > by(t(ests[1:6,]),sigmaaf,var) # covariance matrix INDICES: 0 V1 V2 V3 V4 V5 V1 0.084821949 -0.016838545 0.071401295 -0.009101373 0.073571509 V2 -0.016838545 0.018585512 -0.011153724 0.008925296 -0.010129649 V3 0.071401295 -0.011153724 0.063662960 -0.006943388 0.063777462 V4 -0.009101373 0.008925296 -0.006943388 0.006663085 -0.006893031 V5 0.073571509 -0.010129649 0.063777462 -0.006893031 0.066269353 V6 -0.007043039 0.011674998 -0.005124063 0.007059218 -0.004572478 V6 V1 -0.007043039 V2 0.011674998 V3 -0.005124063 V4 0.007059218 V5 -0.004572478 V6 0.009581218 ------------------------------------------------------------ INDICES: 0.1 V1 V2 V3 V4 V5 V1 0.082687323 -0.012876520 0.074048584 -0.009389201 0.074587397 V2 -0.012876520 0.038435132 -0.008533115 0.024310700 -0.007003113 V3 0.074048584 -0.008533115 0.068878952 -0.007276773 0.068266898 V4 -0.009389201 0.024310700 -0.007276773 0.018666066 -0.006751924 V5 0.074587397 -0.007003113 0.068266898 -0.006751924 0.069183794 V6 -0.007103134 0.032413094 -0.005007336 0.022249239 -0.003984103 V6 V1 -0.007103134 V2 0.032413094 V3 -0.005007336 V4 0.022249239 V5 -0.003984103 V6 0.030042947 ------------------------------------------------------------ INDICES: 0.2 V1 V2 V3 V4 V5 V1 0.084661800 -0.013045761 0.079044360 -0.012088671 0.078513253 V2 -0.013045761 0.061431095 -0.010518373 0.041862346 -0.008313952 V3 0.079044360 -0.010518373 0.075519728 -0.010685311 0.074273808 V4 -0.012088671 0.041862346 -0.010685311 0.032870769 -0.009609443 V5 0.078513253 -0.008313952 0.074273808 -0.009609443 0.074179848 V6 -0.009680606 0.056399825 -0.008309212 0.040221541 -0.006674916 V6 V1 -0.009680606 V2 0.056399825 V3 -0.008309212 V4 0.040221541 V5 -0.006674916 V6 0.054527604 ------------------------------------------------------------ INDICES: 0.5 V1 V2 V3 V4 V5 V6 V1 0.08282837 -0.01691962 0.07992852 -0.01721540 0.07841137 -0.01526536 V2 -0.01691962 0.20596649 -0.01640713 0.15600709 -0.01427926 0.20095915 V3 0.07992852 -0.01640713 0.07777868 -0.01667577 0.07604404 -0.01476439 V4 -0.01721540 0.15600709 -0.01667577 0.12568547 -0.01523737 0.15711023 V5 0.07841137 -0.01427926 0.07604404 -0.01523737 0.07470939 -0.01294329 V6 -0.01526536 0.20095915 -0.01476439 0.15711023 -0.01294329 0.20317047 ------------------------------------------------------------ INDICES: 1 V1 V2 V3 V4 V5 V6 V1 0.08406980 -0.02539715 0.08230381 -0.02421038 0.08037401 -0.02419205 V2 -0.02539715 0.64673026 -0.02640916 0.48726475 -0.02314426 0.61210175 V3 0.08230381 -0.02640916 0.08140619 -0.02463178 0.07894600 -0.02474929 V4 -0.02421038 0.48726475 -0.02463178 0.38557784 -0.02221970 0.48004694 V5 0.08037401 -0.02314426 0.07894600 -0.02221970 0.07715854 -0.02185688 V6 -0.02419205 0.61210175 -0.02474929 0.48004694 -0.02185688 0.60494344 ------------------------------------------------------------ INDICES: 2 V1 V2 V3 V4 V5 V6 V1 0.08593064 -0.04282625 0.08541890 -0.04219130 0.08252589 -0.04570876 V2 -0.04282625 2.02582423 -0.04427003 1.59244705 -0.04105027 1.95754380 V3 0.08541890 -0.04427003 0.08501624 -0.04243319 0.08208882 -0.04597450 V4 -0.04219130 1.59244705 -0.04243319 1.33929955 -0.04000980 1.64373191 V5 0.08252589 -0.04105027 0.08208882 -0.04000980 0.07932332 -0.04309838 V6 -0.04570876 1.95754380 -0.04597450 1.64373191 -0.04309838 2.02410190 ------------------------------------------------------------ INDICES: 5 V1 V2 V3 V4 V5 V6 V1 0.081315716 0.006266427 0.081151829 0.001746343 0.078065081 0.007123342 V2 0.006266427 11.477818942 0.004179657 8.885849744 0.004507188 10.734928425 V3 0.081151829 0.004179657 0.081020838 0.001358534 0.077925696 0.006716007 V4 0.001746343 8.885849744 0.001358534 7.469546235 0.001281145 9.030953724 V5 0.078065081 0.004507188 0.077925696 0.001281145 0.074959157 0.006366405 V6 0.007123342 10.734928425 0.006716007 9.030953724 0.006366405 10.924497606 > # > # 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)) [,1] [,2] [,3] [1,] 0.008791714 0.0006103175 14.405147 [2,] 0.019423743 0.0011976057 16.218813 [3,] 0.028559574 0.0018734293 15.244544 [4,] 0.075692825 0.0064399611 11.753615 [5,] 0.233010100 0.0295878453 7.875197 [6,] 0.589921450 0.0776233421 7.599795 [7,] 3.369416762 0.4753376907 7.088470 > 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)) [,1] [,2] [,3] [1,] 0.0060753946 0.0005634284 10.78290386 [2,] 0.0075371854 0.0010462120 7.20426219 [3,] 0.0066355408 0.0015490292 4.28367704 [4,] 0.0028257134 0.0054465185 0.51881095 [5,] 0.0418814033 0.0195359135 2.14381597 [6,] 0.0007974293 0.0562719886 0.01417098 [7,] 0.5513982563 0.3536124377 1.55932936 > 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)) [,1] [,2] [,3] [1,] -0.002716319 0.0006004971 -4.523451 [2,] -0.011886557 0.0011558839 -10.283522 [3,] -0.021924033 0.0017818667 -12.303969 [4,] -0.072867111 0.0065881906 -11.060262 [5,] -0.191128696 0.0175208004 -10.908674 [6,] -0.589124021 0.0625784560 -9.414167 [7,] -2.818018506 0.3251281022 -8.667410 > # done > rm(list=ls()) > q()