# TODO: Miscellaneous test functions for the VCA-package # # Author: schueta6 ############################################################################### ### load all testdata data(dataEP05A2_1) data(dataEP05A2_2) data(dataEP05A2_3) cat("\n\n***********************************************************************") cat("\nVariance Component Analysis (VCA) - test cases defined in runit.misc.R.") cat("\n***********************************************************************\n\n") # test numerical equivalence of the covariance matrix of estimated variance components against SAS PROC MIXED # for a balanced design when following SAS-code is used: # # proc mixed data=sample1 method=type1 asycov cl; # class lot device day run; # model y=; # random lot device lot*device*day lot*device*day*run; # run; TF001.CovarianceMatrix <- function() { data(VCAdata1) sample1 <- VCAdata1[which(VCAdata1$sample==1),] sample1$device <- gl(3,28,252) # add device variable set.seed(505) sample1$y <- sample1$y + rep(rep(rnorm(3,,.25), c(28,28,28)),3) # add error component according to res1 <- anovaVCA(y~(lot+device)/day/run, sample1) # request A-matrices required for convariance matrix of VCs inf1 <- VCAinference(res1, VarVC=TRUE, constrainCI=FALSE) checkEquals(round(inf1$VCAobj$aov.tab[2, "Var(VC)"], 6), 0.000301) # Var(VC_lot) checkEquals(round(inf1$VCAobj$aov.tab[3, "Var(VC)"], 6), 0.004091) # Var(VC_device) checkEquals(round(inf1$VCAobj$aov.tab[4, "Var(VC)"], 6), 0.000071) # Var(VC_day) checkEquals(round(inf1$VCAobj$aov.tab[5, "Var(VC)"], 6), 0.000084) # Var(VC_run) checkEquals(round(inf1$VCAobj$aov.tab[6, "Var(VC)"], 11), 2.108e-8) # Var(VC_error) } # For fitting the Orthodont data from R-package 'nlme' following SAS PROC MIXED code is used: # # proc mixed data=ortho method=type1; # class subject sex; # model distance = sex sex*age2/solution outp=outp outpm=outpm residual; # random subject subject*age2; # run; # # age2 is equal to age - 11 (the centered version) # test marginal residuals against SAS PROC MIXED implementation # # Here the Orthodont-dataset from package 'nlme' is used (if available) # # edit (schueta 2017-01-23): newly introduced arguement 'order.data' needed to be set to FALSE TF002.marginal_residuals <- function() { data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 fit.R <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, order.data=FALSE) sas.raw <- c(3.384375,0.815625,3.246875,3.678125,-1.115625,-1.684375,-2.753125,-0.821875,0.384375,-1.684375,-1.753125,0.178125,2.884375,3.315625,0.746875,-0.321875,-2.615625,-0.684375,-3.253125,-1.321875,1.884375,1.315625,1.246875,1.178125,-0.615625,-2.184375,-1.253125,-0.821875,1.384375,-2.684375,-1.253125,-1.821875,0.384375,-3.684375,5.246875,-1.321875,4.884375,3.815625,5.246875,4.178125,0.384375,-1.184375,-2.253125,-2.321875,-1.115625,-0.684375,-1.753125,0.678125,-5.615625,0.315625,0.246875,2.178125,-0.115625,1.315625,-0.253125,-1.321875,0.384375,0.315625,0.246875,2.678125,-0.615625,-2.684375,-2.253125,-2.321875,-0.209090909,-2.168181818,-1.627272727,-1.086363636,-0.209090909,-0.668181818,0.8727272727,1.4136363636,-0.709090909,1.8318181818,1.3727272727,1.9136363636,2.2909090909,2.3318181818,1.8727272727,2.4136363636,0.2909090909,0.8318181818,-0.627272727,-0.586363636,-1.209090909,-1.168181818,-2.127272727,-1.586363636,0.2909090909,0.3318181818,-0.127272727,0.9136363636,1.7909090909,0.8318181818,0.3727272727,-0.086363636,-1.209090909,-1.168181818,-1.127272727,-2.586363636,-4.709090909,-3.168181818,-4.127272727,-4.586363636,3.2909090909,2.8318181818,4.8727272727,3.9136363636) sas.stu <- c(1.5050944191,0.37015662,1.473535357,1.6357304998,-0.496139158,-0.764423058,-1.249455871,-0.365503077,0.1709387013,-0.764423058,-0.795623999,0.0792154957,1.2827351328,1.5047362981,0.3389556788,-0.143143791,-1.163217016,-0.310591187,-1.476371806,-0.587862363,0.8380165602,0.5970725556,0.5658716145,0.5239340683,-0.273779871,-0.991338994,-0.568708064,-0.365503077,0.6156572739,-1.218254929,-0.568708064,-0.810221649,0.1709387013,-1.672086801,2.3811990995,-0.587862363,2.1721722779,1.7316522338,2.3811990995,1.858089786,0.1709387013,-0.537507123,-1.022539935,-1.032580936,-0.496139158,-0.310591187,-0.795623999,0.301574782,-2.497372734,0.1432406844,0.1120397432,0.9686526409,-0.051420585,0.5970725556,-0.114876192,-0.587862363,0.1709387013,0.1432406844,0.1120397432,1.1910119272,-0.273779871,-1.218254929,-1.022539935,-1.032580936,-0.094278468,-0.99540562,-0.747075916,-0.489838127,-0.094278468,-0.306760222,0.400666413,0.6374044247,-0.319726978,0.8409821065,0.6302148788,0.8628529351,1.0329640838,1.0705305723,0.8597633446,1.0883014454,0.1311700424,0.3818851749,-0.287978984,-0.264389617,-0.545175489,-0.536308688,-0.976624382,-0.715286637,0.1311700424,0.1523367091,-0.058430519,0.4119559144,0.8075155734,0.3818851749,0.1711179472,-0.038941106,-0.545175489,-0.536308688,-0.51752745,-1.166183658,-2.123315061,-1.454502551,-1.894818245,-2.067977699,1.4838611045,1.3000790381,2.2370541394,1.7646469765) sas.pea <- c(1.4619609019,0.3612064332,1.4379060694,1.5888531685,-0.48192063,-0.745939722,-1.219244704,-0.355028363,0.1660398808,-0.745939722,-0.776386242,0.0769453106,1.245974065,1.4683525887,0.3307599139,-0.139041526,-1.12988114,-0.30308126,-1.440673935,-0.5710152,0.8140003914,0.5826356643,0.552189145,0.5089189843,-0.265933793,-0.967368953,-0.55495701,-0.355028363,0.5980135545,-1.188798184,-0.55495701,-0.787002037,0.1660398808,-1.631656647,2.3236229938,-0.5710152,2.1099214124,1.6897818198,2.3236229938,1.8048400054,0.1660398808,-0.524510491,-0.997815473,-1.002988874,-0.48192063,-0.30308126,-0.776386242,0.2929321475,-2.425802161,0.1397772021,0.1093306829,0.940892658,-0.049946956,0.5826356643,-0.112098548,-0.5710152,0.1660398808,0.1397772021,0.1093306829,1.1568794948,-0.265933793,-1.188798184,-0.997815473,-1.002988874,-0.090321768,-0.960197666,-0.720651498,-0.469280491,-0.090321768,-0.295909972,0.3864946579,0.6106536933,-0.306308605,0.811236183,0.607923889,0.8266405301,0.9896124161,1.0326654141,0.8293531201,1.0426273669,0.1256650687,0.3683777208,-0.277793035,-0.253293654,-0.522295442,-0.517339204,-0.942080729,-0.685267328,0.1256650687,0.1469484897,-0.056363804,0.3946668564,0.7736255792,0.3683777208,0.1650654268,-0.037306817,-0.522295442,-0.517339204,-0.499222266,-1.117241001,-2.0342033,-1.403056128,-1.827797653,-1.981188349,1.4215860898,1.2540946452,2.1579285066,1.6905878775) checkEquals(round(as.numeric(resid(fit.R, "marginal", "raw")), 6), round(sas.raw, 6)) checkEquals(round(as.numeric(resid(fit.R, "marginal", "stud")), 10), sas.stu) checkEquals(round(as.numeric(resid(fit.R, "marginal", "pear")), 10), sas.pea) } # test conditional residuals against SAS PROC MIXED implementation # # Here the Orthodont-dataset from package 'nlme' is used (if available) TF003.conditional_residuals <- function() { data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho$Subject <- factor(as.character(Ortho$Subject)) fit.R <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, order.data=FALSE) sas.raw <- c(1.0554503692,-1.604344233,0.7358611644,1.0760665621,0.2894545536,-0.274141978,-1.33773851,0.5986649577,0.9931804233,-1.056673527,-1.106527478,0.843618572,0.9136978046,1.6799232283,-0.553851348,-1.287625924,-0.816278005,1.0788977187,-1.525926558,0.3692491659,0.5475924658,0.0389662598,0.0303400538,0.0217138478,0.4776426009,-1.099696512,-0.177035626,0.2456252606,2.0163981647,-1.827317063,-0.171032292,-0.51474752,0.4030450154,-3.770492168,5.0559706477,-1.617566536,0.8392188753,-0.210635075,1.2395109744,0.1896570239,1.1967876368,-0.119442428,-0.935672494,-0.751902559,-0.300680854,0.0120393803,-1.175240385,1.1374798497,-4.01735371,1.2731484286,0.5636505667,1.8541527049,-0.246387468,1.3274418148,-0.098728902,-1.024899619,-0.138123401,-0.394116074,-0.650108746,1.5938985808,0.9363555284,-1.00355777,-0.443471069,-0.383384367,0.8329387278,-1.068683204,-0.470305136,0.1280729319,-0.257137352,-0.892383169,0.4723710135,0.8371251962,-1.380761257,0.9565077625,0.2937767823,0.6310458021,0.3127531847,0.3561609268,-0.100431331,0.4429764109,0.142397156,0.79574555,-0.550906056,-0.397557662,0.0545559443,0.1529340124,-0.74868792,-0.150309852,0.0367239257,0.0389039233,-0.458916079,0.5432639185,0.8950900173,0.1034087373,-0.188272543,-0.479953823,-0.027899545,0.1254488494,0.2787972434,-1.067854363,-1.056621665,0.5005286586,-0.442321018,-0.885170694,0.1479608631,-0.418572047,1.5148950435,0.4483621337) sas.stu <- c(1.0130072422,-1.40421926,0.6440702679,1.0327943903,0.277814635,-0.239945666,-1.170869781,0.5745906729,0.9532413755,-0.924864682,-0.968499879,0.8096938975,0.876955014,1.4703705752,-0.484764251,-1.235846255,-0.783452784,0.9443166405,-1.33558336,0.3544004438,0.5255719736,0.034105631,0.026555453,0.0208406626,0.4584350227,-0.962521005,-0.154952304,0.2357478619,1.9353121698,-1.599378588,-0.149697822,-0.494047831,0.3868372511,-3.300163153,4.4252917896,-1.552518871,0.8054711271,-0.184360578,1.0848950914,0.1820302919,1.148660874,-0.104543249,-0.818957247,-0.721666086,-0.288589489,0.010537595,-1.028641578,1.0917380478,-3.855802718,1.114336629,0.4933411206,1.7795911329,-0.236479419,1.1618574896,-0.086413516,-0.983685038,-0.132569005,-0.344954262,-0.569014557,1.5298026822,0.8987015962,-0.878374554,-0.388152744,-0.367967222,0.8058885013,-0.936507499,-0.412137372,0.1239136802,-0.248786649,-0.782012412,0.4139477394,0.8099390114,-1.335920138,0.8382060174,0.2574422042,0.6105521795,0.3025963217,0.3121106213,-0.088009893,0.4285904641,0.1377727158,0.6973270207,-0.482769497,-0.384646719,0.0527842046,0.1340189954,-0.656089521,-0.145428441,0.0355312923,0.0340922509,-0.402156923,0.5256210697,0.8660213873,0.0906190511,-0.164986824,-0.464367011,-0.026993489,0.1099332223,0.2443153482,-1.033175098,-1.02230719,0.4386228215,-0.387614354,-0.856424201,0.1431557379,-0.366802677,1.3275314545,0.433801282) sas.pea <- c(0.805662962,-1.224653252,0.561709108,0.8214000384,0.2209509986,-0.209262363,-1.021143582,0.4569823435,0.7581300884,-0.806596642,-0.844651943,0.6439641857,0.6974581668,1.2823453984,-0.422774515,-0.982890855,-0.623094154,0.8235611613,-1.16479424,0.281861076,0.4179968863,0.0297443378,0.0231596466,0.0165749555,0.3646016562,-0.839437623,-0.135137616,0.1874945339,1.5391887345,-1.394856377,-0.130555057,-0.392925166,0.3076586549,-2.878151335,3.8594029686,-1.234746308,0.6406057401,-0.160785276,0.9461629957,0.1447719799,0.9135507463,-0.091174671,-0.714232232,-0.573954077,-0.229520435,0.009190089,-0.897102959,0.8682789942,-3.066589566,0.9718396658,0.4302546083,1.4153409804,-0.188076354,1.0132837466,-0.075363297,-0.782342484,-0.105434525,-0.300842875,-0.496251225,1.2166797125,0.7147536166,-0.766051488,-0.338517305,-0.29265098,0.6358118792,-0.815764058,-0.359000707,0.097762643,-0.196282124,-0.681187944,0.3605776652,0.6390075601,-1.053984381,0.7301365366,0.2242503101,0.4816997984,0.2387356758,0.2718703556,-0.076662822,0.3381397153,0.1086968349,0.6074209982,-0.420526268,-0.303469962,0.0416445008,0.1167399937,-0.571500228,-0.114736878,0.028032684,0.0296967541,-0.350307033,0.4146927516,0.6832541782,0.0789355824,-0.143715156,-0.366365894,-0.021296719,0.0957595871,0.2128158931,-0.815131373,-0.806557054,0.3820714014,-0.33763943,-0.675682405,0.1129438112,-0.319510993,1.1563734909,0.3422508298) checkEquals(round(as.numeric(resid(fit.R, "conditional", "raw")), 10), sas.raw) checkEquals(round(as.numeric(resid(fit.R, "conditional", "stud")), 10), sas.stu) checkEquals(round(as.numeric(resid(fit.R, "conditional", "pear")), 10), sas.pea) } # test different DDFM-methods for test of linear hypotheses via 'test.fixef' or 'test.lsmeans' # # ddfm="satterthwaite" can only be tested when using the same variance-convariance matrix of VCs # from SAS PROC MIXED --> impute this for testing TF004.ddfm.LMM.unbalanced <- function() { data(dataEP05A2_2) dat2ub <- dataEP05A2_2[-c(11,12,23,32,40,41,42),] fit2ub <- anovaMM(y~day/(run), dat2ub, VarVC.method="scm") L <- getL(fit2ub, c("day1-day2", "day2-day3", "day3-day6", "day14-day20"), "fixef") # ddfm="contain" (default) tst.con <- test.fixef(fit2ub, L=L, ddfm="contain") checkEquals(as.numeric(tst.con[,"DF"]), c(18, 18, 18, 18)) # ddfm="residual" tst.res <- test.fixef(fit2ub, L=L, ddfm="residual") checkEquals(as.numeric(tst.res[,"DF"]), c(53, 53, 53, 53)) # ddfm="satterthwaite" tst.satt <- test.fixef(fit2ub, L=L, ddfm="satt") checkEquals(as.numeric(tst.satt[,"DF"]), c(17.83, 17.83, 19.6, 17.83), tolerance=1) # relax equivalence threshold } # test different DDFM-methods for test of linear hypotheses via 'test.fixef' or 'test.lsmeans' # # ddfm="satterthwaite" can only be tested when using the same variance-convariance matrix of VCs # from SAS PROC MIXED --> impute this for testing # # Note: I used the variance-covariance matrix of variance components as displayed in SAS V9.3 64-Bit. # I could not figure out (in a reasonable amount of time) how to get this with full precision. # Thus, numeric differences between my implementation and that of SAS PROC MIXED are larger due # to error propagation. The numeric tolerance has been chosen larger than usual. TF005.ddfm.LMM.unbalanced <- function() { data(VCAdata1) datS2 <- VCAdata1[VCAdata1$sample == 2, ] datS2ub <- datS2[-c(15, 32, 33, 60, 62, 63, 64, 65, 74),] fitS2ub <- anovaMM(y~(lot+device)/(day)/(run), datS2ub) L <- getL(fitS2ub, c("lot1-lot2", "device1-device2"), "fixef") # ddfm="contain" (default) tst.con <- test.fixef(fitS2ub, L=L, ddfm="contain") checkEquals(as.numeric(tst.con[,"DF"]), c(58, 58)) # ddfm="residual" tst.res <- test.fixef(fitS2ub, L=L, ddfm="residual") checkEquals(as.numeric(tst.res[,"DF"]), c(238, 238)) # ddfm="satterthwaite" fit <- fitS2ub fit$VarCov <- matrix(c(0.004773, -0.00009, -5.14E-6, -0.00009, 0.000211, -0.00002, -5.14E-6, -0.00002, 0.000035), 3, 3) # extracted from SAS PROC MIXED output using options "asycov" in the proc mixed statement tst.satt <- test.fixef(fit, L=L, ddfm="satt") checkEquals(as.numeric(round(tst.satt[,"DF"], c(1,1))), c(56.4, 55.9), tolerance=0.1) # SAS-values extracted not sufficiently precise --> increase tolerance } # Test function 'stepwiseVCA' against hand-calculated values to verify the algorithm. TF006.stepwiseVCA.fully_nested <- function() { data(VCAdata1) datS7L1 <- VCAdata1[VCAdata1$sample == 7 & VCAdata1$lot == 1, ] fit0 <- anovaVCA(y~device/day/run, datS7L1) Ci <- getMat(fit0, "Ci.MS") tab <- fit0$aov.tab sw.res <- stepwiseVCA(fit0) nr <- nrow(fit0$aov.tab) for(i in 1:length(sw.res)) { checkEquals(sum(fit0$aov.tab[(nr-i):nr, "VC"]), sw.res[[i]]$aov.tab[1,"VC"]) # total variance correct? tmpTotDF <- VCA:::SattDF(tab[(nr-i):nr, "MS"], Ci[(nr-i-1):(nr-1), (nr-i-1):(nr-1)], tab[(nr-i):nr, "DF"], "total") checkEquals(tmpTotDF, sw.res[[i]]$aov.tab[1,"DF"]) # total DFs correct? } } # Tests for function 'getL'. TF007.getL.simple_contrasts <- function() { data(VCAdata1) dat1 <- VCAdata1[VCAdata1$sample==1,] fit <- anovaMM(y~(lot+device)/day/(run), dat1) L <- getL(fit, c("lot1-lot2", "1lot1-1*lot2")) checkEquals(L[1,], L[2,]) tmp <- rep(0,70) tmp[2:3] <- c(1,-1) checkEquals(as.numeric(L[1,]), tmp) } TF008.getL.complex_contrasts <- function() { data(VCAdata1) dat1 <- VCAdata1[VCAdata1$sample==1,] fit <- anovaMM(y~(lot+device)/day/run, dat1) L <- getL(fit, c("lot1:device1:day1-lot1:device1:day2", "lot2:device1:day2-lot3:device1:day4")) L0 <- matrix(0, nrow=2, ncol=ncol(L)) colnames(L0) <- colnames(L) rownames(L0) <- rownames(L) L0[1,c("lot1:device1:day1", "lot1:device1:day2")] <- c(1, -1) L0[2,c("lot2:device1:day2", "lot3:device1:day4")] <- c(1, -1) checkEquals(L, L0) } TF009.getL.complex_contrasts <- function() { data(VCAdata1) dat1 <- VCAdata1[VCAdata1$sample==1,] fit <- anovaMM(y~(lot+device)/day/run, dat1) L <- getL(fit, c("Custom Linear Hypothesis"="0.25lot1-.75*lot2")) checkEquals(as.numeric(L[2:3]), c(.25, -.75)) } # test function ranef and implicitely all intermediate results dataEP05A2_1 balanced (complete) # # SAS-reference (SAS 9.3 TS Level 1M2 X64_7PRO platform) results obtained using following SAS-code: # # proc mixed data=ep5_1 method=type1; # class day run; # model y =; # random day day*run/G V solution; # ods output G=G CovParms=CovParms V=V SolutionR=SolutionR; # run; TF010.ranef.balanced.nested <- function() { old.opt <- options("scipen" = 21) SASre <- c(0.000201,0.002605,-0.00054,-0.00182,-0.00119,0.000152,0.001162,-0.00122, 0.000293,-0.00087,-0.00009,0.000457,-0.00079,-0.00031,0.001124,0.00062,0.00058,0.000421,-0.00096,0.000188,0.23,1.7866,0.4298,-2.2927,-0.5359,0.3628,0.1836,-1.8719,0.7564,-1.6238,-0.6192,0.7676,-0.3797,-0.7438,-0.09473,0.6637,0.1288,0.5929,-1.4274,0.5381,0.01714,1.4101,-1.0957,0.05366,-0.9255,-0.1757,1.242,0.3798,-0.3971,0.5512,0.5044,-0.2064,-0.5919,0.3678,1.4737,0.09764,0.5825,-0.07603,0.245,-0.3075) digits <- nchar(gsub(".*\\.", "", as.character(SASre))) fit <- anovaVCA(y~day/run, dataEP05A2_1) fit <- solveMME(fit) checkEquals(as.numeric(round(ranef(fit)[,1], digits)), SASre) options(old.opt) } # test function ranef and implicitely all intermediate results dataEP05A2_2 balanced (complete) TF011.ranef.balanced.nested <- function() { old.opt <- options("scipen" = 21) SASre <- c(-1.3614,0.9946,-0.3492,2.2449,-1.0212,-0.2419,-0.09628,-1.1657,1.4828,-0.08096,-0.4642,-0.473,-0.2088,0.4715,1.0583,-0.17,-0.07924,-0.1185,0.4657,-0.8874,-1.536,-0.3423,-0.06817,0.5967,0.7873,-0.3391,-0.3554,0.1041,0.07095,-0.03257,-2.3438,1.1429,-0.4562,0.2562,1.6627,-0.9192,-1.2352,-0.6738,0.7737,-0.4168,-0.5393,1.8585,-0.4641,2.8256,-2.3442,-0.02973,0.2087,-1.8812,2.1895,-0.09085,1.6362,-1.864,0.1378,0.4626,-0.04928,0.6601,1.1144,0.4931,-0.06373,-0.9361) digits <- nchar(gsub(".*\\.", "", as.character(SASre))) fit <- anovaVCA(y~day/run, dataEP05A2_2) fit <- solveMME(fit) checkEquals(as.numeric(round(ranef(fit)[,1], digits)), SASre) options(old.opt) } # test function ranef and implicitely all intermediate results dataEP05A2_3 balanced (complete) TF012.ranef.balanced.nested <- function() { old.opt <- options("scipen" = 21) SASre <- c(0.6753,-2.1525,0.595,4.8068,1.9164,1.6353,-2.7515,-3.2263,0.2208,1.6819,0.005058,4.8909,0.9292,3.971,-0.2123,-5.0957,-3.8365,-1.7749,-1.7071,-0.5711,0.5098,-0.07488,1.0364,3.7443,0.2138,-1.9382,-2.065,-2.4335,-0.3697,0.008577,0.4418,2.8435,-1.7807,2.4792,2.1879,-0.1801,-1.2571,-0.6198,0.1537,-0.6488,-0.1216,-1.1625,-0.6944,-0.981,0.8879,2.8783,0.4832,0.5788,0.4966,0.9583,-0.4389,-0.03192,2.3149,-0.1964,-2.31,-2.7492,-0.9484,-0.4006,-1.1351,0.3206) digits <- nchar(gsub(".*\\.", "", as.character(SASre))) fit <- anovaVCA(y~day/run, dataEP05A2_3) fit <- solveMME(fit) checkEquals(as.numeric(round(ranef(fit)[,1], digits)), SASre) options(old.opt) } f <- function(x) { nam <- unlist(strsplit(as.character(x[1]), "\\*")) if(x[3] == "_") nam <- paste(nam, sub(" ", "", x[2]), sep="") else nam <- paste(nam, gsub(" ", "", x[2:3]), sep="") nam <- paste(nam, collapse=":") return(nam) } # test function ranef and implicitely all intermediate results dataEP05A2_1 unbalanced # deleting observations c(3, 13, 21, 22, 50) TF013.ranef.unbalanced.nested <- function() { old.opt <- options("scipen" = 21) SASre <- c(-0.04937,-0.6243,0.1375,0.3915,0.2943,0.03382,-0.275,0.3004,-0.06473,0.2177,0.02878,-0.1045,0.07342,0.08029,-0.2658,-0.144,-0.1341,-0.09579,0.2393,-0.03934,0.2539,2.2287,0.3491,-2.3509,-0.7547,0.3425,-2.1535,0.813,-1.8442,-0.6803,0.8489,0.2013,-0.8418,0.04628,0.7644,0.1997,0.6611,-1.6523,0.5696,0.06753,1.8355,-1.244,-0.1982,-1.1615,-0.2202,1.4478,0.198,-0.3916,0.4271,0.493,-0.1683,-0.6793,0.3191,1.6842,0.1732,0.6735,-0.03746,0.09424,-0.3134) digits <- nchar(gsub(".*\\.", "", as.character(SASre))) fit <- anovaVCA(y~day/run, dataEP05A2_1[-c(3,13,21,22,50), ], NegVC=TRUE) fit <- solveMME(fit) checkEquals(as.numeric(round(ranef(fit)[,1], digits)), SASre) options(old.opt) } # test function ranef and implicitely all intermediate results dataEP05A2_2 unbalanced # deleting observations c(8, 9, 12, 32, 35, 51, 52, 53, 67, 68, 73) TF014.ranef.unbalanced.nested <- function() { old.opt <- options("scipen" = 21) SASre <- c(-1.1456,0.6634,0.463,1.9848,-0.8504,-0.1739,-0.04745,-0.9272,0.9797,-0.03415,-0.3668,-0.3744,-0.2066,0.329,0.9548,-0.1114,-0.4822,-0.06678,0.1464,-0.7342,-1.6664,-0.08874,1.0566,0.8392,0.7674,-0.3348,-0.3397,0.01771,0.4464,-0.00412,-2.4289,1.1808,-0.4149,0.03393,1.8452,-0.9296,-0.9685,-0.6712,0.1018,-0.4684,-0.6344,1.4211,-0.1266,3.1471,-2.4752,-0.01445,0.2444,-1.8798,1.5212,-0.06447,1.6922,-1.9328,0.6268,0.07245,0.7058,0.5371,0.1922,-1.0062) digits <- nchar(gsub(".*\\.", "", as.character(SASre))) fit <- anovaVCA(y~day/run, dataEP05A2_2[-c(8, 9, 12, 32, 35, 51, 52, 53, 67, 68, 73), ], NegVC=TRUE) fit <- solveMME(fit) checkEquals(as.numeric(round(ranef(fit)[,1], digits)), SASre) options(old.opt) } # test function ranef and implicitely all intermediate results dataEP05A2_2 unbalanced # deleting observations c(1, 11, 21, 31, 41, 51, 61, 71) TF015.ranef.unbalanced.nested <- function() { old.opt <- options("scipen" = 21) SASre <- c(0.6958,-2.1514,0.3809,4.7516,1.8846,1.4953,-2.7456,-3.3254,0.2027,1.652,0.1019,4.8351,0.9999,3.9225,-0.2269,-5.2604,-3.8218,-1.0984,-1.7096,-0.5828,0.5248,-0.08342,1.068,3.5725,0.2055,-1.9519,-1.9777,-2.2817,-0.3553,0.009593,0.5286,2.7164,-1.7359,2.3667,2.0749,-0.3834,-1.2134,-0.899,0.1355,-0.6235,-0.1364,-1.1175,-0.8554,-0.9202,0.8465,2.7865,0.4451,0.4255,0.4684,0.9125,-0.4717,-0.01747,2.2941,-0.1771,-2.2015,-2.553,-0.9199,0.2859,-1.0898,0.2982) digits <- nchar(gsub(".*\\.", "", as.character(SASre))) fit <- anovaVCA(y~day/run, dataEP05A2_3[-c(1, 11, 21, 31, 41, 51, 61, 71 ), ], NegVC=TRUE) fit <- solveMME(fit) checkEquals(as.numeric(round(ranef(fit)[,1], digits)), SASre) options(old.opt) } TF016.as.matrix.anovaVCA <- function() { data(dataEP05A2_2) fit <- anovaVCA(y~day/run, dataEP05A2_2) checkEquals(c(as.matrix(fit)), c(fit$aov.tab)) } TF017.as.matrix.anovaMM <- function() { data(dataEP05A2_2) fit <- anovaMM(y~day/(run), dataEP05A2_2) checkEquals(c(as.matrix(fit)), c(fit$aov.tab)) } #' Test generic Satterthwaite approximation for total DF against Neter-Wasserman formula for fully-nested designs. #' #' @param input (data.frame) with components 'VC' (variance component) estimates and 'N' (number of levels per factor level). #' Row 1 should include the residual error level (repeatability), i.e. in the 'N' comonent the number of replicates per factor level #' of the next higher order. NWformula <- function(input=NULL) { if( is.null(input) ) stop("'input' is NULL, cannot compute Satterthwaite approximation of total degrees of freedom!") out <- input out$DF <- NA out$M <- NA out$a <- NA # compute all terms required for Neter & Wassermann formula term1 <- 1 # used for computing 'a'-terms for(i in 1:nrow(out)) { # compute DF per level out$DF[i] <- out$N[i]-1 for(j in (i+1):nrow(out)) { if(j > nrow(out)) break else out$DF[i] <- out$DF[i] * out$N[j] } # compute 'M'-term if(i == 1) out$M[i] <- out$VC[i] else { n <- 1 for(j in 1:(i-1)) n <- n * out$N[j] # i=2: n=n1 , i=3: n=n1*n2 out$M[i] <- out$VC[i]*n + out$M[i-1] # i=2: M2 = n1*VC2+M1 i=3: M3=n1*n2*VC3 + M2 } # compute 'a'-term if(i < nrow(out)) { out$a[i] <- term1 * (1-1/out$N[i]) term1 <- term1*(1/out$N[i]) } else out$a[i] <- term1 # second term for computing coefficient a_i is not used } DFtotal <- eval(parse(text=paste(paste(out$N, collapse="*"), "-1", sep=""))) # total residual DF numer <- 0 denom <- 0 for(i in 1:nrow(out)) { numer <- numer + out$a[i] * out$M[i] denom <- denom + (out$a[i] * out$M[i])^2 / out$DF[i] } DFsatt <- numer^2 / denom return(DFsatt) } data(VCAdata1) datS1 <- VCAdata1[VCAdata1$sample==1,] TF018.SattDF_total_vs_Neter_Wasserman <- function() { fit <- anovaVCA(y~day/run, datS1) input <- data.frame(VC=rev(fit$aov.tab[-1,"VC"]), N=c(6,2,21)) NWdf <- NWformula(input) checkEquals(NWdf, fit$aov.tab[1,"DF"]) } # deeper nesting-structure TF019.SattDF_total_vs_Neter_Wasserman <- function() { fit <- anovaVCA(y~device/day/run, datS1) input <- data.frame(VC=rev(fit$aov.tab[-1,"VC"]), N=c(6,2,7,3)) NWdf <- NWformula(input) checkEquals(NWdf, fit$aov.tab[1,"DF"]) } # deepest nesting structure, use slightly larger precision, due to # the deep nesting structure resulting in many multiplications/division # which results in more severe error propagation TF020.SattDF_total_vs_Neter_Wasserman <- function() { fit <- anovaVCA(y~lot/device/day/run, datS1) input <- data.frame(VC=rev(fit$aov.tab[-1,"VC"]), N=c(2,2,7,3,3)) NWdf <- NWformula(input) checkEquals(NWdf, fit$aov.tab[1,"DF"], tol=1e-7) } # test compuation of LS Means against SAS PROX MIXED # specifically test: LS Means estimates, standard errors, # DFs (method="containment"), t-Value # # proc mixed data=vcaS1 method=type1; # class lot device day run; # model y= lot device; # random lot*device*day lot*device*day*run; # lsmeans lot; # lsmeans device; # run; TF021.lsmeans.balanced <- function() { data(VCAdata1) datS1 <- VCAdata1[VCAdata1$sample == 1, ] fit <- anovaMM(y~(lot+device)/(day)/(run), datS1) lsm <- lsmeans(fit, type="c", ddfm="cont") checkEquals(round(as.numeric(lsm[,"Estimate"]), 4), c(2.8138, 2.5507, 2.6719, 2.7887, 2.7136, 2.5341)) checkEquals(round(as.numeric(lsm[,"SE"]), 5), rep(0.04266, 6)) checkEquals(round(as.numeric(lsm[,"DF"]), 5), rep(58, 6)) checkEquals(round(as.numeric(lsm[,"t Value"]), 2), c(65.96, 59.79, 62.64, 65.37, 63.61, 59.41)) } TF022.lsmeans.unbalanced <- function() { data(VCAdata1) datS1 <- VCAdata1[VCAdata1$sample == 1, ] datS1ub <- datS1[-c(1,11,12,20,55,56,57,103,121,122,179),] fit <- anovaMM(y~(lot+device)/(day)/(run), datS1ub) lsm <- lsmeans(fit) checkEquals(round(as.numeric(lsm[,"Estimate"]), 4), c(2.8062, 2.5549, 2.6718, 2.7778, 2.7207, 2.5343)) checkEquals(round(as.numeric(lsm[,"SE"]), 5), c(0.04106, 0.04062, 0.04019, 0.04063, 0.04105, 0.04019)) } # test different DDFM-methods for test of linear hypotheses via 'test.fixef' or 'test.lsmeans' # # this model includes (numeric) covariate # # ddfm="satterthwaite" can only be tested when using the same variance-convariance matrix of VCs # from SAS PROC MIXED --> impute this for testing TF023.ddfm_fixef.all_methods.Orthodont.balanced <- function() { data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 fit.Ortho <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho) L <- getL(fit.Ortho, "SexMale-SexFemale", "fixef") # ddfm="contain" (default) tst.con <- test.fixef(fit.Ortho, L=L, ddfm="contain") checkEquals(as.numeric(tst.con[,"DF"]), 54) # ddfm="residual" tst.res <- test.fixef(fit.Ortho, L=L, ddfm="residual") checkEquals(as.numeric(tst.res[,"DF"]), 104) # ddfm="satterthwaite" fit <- fit.Ortho fit$VarCov <- matrix(c(1.1494, 0.001364, -0.02727, 0.001364, 0.001393, -0.00545, -0.02727, -0.00545, 0.1091), 3, 3) tst.satt <- test.fixef(fit, L=L, ddfm="satt") checkEquals(as.numeric(tst.satt[,"DF"]), 25, tolerance=0.1) } TF024.ddfm_lsmeans.all_methods.Orthodont.unbalanced <- function() { data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho.ub <- Ortho[-c(5,7,23,33,51,72,73,74,90),] fit.Ortho <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho.ub) L <- getL(fit.Ortho, "SexMale-SexFemale", "lsmeans") # ddfm="contain" (default) tst.con <- test.lsmeans(fit.Ortho, L=L, ddfm="contain") checkEquals(as.numeric(tst.con[,"DF"]), 45) # ddfm="residual" tst.res <- test.lsmeans(fit.Ortho, L=L, ddfm="residual") checkEquals(as.numeric(tst.res[,"DF"]), 95) # ddfm="satterthwaite" fit <- fit.Ortho fit$VarCov <- matrix(c(1.4381, 0.002238, -0.05008, 0.002238, 0.001532, -0.00800, -0.05008, -0.00800, 0.1572), 3, 3) tst.satt <- test.lsmeans(fit, L=L, ddfm="satt") checkEquals(round(as.numeric(tst.satt[,"DF"]),2), 23.35, tolerance=0.1) # allow larger numeric tolerance due to systematic difference between R and SAS implementation # of the Satterthwaite approximation of the denominator degrees of freedom } # proc mixed data=dats1cov method=type1; # class device day run; # model y=device day*cov/solution; # random device*run; # lsmeans device; # run; TF025.lsmeans.including_covariate <- function() { data(VCAdata1) datS1 <- VCAdata1[VCAdata1$sample == 1, ] set.seed(20140608) datS1$cov <- 20 + rnorm(252) fit <- anovaMM(y~device+day:cov+device:(run), datS1) lsm <- lsmeans(fit, "device") checkEquals( round(as.numeric(lsm[,"Estimate"]), 4), c(2.8773, 1.8635, 3.2976) ) } # proc mixed data=datS1cov method=type1; # class lot device day run; # model y = lot device day*cov lot*device*day/solution; # random lot*device*day*run; # lsmeans lot; # lsmeans device; # run; TF026.lsmeans.complex_model.including_covariate <- function() { data(VCAdata1) datS1 <- VCAdata1[VCAdata1$sample == 1, ] set.seed(20140608) datS1$cov <- 20 + rnorm(252) fit <- anovaMM(y~lot+device+day:cov+lot:device:day+lot:device:day:(run), datS1) lsm <- lsmeans(fit, c("lot", "device")) checkEquals( round(as.numeric(lsm[,"Estimate"]), 4), c(2.8143, 2.5496, 2.6720, 2.7762, 2.7380, 2.5217) ) } # check whether result of by-processed samples are identical to # separately computed results. TF027.anovaVCA.by_processing <- function() { data(CA19_9) fit.lst <- anovaVCA(result~site/day, CA19_9, by="sample") samples <- gsub("sample\\.", "",names(fit.lst)) for(i in 1:length(fit.lst)) { tmp.fit <- anovaVCA(result~site/day, CA19_9[CA19_9$sample == samples[i],]) checkEquals(fit.lst[[i]]$aov.tab, tmp.fit$aov.tab) } } # check whether result of by-processed samples are identical to # separately computed results. Use artifical mixed model with site # as fixed factor TF028.anovaMM.by_processing <- function() { data(CA19_9) fit.lst <- anovaMM(result~site/(day), CA19_9, by="sample") samples <- gsub("sample\\.", "",names(fit.lst)) for(i in 1:length(fit.lst)) { tmp.fit <- anovaMM(result~site/(day), CA19_9[CA19_9$sample == samples[i],]) print(checkEquals(fit.lst[[i]]$aov.tab, tmp.fit$aov.tab)) } } # check whether function VCAinference correctly handles # the list-object returned by function 'anovaVCA' when by-processing # is used TF029.anovaVCA.by_processing <- function() { data(CA19_9) fit.lst <- anovaVCA(result~site/day, CA19_9, by="sample") samples <- gsub("sample\\.", "",names(fit.lst)) inf.lst <- VCAinference(fit.lst) for(i in 1:length(fit.lst)) { tmp.fit <- anovaVCA(result~site/day, CA19_9[CA19_9$sample == samples[i],]) tmp.inf <- VCAinference(tmp.fit) checkEquals(inf.lst[[i]]$VCAobj$aov.tab, tmp.inf$VCAobj$aov.tab) checkEquals(inf.lst[[i]]$ConfInt, tmp.inf$ConfInt) } } # check whether function VCAinference correctly handles # the list-object returned by function 'anovaVCA' when by-processing # and vector-arguments, e.g. TF030.anovaVCA.by_processing <- function() { data(CA19_9) fit.lst <- anovaVCA(result~site/day, CA19_9, by="sample") samples <- gsub("sample\\.", "",names(fit.lst)) total.specs <- c(1, 3, 5, 30, 80, 200) error.specs <- c(.5, 2, 2, 5, 50, 75) inf.lst <- VCAinference(fit.lst, total.claim=total.specs, error.claim=error.specs) for(i in 1:length(fit.lst)) { tmp.fit <- anovaVCA(result~site/day, CA19_9[CA19_9$sample == samples[i],]) tmp.inf <- VCAinference(tmp.fit, total.claim=total.specs[i], error.claim=error.specs[i]) checkEquals(inf.lst[[i]]$VCAobj$aov.tab, tmp.inf$VCAobj$aov.tab) } } # check whether linear hypothesis of LSMeans including constrained fixed effects does return # a value and does not give a warning any more. TF031.test.lsmeans <- function() { data(dataEP05A2_1) fit <- anovaMM(y~day/(run), dataEP05A2_1) lc.mat <- getL(fit, "day19-day20", "lsm") # day20 is contrained to 0 res <- test.lsmeans(fit, lc.mat) checkEquals(as.numeric(round(res, 6)), c(-1.220903, 20, 1.523546, -0.801356, 0.432343)) } # check equality of residuals for balanced models fitted once by ANOVA and once by REML TF032.ANOVA_vs_REML.residuals.raw <- function() { data(dataEP05A2_1) fit1.aov <- anovaVCA(y~day/run, dataEP05A2_1) fit1.reml <- remlVCA(y~day/run, dataEP05A2_1) checkEquals(round(resid(fit1.aov), 4), round(resid(fit1.reml), 4)) data(dataEP05A2_2) fit2.aov <- anovaVCA(y~day/run, dataEP05A2_2) fit2.reml <- remlVCA(y~day/run, dataEP05A2_2) checkEquals(round(resid(fit2.aov), 6), round(resid(fit2.reml), 6)) data(dataEP05A2_3) fit3.aov <- anovaVCA(y~day/run, dataEP05A2_3) fit3.reml <- remlVCA(y~day/run, dataEP05A2_3) checkEquals(round(resid(fit3.aov), 4), round(resid(fit3.reml), 4)) data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho$Subject <- factor(as.character(Ortho$Subject)) fit.anovaMM <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho) fit.remlMM <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, cov=FALSE) checkEquals(round(resid(fit.anovaMM), 4), round(resid(fit.remlMM), 4)) } TF033.ANOVA_vs_REML.residuals.studentized <- function() { data(dataEP05A2_1) fit1.aov <- anovaVCA(y~day/run, dataEP05A2_1) fit1.reml <- remlVCA(y~day/run, dataEP05A2_1) checkEquals(round(resid(fit1.aov, mode="student"), 4), round(resid(fit1.reml, mode="student"), 4)) data(dataEP05A2_2) fit2.aov <- anovaVCA(y~day/run, dataEP05A2_2) fit2.reml <- remlVCA(y~day/run, dataEP05A2_2) checkEquals(round(resid(fit2.aov, mode="student"), 5), round(resid(fit2.reml, mode="student"), 5)) data(dataEP05A2_3) fit3.aov <- anovaVCA(y~day/run, dataEP05A2_3) fit3.reml <- remlVCA(y~day/run, dataEP05A2_3) checkEquals(round(resid(fit3.aov, mode="student"), 4), round(resid(fit3.reml, mode="student"), 4)) data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho$Subject <- factor(as.character(Ortho$Subject)) fit.anovaMM <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho) fit.remlMM <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, cov=FALSE) checkEquals(round(resid(fit.anovaMM, mode="student"), 4), round(resid(fit.remlMM, mode="student"), 4)) } TF034.ANOVA_vs_REML.residuals.pearson <- function() { data(dataEP05A2_1) fit1.aov <- anovaVCA(y~day/run, dataEP05A2_1) fit1.reml <- remlVCA(y~day/run, dataEP05A2_1) checkEquals(round(resid(fit1.aov, mode="pearson"), 4), round(resid(fit1.reml, mode="pearson"), 4)) data(dataEP05A2_2) fit2.aov <- anovaVCA(y~day/run, dataEP05A2_2) fit2.reml <- remlVCA(y~day/run, dataEP05A2_2) checkEquals(round(resid(fit2.aov, mode="pearson"), 6), round(resid(fit2.reml, mode="pearson"), 6)) data(dataEP05A2_3) fit3.aov <- anovaVCA(y~day/run, dataEP05A2_3) fit3.reml <- remlVCA(y~day/run, dataEP05A2_3) checkEquals(round(resid(fit3.aov, mode="pearson"), 4), round(resid(fit3.reml, mode="pearson"), 4)) data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho$Subject <- factor(as.character(Ortho$Subject)) fit.anovaMM <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho) fit.remlMM <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, cov=FALSE) checkEquals(round(resid(fit.anovaMM, mode="pearson"), 4), round(resid(fit.remlMM, mode="pearson"), 4)) } TF035.ANOVA_vs_REML.residuals.standardized <- function() { data(dataEP05A2_1) fit1.aov <- anovaVCA(y~day/run, dataEP05A2_1) fit1.reml <- remlVCA(y~day/run, dataEP05A2_1) checkEquals(round(resid(fit1.aov, mode="standard"), 4), round(resid(fit1.reml, mode="standard"), 4)) data(dataEP05A2_2) fit2.aov <- anovaVCA(y~day/run, dataEP05A2_2) fit2.reml <- remlVCA(y~day/run, dataEP05A2_2) checkEquals(round(resid(fit2.aov, mode="standard"), 6), round(resid(fit2.reml, mode="standard"), 6)) data(dataEP05A2_3) fit3.aov <- anovaVCA(y~day/run, dataEP05A2_3) fit3.reml <- remlVCA(y~day/run, dataEP05A2_3) checkEquals(round(resid(fit3.aov, mode="standard"), 4), round(resid(fit3.reml, mode="standard"), 4)) data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho$Subject <- factor(as.character(Ortho$Subject)) fit.anovaMM <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho) fit.remlMM <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, cov=FALSE) checkEquals(round(resid(fit.anovaMM, mode="standard"), 4), round(resid(fit.remlMM, mode="standard"), 4)) } # check equality of random effects for balanced models fitted once by ANOVA and once by REML TF036.ANOVA_vs_REML.ranef.raw <- function() { data(dataEP05A2_1) fit1.aov <- anovaVCA(y~day/run, dataEP05A2_1) fit1.reml <- remlVCA(y~day/run, dataEP05A2_1) re.aov <- ranef(fit1.aov, mode="raw") re.reml <- ranef(fit1.reml, mode="raw") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(dataEP05A2_2) fit2.aov <- anovaVCA(y~day/run, dataEP05A2_2) fit2.reml <- remlVCA(y~day/run, dataEP05A2_2) re.aov <- ranef(fit2.aov, mode="raw") re.reml <- ranef(fit2.reml, mode="raw") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(dataEP05A2_3) fit3.aov <- anovaVCA(y~day/run, dataEP05A2_3) fit3.reml <- remlVCA(y~day/run, dataEP05A2_3) re.aov <- ranef(fit3.aov, mode="raw") re.reml <- ranef(fit3.reml, mode="raw") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho$Subject <- factor(as.character(Ortho$Subject)) fit.anovaMM <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho) fit.remlMM <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, cov=FALSE) re.aov <- ranef(fit.anovaMM, mode="raw") re.reml <- ranef(fit.remlMM, mode="raw") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) } TF037.ANOVA_vs_REML.ranef.studentized <- function() { data(dataEP05A2_1) fit1.aov <- anovaVCA(y~day/run, dataEP05A2_1) fit1.reml <- remlVCA(y~day/run, dataEP05A2_1) re.aov <- ranef(fit1.aov, mode="student") re.reml <- ranef(fit1.reml, mode="student") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(dataEP05A2_2) fit2.aov <- anovaVCA(y~day/run, dataEP05A2_2) fit2.reml <- remlVCA(y~day/run, dataEP05A2_2) re.aov <- ranef(fit2.aov, mode="student") re.reml <- ranef(fit2.reml, mode="student") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(dataEP05A2_3) fit3.aov <- anovaVCA(y~day/run, dataEP05A2_3) fit3.reml <- remlVCA(y~day/run, dataEP05A2_3) re.aov <- ranef(fit3.aov, mode="student") re.reml <- ranef(fit3.reml, mode="student") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho$Subject <- factor(as.character(Ortho$Subject)) fit.anovaMM <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho) fit.remlMM <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, cov=FALSE) re.aov <- ranef(fit.anovaMM, mode="student") re.reml <- ranef(fit.remlMM, mode="student") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) } TF038.ANOVA_vs_REML.ranef.standardized <- function() { data(dataEP05A2_1) fit1.aov <- anovaVCA(y~day/run, dataEP05A2_1) fit1.reml <- remlVCA(y~day/run, dataEP05A2_1) re.aov <- ranef(fit1.aov, mode="standard") re.reml <- ranef(fit1.reml, mode="standard") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(dataEP05A2_2) fit2.aov <- anovaVCA(y~day/run, dataEP05A2_2) fit2.reml <- remlVCA(y~day/run, dataEP05A2_2) re.aov <- ranef(fit2.aov, mode="standard") re.reml <- ranef(fit2.reml, mode="standard") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(dataEP05A2_3) fit3.aov <- anovaVCA(y~day/run, dataEP05A2_3) fit3.reml <- remlVCA(y~day/run, dataEP05A2_3) re.aov <- ranef(fit3.aov, mode="standard") re.reml <- ranef(fit3.reml, mode="standard") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age - 11 Ortho$Subject <- factor(as.character(Ortho$Subject)) fit.anovaMM <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho) fit.remlMM <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho, cov=FALSE) re.aov <- ranef(fit.anovaMM, mode="standard") re.reml <- ranef(fit.remlMM, mode="standard") re.reml <- re.reml[rownames(re.aov),,drop=FALSE] checkEquals(round(re.aov[,1], 4), round(re.reml[,1], 4)) } TF039.remlVCA.by_processing <- function() { data(CA19_9) fit.lst <- remlVCA(result~site/day, CA19_9, by="sample") samples <- gsub("sample\\.", "",names(fit.lst)) total.specs <- c(1, 3, 5, 30, 80, 200) error.specs <- c(.5, 2, 2, 5, 50, 75) inf.lst <- VCAinference(fit.lst, total.claim=total.specs, error.claim=error.specs) for(i in 1:length(fit.lst)) { tmp.fit <- remlVCA(result~site/day, CA19_9[CA19_9$sample == samples[i],]) tmp.inf <- VCAinference(tmp.fit, total.claim=total.specs[i], error.claim=error.specs[i]) checkEquals(inf.lst[[i]]$VCAobj$aov.tab, tmp.inf$VCAobj$aov.tab) checkEquals(inf.lst[[i]]$ConfInt, tmp.inf$ConfInt) } } TF040.remlMM.by_processing <- function() { data(CA19_9) fit.lst <- remlMM(result~(site)/day, CA19_9, by="sample") samples <- gsub("sample\\.", "",names(fit.lst)) total.specs <- c(1, 3, 5, 30, 80, 200) error.specs <- c(.5, 2, 2, 5, 50, 75) inf.lst <- VCAinference(fit.lst, total.claim=total.specs, error.claim=error.specs) for(i in 1:length(fit.lst)) { tmp.fit <- remlMM(result~(site)/day, CA19_9[CA19_9$sample == samples[i],]) tmp.inf <- VCAinference(tmp.fit, total.claim=total.specs[i], error.claim=error.specs[i]) print(checkEquals(inf.lst[[i]]$VCAobj$aov.tab, tmp.inf$VCAobj$aov.tab)) checkEquals(inf.lst[[i]]$ConfInt, tmp.inf$ConfInt) } } # check whether linear hypothesis of LSMeans including constrained fixed effects does return # a value and does not give a warning any more. TF041.REML.test.lsmeans <- function() { data(dataEP05A2_1) fit <- remlMM(y~day/(run), dataEP05A2_1) lc.mat <- getL(fit, "day19-day20", "lsm") # day20 is contrained to 0 res <- test.lsmeans(fit, lc.mat) checkEquals(round(as.numeric(res),7), c( -1.2209032, 20, 1.5235456, -0.8013565, 0.4323434)) } # check whether unordered data lead to equal results as ordered data regarding balancedness or not TF042.balancedness.ordered_vs_unordered <- function() { data(dataEP05A2_1) dat <- dataEP05A2_1 dat <- dat[sample(1:nrow(dat)),] # no fixed seed required, has to work with any permutation fit1 <- anovaVCA(y~day/run, dat) fit2 <- remlVCA(y~day/run, dat) fit3 <- anovaMM(y~day/(run), dat) fit4 <- remlMM(y~day/(run), dat) checkEquals(fit1$balanced, "balanced") checkEquals(fit2$balanced, "balanced") checkEquals(fit3$balanced, "balanced") checkEquals(fit4$balanced, "balanced") } # check LS Means evaluation at specific values of covariates and/or for different weighting of factor-variables #load(file=file.path(R.home(), "library/VCA/tests/runit/UnitTests/LSMeans_Data.RData")) data(LSMeans_Data) fit.vca <- remlMM(y~snp+time+snp:time+sex+(id)+(id):time, LSMeans_Data, VarVC=FALSE) # results differ from SAS PROC MIXED results after 2nd decimal place because covariance parameters # are slightly different as well TF043.LSMeans.atCovarLevel <- function() { lsm <- lsmeans(fit.vca, var="snp", at=list(time=1:4)) checkEquals(as.numeric(round(lsm[-c(1,2),"Estimate"],2)), c(4.89, 5.01, 5.80, 5.92, 6.71, 6.83, 7.62, 7.75)) } # test whether specifying a non-existing covariable leads to result with only two (original) # LS Means. TF044.LSMeans.atCovarLevel <- function() { lsm <- lsmeans(fit.vca, var="snp", at=list(tim=1:4)) checkEquals(nrow(lsm), 2) } # test whether specifying a weighting scheme where elements add to something !=1 leads to result # with only two (original) LS Means. TF045.LSMeans.atCovarLevel <- function() { lsm1 <- lsmeans(fit.vca, var="snp", at=list(sex=c(Male=.3, Female=.6))) checkEquals(nrow(lsm1), 2) lsm2 <- lsmeans(fit.vca, var="snp", at=list(sex=c(Male=.5, Female=.6))) checkEquals(nrow(lsm2), 2) } # test whether anovaMM and remlMM only accept correctly specified terms in the model formula TF046.model.terms <- function() { data(dataEP05A2_1) checkException(anovaMM(y~(a)+.b:y, dataEP05A2_1)) checkException(remlMM(y~(a)+.b:y, dataEP05A2_1)) fit0 <- anovaVCA(y~day/run, dataEP05A2_1) dat <- dataEP05A2_1 dat$day.var <- dat$day checkEquals(as.numeric(anovaVCA(y~day.var/run, dat)$aov.tab[,"VC"]), as.numeric(fit0$aov.tab[,"VC"]), tolerance=1e-8) checkEquals(as.numeric(anovaMM(y~(day.var)/(run), dat)$aov.tab[,"VC"]), as.numeric(fit0$aov.tab[,"VC"]), tolerance=1e-8) checkEquals(as.numeric(remlMM(y~(day.var)/(run), dat)$aov.tab[,"VC"]), as.numeric(fit0$aov.tab[,"VC"]), tolerance=1e-7) checkEquals(as.numeric(remlVCA(y~day.var/run, dat)$aov.tab[,"VC"]), as.numeric(fit0$aov.tab[,"VC"]), tolerance=1e-7) } # test whether function as.matrix.VCA correctly works for models fitted via REML TF047.as.matrix.VCA.REML <- function() { data(dataEP05A2_1) fit1 <- remlVCA(y~day/run, dataEP05A2_1) mat1 <- as.matrix(fit1) checkEquals( fit1$aov.tab[,"DF"], as.numeric(mat1[,"DF"])) checkEquals( fit1$aov.tab[,"VC"], as.numeric(mat1[,"VC"])) checkEquals( fit1$aov.tab[,"%Total"], as.numeric(mat1[,"%Total"])) checkEquals( fit1$aov.tab[,"SD"], as.numeric(mat1[,"SD"])) checkEquals( fit1$aov.tab[,"CV[%]"], as.numeric(mat1[,"CV[%]"])) checkEquals( fit1$aov.tab[,"Var(VC)"], as.numeric(mat1[,"Var(VC)"])) fit2 <- remlMM(y~day/(run), dataEP05A2_1) mat2 <- as.matrix(fit2) checkEquals( fit2$aov.tab[,"DF"], as.numeric(mat2[,"DF"])) checkEquals( fit2$aov.tab[,"VC"], as.numeric(mat2[,"VC"])) checkEquals( fit2$aov.tab[,"%Total"], as.numeric(mat2[,"%Total"])) checkEquals( fit2$aov.tab[,"SD"], as.numeric(mat2[,"SD"])) checkEquals( fit2$aov.tab[,"CV[%]"], as.numeric(mat2[,"CV[%]"])) checkEquals( fit2$aov.tab[,"Var(VC)"], as.numeric(mat2[,"Var(VC)"])) } # test whether function as.matrix.VCAinference correctly works TF048.as.matrix.VCAinference.REML <- function() { data(dataEP05A2_2) fit.reml <- remlVCA(y~day/run, dataEP05A2_2) inf.reml <- VCAinference(fit.reml) VC.mat <- as.matrix(inf.reml, what="VC", digits=12) SD.mat <- as.matrix(inf.reml, what="SD", digits=12) CV.mat <- as.matrix(inf.reml, what="CV", digits=12) # check individual matrices checkEquals(as.numeric(fit.reml$aov.tab[,"VC"]), as.numeric(VC.mat[,"Estimate"])) checkEquals(as.numeric(inf.reml$ConfInt$VC$OneSided[,"LCL"]), as.numeric(VC.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.reml$ConfInt$VC$OneSided[,"UCL"]), as.numeric(VC.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.reml$ConfInt$VC$TwoSided[,"LCL"]), as.numeric(VC.mat[,"CI LCL"])) checkEquals(as.numeric(inf.reml$ConfInt$VC$TwoSided[,"UCL"]), as.numeric(VC.mat[,"CI UCL"])) checkEquals(as.numeric(fit.reml$aov.tab[,"SD"]), as.numeric(SD.mat[,"Estimate"])) checkEquals(as.numeric(inf.reml$ConfInt$SD$OneSided[,"LCL"]), as.numeric(SD.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.reml$ConfInt$SD$OneSided[,"UCL"]), as.numeric(SD.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.reml$ConfInt$SD$TwoSided[,"LCL"]), as.numeric(SD.mat[,"CI LCL"])) checkEquals(as.numeric(inf.reml$ConfInt$SD$TwoSided[,"UCL"]), as.numeric(SD.mat[,"CI UCL"])) checkEquals(as.numeric(fit.reml$aov.tab[,"CV[%]"]), as.numeric(CV.mat[,"Estimate"])) checkEquals(as.numeric(inf.reml$ConfInt$CV$OneSided[,"LCL"]), as.numeric(CV.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.reml$ConfInt$CV$OneSided[,"UCL"]), as.numeric(CV.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.reml$ConfInt$CV$TwoSided[,"LCL"]), as.numeric(CV.mat[,"CI LCL"])) checkEquals(as.numeric(inf.reml$ConfInt$CV$TwoSided[,"UCL"]), as.numeric(CV.mat[,"CI UCL"])) # check list of matrices mat.list <- as.matrix(inf.reml, digits=12) checkEquals(VC.mat, mat.list[[1]]) checkEquals(SD.mat, mat.list[[2]]) checkEquals(CV.mat, mat.list[[3]]) } TF049.as.matrix.VCAinference.ANOVA <- function() { data(dataEP05A2_3) fit.anova <- anovaVCA(y~day/run, dataEP05A2_3) inf.anova <- VCAinference(fit.anova, VarVC=FALSE) # there should be many NAs indicating missing values VC.mat <- as.matrix(inf.anova, what="VC", digits=12) SD.mat <- as.matrix(inf.anova, what="SD", digits=12) CV.mat <- as.matrix(inf.anova, what="CV", digits=12) # check individual matrices for missing values for indermediate variance components checkEquals(as.numeric(fit.anova$aov.tab[,"VC"]), as.numeric(VC.mat[,"Estimate"])) checkEquals(as.numeric(inf.anova$ConfInt$VC$OneSided[,"LCL"]), as.numeric(VC.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$VC$OneSided[,"UCL"]), as.numeric(VC.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.anova$ConfInt$VC$TwoSided[,"LCL"]), as.numeric(VC.mat[,"CI LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$VC$TwoSided[,"LCL"]), as.numeric(VC.mat[,"CI LCL"])) checkEquals(as.numeric(fit.anova$aov.tab[,"SD"]), as.numeric(SD.mat[,"Estimate"])) checkEquals(as.numeric(inf.anova$ConfInt$SD$OneSided[,"LCL"]), as.numeric(SD.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$SD$OneSided[,"UCL"]), as.numeric(SD.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.anova$ConfInt$SD$TwoSided[,"LCL"]), as.numeric(SD.mat[,"CI LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$SD$TwoSided[,"LCL"]), as.numeric(SD.mat[,"CI LCL"])) checkEquals(as.numeric(fit.anova$aov.tab[,"CV[%]"]), as.numeric(CV.mat[,"Estimate"])) checkEquals(as.numeric(inf.anova$ConfInt$CV$OneSided[,"LCL"]), as.numeric(CV.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$CV$OneSided[,"UCL"]), as.numeric(CV.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.anova$ConfInt$CV$TwoSided[,"LCL"]), as.numeric(CV.mat[,"CI LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$CV$TwoSided[,"LCL"]), as.numeric(CV.mat[,"CI LCL"])) # check for VarVC=TRUE inf.anova <- VCAinference(fit.anova, VarVC=TRUE) # there should be many NAs indicating missing values VC.mat <- as.matrix(inf.anova, what="VC", digits=12) SD.mat <- as.matrix(inf.anova, what="SD", digits=12) CV.mat <- as.matrix(inf.anova, what="CV", digits=12) # check individual matrices for missing values for indermediate variance components checkEquals(as.numeric(fit.anova$aov.tab[,"VC"]), as.numeric(VC.mat[,"Estimate"])) checkEquals(as.numeric(inf.anova$ConfInt$VC$OneSided[,"LCL"]), as.numeric(VC.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$VC$OneSided[,"UCL"]), as.numeric(VC.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.anova$ConfInt$VC$TwoSided[,"LCL"]), as.numeric(VC.mat[,"CI LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$VC$TwoSided[,"LCL"]), as.numeric(VC.mat[,"CI LCL"])) checkEquals(as.numeric(fit.anova$aov.tab[,"SD"]), as.numeric(SD.mat[,"Estimate"])) checkEquals(as.numeric(inf.anova$ConfInt$SD$OneSided[,"LCL"]), as.numeric(SD.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$SD$OneSided[,"UCL"]), as.numeric(SD.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.anova$ConfInt$SD$TwoSided[,"LCL"]), as.numeric(SD.mat[,"CI LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$SD$TwoSided[,"LCL"]), as.numeric(SD.mat[,"CI LCL"])) checkEquals(as.numeric(fit.anova$aov.tab[,"CV[%]"]), as.numeric(CV.mat[,"Estimate"])) checkEquals(as.numeric(inf.anova$ConfInt$CV$OneSided[,"LCL"]), as.numeric(CV.mat[,"One-Sided LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$CV$OneSided[,"UCL"]), as.numeric(CV.mat[,"One-Sided UCL"])) checkEquals(as.numeric(inf.anova$ConfInt$CV$TwoSided[,"LCL"]), as.numeric(CV.mat[,"CI LCL"])) checkEquals(as.numeric(inf.anova$ConfInt$CV$TwoSided[,"LCL"]), as.numeric(CV.mat[,"CI LCL"])) } TF050.lmerSummary <- function() { data(VCAdata1) fit0 <- remlVCA(y~(device+lot)/day/run, subset(VCAdata1, sample==5)) fit1 <- lme4:::lmer(y~(1|device)+(1|lot)+(1|device:lot:day)+(1|device:lot:day:run), subset(VCAdata1, sample==5), control=lme4:::lmerControl(optimizer="bobyqa")) sum1 <- lmerSummary(fit1, tab.only=TRUE) sum1 <- sum1[rownames(fit0$aov.tab),] checkEquals(as.numeric(fit0$aov.tab[,"VC"]), as.numeric(sum1[,"VC"]), tolerance=1e-5) checkEquals(as.numeric(fit0$aov.tab[,"%Total"]), as.numeric(sum1[,"%Total"]), tolerance=1e-5) checkEquals(as.numeric(fit0$aov.tab[,"SD"]), as.numeric(sum1[,"SD"]), tolerance=1e-6) checkEquals(as.numeric(fit0$aov.tab[,"CV[%]"]), as.numeric(sum1[,"CV[%]"]), tolerance=1e-6) } TF051.Scale.reScale.anovaVCA <- function() { data(VCAdata1) scaled.fits <- Scale("anovaVCA", y~(device+lot)/day/run, VCAdata1, by="sample") reference.fits <- anovaVCA(y~(device+lot)/day/run, VCAdata1, by="sample") scaled.infs <- VCAinference(scaled.fits, VarVC=TRUE) reference.infs <- VCAinference(reference.fits, VarVC=TRUE) rescaled.infs <- reScale(scaled.infs) for(i in 1:length(scaled.infs)) { scl <- rescaled.infs[[i]] ref <- reference.infs[[i]] sfac <- 1e12/scl$VCAobj$scale # tolerance depending on scaling-factor, i.e. much scaling reduces required accuracy due to less precise reference results that can be generated tol <- 1/sfac checkEquals(ref$VCAobj$aov.tab[,"VC"], scl$VCAobj$aov.tab[,"VC"], tolerance=tol) checkEquals(diag(ref$vcovVC), diag(scl$vcovVC), tolerance=tol) } } #TF052.Scale.reScale.anovaVCA <- function() #{ # data(VCAdata1) # # scaled.fits <- Scale("anovaVCA", y~(device+lot)/day/run, VCAdata1, by="sample") # reference.fits <- anovaVCA(y~(device+lot)/day/run, VCAdata1, by="sample") # scaled.infs <- VCAinference(scaled.fits, VarVC=TRUE) # reference.infs <- VCAinference(reference.fits, VarVC=TRUE) # rescaled.infs <- reScale(scaled.infs) # # for(i in 1:length(scaled.infs)) # { # scl <- rescaled.infs[[i]] # ref <- reference.infs[[i]] # # tol <- 1e-10 # # checkEquals(ref$VCAobj$aov.tab[,"SD"], scl$VCAobj$aov.tab[,"SD"], tolerance=tol) # checkEquals(diag(ref$vcovVC), diag(scl$vcovVC), tolerance=tol) # } #} TF052.Scale.reScale.anovaMM <- function() { data(VCAdata1) scaled.fits <- Scale("anovaMM", y~((device)+(lot))/(day)/(run), VCAdata1, by="sample") reference.fits <- anovaMM( y~((device)+(lot))/(day)/(run), VCAdata1, by="sample") scaled.infs <- VCAinference(scaled.fits, VarVC=TRUE) reference.infs <- VCAinference(reference.fits, VarVC=TRUE) rescaled.infs <- reScale(scaled.infs) for(i in 1:length(scaled.infs)) { scl <- rescaled.infs[[i]] ref <- reference.infs[[i]] tol <- 1e-10 checkEquals(ref$VCAobj$aov.tab[,"SD"], scl$VCAobj$aov.tab[,"SD"], tolerance=tol) checkEquals(diag(ref$vcovVC), diag(scl$vcovVC), tolerance=tol) } } TF053.Scale.reScale.remlVCA <- function() { data(VCAdata1) scaled.fits <- Scale("remlVCA", y~(device+lot)/day/run, VCAdata1, by="sample") reference.fits <- remlVCA(y~(device+lot)/day/run, VCAdata1, by="sample") scaled.infs <- VCAinference(scaled.fits, VarVC=TRUE) reference.infs <- VCAinference(reference.fits, VarVC=TRUE) rescaled.infs <- reScale(scaled.infs) for(i in 1:length(scaled.infs)) { scl <- rescaled.infs[[i]] ref <- reference.infs[[i]] tol <- 1e-5 checkEquals(ref$VCAobj$aov.tab[,"SD"], scl$VCAobj$aov.tab[,"SD"], tolerance=tol) checkEquals(diag(ref$vcovVC), diag(scl$vcovVC), tolerance=tol) } } TF054.Scale.reScale.remlMM <- function() { data(VCAdata1) scaled.fits <- Scale("remlMM", y~((device)+(lot))/(day)/(run), VCAdata1, by="sample") reference.fits <- remlMM( y~((device)+(lot))/(day)/(run), VCAdata1, by="sample") scaled.infs <- VCAinference(scaled.fits, VarVC=TRUE) reference.infs <- VCAinference(reference.fits, VarVC=TRUE) rescaled.infs <- reScale(scaled.infs) for(i in 1:length(scaled.infs)) { scl <- rescaled.infs[[i]] ref <- reference.infs[[i]] tol <- 1e-5 checkEquals(ref$VCAobj$aov.tab[,"SD"], scl$VCAobj$aov.tab[,"SD"], tolerance=tol) checkEquals(diag(ref$vcovVC), diag(scl$vcovVC), tolerance=tol) } } # ensure that fitVCA gives same results as the functions it calls TF055.fitVCA.anovaVCA.remlVCA <- function() { data(VCAdata1) sgnf <- 5 # number of significant digits, for(i in 1:10) # over samples { tmpData <- subset(VCAdata1, sample==i) fit0.anova <- anovaVCA(y~(device+lot)/day/run, tmpData) fit0.reml <- remlVCA( y~(device+lot)/day/run, tmpData) fit1.anova <- fitVCA( y~(device+lot)/day/run, tmpData, method="anova", scale=(i%%2==0)) # alternate whether to scale data or not fit1.reml <- fitVCA( y~(device+lot)/day/run, tmpData, method="reml", scale=(i%%2!=0)) cat("\nsample",i,":\n") print(checkEquals(round(fit0.anova$aov.tab[,"VC"], sgnf), round(fit1.anova$aov.tab[,"VC"], sgnf))) print(checkEquals(round(fit0.reml$aov.tab[ ,"VC"], sgnf), round(fit1.reml$aov.tab[, "VC"], sgnf))) } } TF056.orderData.remlVCA <- function() { data(MLrepro) MLrepro.ord <- with(MLrepro, MLrepro[order(Lab, Lot, Day, Run),]) # manually order factors opt.new <- options(scipen=12) M1.ref <- remlVCA(Result ~ (Lab + Lot)/Day/Run, MLrepro.ord, VarVC = TRUE) M1 <- remlVCA(Result ~ (Lab + Lot)/Day/Run, MLrepro, VarVC = TRUE) checkEquals(M1$aov.tab, M1.ref$aov.tab) options(opt.new) } TF057.orderData.remlMM <- function() { data(MLrepro) MLrepro.ord <- with(MLrepro, MLrepro[order(Lab, Lot, Day, Run),]) # manually order factors opt.new <- options(scipen=12) M1.ref <- remlMM(Result ~ ((Lab) + (Lot))/(Day)/(Run), MLrepro.ord, VarVC = TRUE, order.data = FALSE) M1 <- remlMM(Result ~ ((Lab) + (Lot))/(Day)/(Run), MLrepro, VarVC = TRUE, order.data = TRUE) checkEquals(M1$aov.tab, M1.ref$aov.tab) options(opt.new) } TF058.predict.anovaMM <- function() { # fit LMM with fixed lot and device effects and test for lot-differences data(VCAdata1) VCAdata1_sample5 <- subset(VCAdata1, sample==5) fitS5 <- fitLMM(y~(device+lot)/(day)/(run), VCAdata1_sample5, "anova") # fitted values including fixed and random effects pred <- predict(fitS5) # SAS 9.4 reference ## proc mixed data=VCAdata1_sample5; ## class device lot day run; ## model y=device lot/ s outp=outp; ## random device*lot*day device*lot*day*run; ## run; SASref <- c(17.899050479, 17.899050479, 18.019161777, 18.019161777, 17.851799074, 17.851799074, 17.844764076, 17.844764076, 16.814900319, 16.814900319, 18.099103148, 18.099103148, 18.323024006, 18.323024006, 18.402690956, 18.402690956, 18.02856588, 18.02856588, 17.414336609, 17.414336609, 18.01615201, 18.01615201, 18.150495477, 18.150495477, 17.810710437, 17.810710437, 17.541187839, 17.541187839, 17.891970672, 17.891970672, 17.477222676, 17.477222676, 16.936743521, 16.936743521, 17.009863088, 17.009863088, 16.996838687, 16.996838687, 17.742573591, 17.742573591, 17.181168314, 17.181168314, 17.913059433, 17.913059433, 17.886696332, 17.886696332, 17.920933794, 17.920933794, 17.453001574, 17.453001574, 17.908328271, 17.908328271, 17.168513101, 17.168513101, 17.208107861, 17.208107861, 18.350257975, 18.350257975, 18.454610947, 18.454610947, 18.05353796, 18.05353796, 18.148658917, 18.148658917, 17.782646776, 17.782646776, 18.191321439, 18.191321439, 17.747847363, 17.747847363, 18.160575146, 18.160575146, 17.854959068, 17.854959068, 17.94836637, 17.94836637, 17.778571123, 17.778571123, 17.981690354, 17.981690354, 17.721834759, 17.721834759, 18.210214242, 18.210214242, 17.426262307, 17.426262307, 17.516597821, 17.516597821, 17.995705373, 17.995705373, 17.839060278, 17.839060278, 17.870866129, 17.870866129, 17.455708536, 17.455708536, 17.484376436, 17.484376436, 17.633759632, 17.633759632, 18.059467054, 18.059467054, 18.109651477, 18.109651477, 17.763661497, 17.763661497, 17.678617607, 17.678617607, 17.649497267, 17.649497267, 17.671977329, 17.671977329, 16.844471453, 16.844471453, 16.304322063, 16.304322063, 17.155502328, 17.155502328, 16.959623965, 16.959623965, 17.248848449, 17.248848449, 17.203171369, 17.203171369, 16.377849768, 16.377849768, 16.840993943, 16.840993943, 17.598257565, 17.598257565, 16.993849474, 16.993849474, 15.847775522, 15.847775522, 16.139327587, 16.139327587, 17.235607942, 17.235607942, 17.194286781, 17.194286781, 18.151865586, 18.151865586, 18.006728937, 18.006728937, 17.932605325, 17.932605325, 17.993503934, 17.993503934, 17.972653873, 17.972653873, 17.938784406, 17.938784406, 17.656457414, 17.656457414, 18.171206871, 18.171206871, 18.056473288, 18.056473288, 18.12599667, 18.12599667, 18.195670704, 18.195670704, 18.531011728, 18.531011728, 18.067735776, 18.067735776, 18.445249298, 18.445249298, 17.374363044, 17.374363044, 17.467161232, 17.467161232, 17.095178243, 17.095178243, 17.74331252, 17.74331252, 17.470110835, 17.470110835, 17.734284011, 17.734284011, 17.575357936, 17.575357936, 18.042970992, 18.042970992, 17.151979443, 17.151979443, 17.175104385, 17.175104385, 17.442102655, 17.442102655, 16.823553712, 16.823553712, 17.519946213, 17.519946213, 17.338524873, 17.338524873, 17.147467912, 17.147467912, 16.713978566, 16.713978566, 16.472854882, 16.472854882, 16.65792692, 16.65792692, 15.675853935, 15.675853935, 16.380109642, 16.380109642, 16.915108884, 16.915108884, 17.599223625, 17.599223625, 16.92119403, 16.92119403, 16.700892957, 16.700892957, 15.747910906, 15.747910906, 15.559174464, 15.559174464, 17.193351366, 17.193351366, 17.602150973, 17.602150973, 17.502256685, 17.502256685, 17.701333304, 17.701333304, 17.407119593, 17.407119593, 17.591792282, 17.591792282, 17.188676031, 17.188676031, 17.659079358, 17.659079358, 17.552768564, 17.552768564, 17.347926526, 17.347926526, 16.295013278, 16.295013278, 16.52279957, 16.52279957, 16.232855875, 16.232855875, 15.900764272, 15.900764272, 16.411652347, 16.411652347, 17.135995502, 17.135995502) checkEqualsNumeric(as.vector(pred), SASref) } TF059.predict.remlMM.newdata <- function() { # fit LMM with fixed lot and device effects and test for lot-differences data(VCAdata1) datS5 <- subset(VCAdata1, sample==5) # add continous covariate for changes in newdata set.seed(1) datS5$cov <- round(rnorm(nrow(datS5), 65, 8),1) fitS5 <- fitLMM(y~cov+device+lot+(day)+(run), datS5, "reml") # select 2 observations and edit covariate to calculate predictions for newdata <- datS5[c(1,62),] newdata$cov <- newdata$cov + 3.5 # fitted values including fixed and random effects pred <- predict(fitS5, newdata) # SAS 9.4 reference ## proc mixed data=dats5 order=data; ## class device lot day run; ## model y=cov device lot; ## random day run; ## ods output Estimates=pred_newdata; ## estimate 'newdata_obs1' ## intercept 1 ## cov 63.5 ## device 1 0 0 ## lot 1 0 0 | ## day 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## run 1 0; ## estimate 'newdata_obs2' ## intercept 1 ## cov 68.2 ## device 0 0 1 ## lot 1 0 0 | ## day 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 ## run 1 0; ## run; SASref <- c(18.12473221,17.28640975) checkEqualsNumeric(as.vector(pred), SASref, tolerance=1e-6) } TF060.predict.remlMM.newdata <- function() { data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age-11 fit.fitLMM <- fitLMM(distance~Sex*age2+(Subject)*age2, Ortho, "reml") # user-specified 'newdata' newdata <- Ortho[c(45,75),] newdata$age2 <- newdata$age2 + 5 # fitted values including fixed and random effects pred <- predict(fit.fitLMM, newdata) # SAS 9.4 reference ## proc mixed data=Ortho order=data; ## class sex subj; ## model distance=sex age2 sex*age2; ## random intercept age2 / subject=subj type=un; ## ods output Estimates=pred_newdata; ## estimate 'newdata_obs1' ## intercept 1 ## sex 0 1 ## age2 2 ## sex*age2 0 2 | ## intercept 1 ## age2 2 / subject 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0; ## estimate 'newdata_obs2' ## intercept 1 ## sex 1 0 ## age2 6 ## sex*age2 6 0 | ## intercept 1 ## age2 6 / subject 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ## run; SASref <- c(26.0117,27.2067) checkEqualsNumeric(as.vector(pred), SASref, tolerance=1e-4) } TF061.predict.remlMM.restriction <- function() { data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age-11 fit.fitLMM <- fitLMM(distance~Sex*age2+(Subject)*age2, Ortho, "reml") # restrict to fixed effects only pred <- predict(fit.fitLMM, re=NA) # SAS 9.4 reference ## proc mixed data=Ortho order=data; ## class sex subj; ## model distance=sex age2 sex*age2 / outpm=outpm; ## random intercept age2 / subject=subj type=un; ## run; SASref <- c(21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 21.209090909, 22.168181818, 23.127272727, 24.086363636, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875, 22.615625, 24.184375, 25.753125, 27.321875) checkEqualsNumeric(as.vector(pred), SASref, tolerance=1e-6) } #TODO: More test functions for fitLMM necessary? TF062.comparison.fitLMM.REML <- function() { data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age-11 fit.fitLMM <- fitLMM(distance~Sex*age2+(Subject)*age2, Ortho, "reml") fit.remlMM <- remlMM(distance~Sex*age2+(Subject)*age2, Ortho) checkEqualsNumeric(unlist(fit.fitLMM$aov.tab), unlist(fit.remlMM$aov.tab), tolerance=1e-7) } TF063.comparison.fitLMM.ANOVA <- function() { data(Orthodont) Ortho <- Orthodont Ortho$age2 <- Ortho$age-11 fit.fitLMM <- fitLMM(distance~Sex*age2+(Subject)*age2, Ortho) fit.anovaMM <- anovaMM(distance~Sex*age2+(Subject)*age2, Ortho) checkEqualsNumeric(unlist(fit.fitLMM$aov.tab), unlist(fit.anovaMM$aov.tab), tolerance=1e-8) } TF064.Scale.reScale.response <- function() { data(VCAdata1) sample1 <- VCAdata1[VCAdata1$sample==1,] scaled.fit <- Scale("remlMM", y~((device)+(lot))/(day)/(run), sample1) scaled.inf <- VCAinference(scaled.fit, VarVC=TRUE) rescaled.inf <- reScale(scaled.inf) # order input data accroding to model fit and obtain response variable ref.response <- orderData(sample1 , y~((device)+(lot))/(day)/(run))$y checkEqualsNumeric(rescaled.inf$VCAobj$data$y, ref.response, tolerance=1e-10) } TF065.comparison.fitVCA.REML <- function() { data(dataEP05A2_2) fit.fitVCA <- fitVCA(y~day/run, dataEP05A2_2, "reml") fit.remlVCA <- remlVCA(y~day/run, dataEP05A2_2) checkEqualsNumeric(unlist(fit.fitVCA$aov.tab), unlist(fit.remlVCA$aov.tab), tolerance=1e-8) } TF066.comparison.fitVCA.ANOVA <- function() { data(dataEP05A2_2) fit.fitVCA <- fitVCA(y~day/run, dataEP05A2_2, "anova") fit.anovaVCA <- anovaVCA(y~day/run, dataEP05A2_2) checkEqualsNumeric(fit.fitVCA$aov.tab, fit.anovaVCA$aov.tab, tolerance=1e-8) } TF067.checkException.NoReScaling <- function() { data(VCAdata1) datS511 <- subset(VCAdata1, sample==5 & lot == 1 & device == 1) fit <- Scale("anovaVCA", y~day/run, datS511) options(warn=2) checkException(ranef(fit)) checkException(resid(fit)) options(warn=0) } TF068.check.Scale.reScale.anovaVCA <- function() { data(VCAdata1) datS5 <- subset(VCAdata1, sample == 5) datHV <- datS5 # HV = huge value datHV$y <- datHV$y * 1e7 # make response variable really big numbers #checkException(anovaVCA(y~(device+lot)/day/run, datHV)) # should generate an error fit.anovaVCA <- Scale("anovaVCA", y~(device+lot)/day/run, datHV) # should not generate an error fit.anovaVCA <- reScale(fit.anovaVCA) # check equivalence of scaling against not scaling # the null model (reference) without scaling fit0 <- anovaVCA(y~((device)+(lot))/day/run, datS5) VarCov0 <- vcovVC(fit0) # the model to compare to the null model fit1 <- Scale("anovaVCA", y~((device)+(lot))/day/run, datS5) fit1 <- reScale(fit1, VarVC=TRUE) VarCov1 <- vcovVC(fit1) tol <- 1e-6 digits <- log10(1/tol) checkEquals(round(VarCov0, digits=digits), round(VarCov1, digits=digits), tolerance=tol) # check equality of raw random effects checkEquals(round(ranef(fit0),digits=digits), round(ranef(fit1),digits=digits), tolerance=tol) # check equality of raw studentized random effects checkEquals(round(ranef(fit0, mode="student"),digits=digits), round(ranef(fit1, mode="student"),digits=digits), tolerance=tol) # check equality of raw standardized random effects checkEquals(round(ranef(fit0, mode="standard"),digits=digits), round(ranef(fit1, mode="standard"),digits=digits), tolerance=tol) # check equality of raw conditional residuals checkEquals(round(resid(fit0, term="cond"),digits=digits), round(resid(fit1, term="cond"),digits=digits), tolerance=tol) # check equality of raw studentized conditional residuals checkEquals(round(resid(fit0, term="cond", mode="student"),digits=digits), round(resid(fit1, term="cond", mode="student"),digits=digits), tolerance=tol) # check equality of raw pearson-type transformed conditional residuals checkEquals(round(resid(fit0, term="cond", mode="pearson"),digits=digits), round(resid(fit1,term="cond", mode="pearson"),digits=digits), tolerance=tol) # check equality of raw marginal residuals checkEquals(round(resid(fit0, term="marginal"),digits=digits), round(resid(fit1, term="marginal"),digits=digits), tolerance=tol) # check equality of studentized marginal residuals checkEquals(round(resid(fit0, term="marginal", mode="student"),digits=digits), round(resid(fit1, term="marginal", mode="student"),digits=digits), tolerance=tol) # check equality of pearson-type transformed marginal residuals checkEquals(round(resid(fit0, term="marginal", mode="pearson"),digits=digits), round(resid(fit1,term="marginal", mode="pearson"),digits=digits), tolerance=tol) } TF069.check.Scale.reScale.remlVCA <- function() { data(VCAdata1) datS5 <- subset(VCAdata1, sample == 5) datHV <- datS5 # HV = huge value datHV$y <- datHV$y * 1e7 # make response variable really big numbers #checkException(remlVCA(y~(device+lot)/day/run, datHV)) # should generate an error fit.remlVCA <- Scale("remlVCA", y~(device+lot)/day/run, datHV) # should not generate an error fit.remlVCA <- reScale(fit.remlVCA) # check equivalence of scaling against not scaling # the null model (reference) without scaling fit0 <- remlVCA(y~((device)+(lot))/day/run, datS5) VarCov0 <- vcovVC(fit0) # the model to compare to the null model fit1 <- Scale("remlVCA", y~((device)+(lot))/day/run, datS5) fit1 <- reScale(fit1, VarVC=TRUE) VarCov1 <- vcovVC(fit1) tol <- 1e-4 # lower precision due to differences in the results of the numerical optimization algorithm applied to data on different scales digits <- log10(1/tol) checkEquals(round(VarCov0, 6), round(VarCov1, 6), tolerance=1e-6) # check equality of raw random effects checkEquals(round(ranef(fit0),digits=digits), round(ranef(fit1),digits=digits), tolerance=tol) # check equality of raw studentized random effects checkEquals(round(ranef(fit0, mode="student"),digits=digits), round(ranef(fit1, mode="student"),digits=digits), tolerance=tol) # check equality of raw standardized random effects checkEquals(round(ranef(fit0, mode="standard"),digits=digits), round(ranef(fit1, mode="standard"),digits=digits), tolerance=tol) # check equality of raw conditional residuals checkEquals(round(resid(fit0, term="cond"),digits=digits), round(resid(fit1, term="cond"),digits=digits), tolerance=tol) # check equality of raw studentized conditional residuals checkEquals(round(resid(fit0, term="cond", mode="student"),digits=digits), round(resid(fit1, term="cond", mode="student"),digits=digits), tolerance=tol) # check equality of raw pearson-type transformed conditional residuals checkEquals(round(resid(fit0, term="cond", mode="pearson"),digits=digits), round(resid(fit1,term="cond", mode="pearson"),digits=digits), tolerance=tol) # check equality of raw marginal residuals checkEquals(round(resid(fit0, term="marginal"),digits=digits), round(resid(fit1, term="marginal"),digits=digits), tolerance=tol) # check equality of studentized marginal residuals checkEquals(round(resid(fit0, term="marginal", mode="student"),digits=digits), round(resid(fit1, term="marginal", mode="student"),digits=digits), tolerance=tol) # check equality of pearson-type transformed marginal residuals checkEquals(round(resid(fit0, term="marginal", mode="pearson"),digits=digits), round(resid(fit1,term="marginal", mode="pearson"),digits=digits), tolerance=tol) } TF070.check.Scale.reScale.anovaMM <- function() { data(VCAdata1) datS5 <- subset(VCAdata1, sample == 5) datHV <- datS5 # HV = huge value datHV$y <- datHV$y * 1e7 # make response variable really big numbers #checkException(anovaMM(y~(device+lot)/(day)/(run), datHV)) # should generate an error fit.anovaMM <- Scale("anovaMM", y~(device+lot)/(day)/(run), datHV) # should not generate an error fit.anovaMM <- reScale(fit.anovaMM) # check equivalence of scaling against not scaling # the null model (reference) without scaling fit0 <- anovaMM(y~(device+lot)/(day)/(run), datS5) VarCov0 <- vcovVC(fit0) # the model to compare to the null model fit1 <- Scale("anovaMM", y~(device+lot)/(day)/(run), datS5) fit1 <- reScale(fit1, VarVC=TRUE) VarCov1 <- vcovVC(fit1) tol <- 1e-6 # lower precision due to differences in the results of the numerical optimization algorithm applied to data on different scales digits <- log10(1/tol) checkEquals(round(VarCov0, 6), round(VarCov1, 6), tolerance=1e-6) # check equality of raw random effects checkEquals(round(ranef(fit0),digits=digits), round(ranef(fit1),digits=digits), tolerance=tol) # check equality of raw studentized random effects checkEquals(round(ranef(fit0, mode="student"),digits=digits), round(ranef(fit1, mode="student"),digits=digits), tolerance=tol) # check equality of raw standardized random effects checkEquals(round(ranef(fit0, mode="standard"),digits=digits), round(ranef(fit1, mode="standard"),digits=digits), tolerance=tol) # check equality of raw conditional residuals checkEquals(round(resid(fit0, term="cond"),digits=digits), round(resid(fit1, term="cond"),digits=digits), tolerance=tol) # check equality of raw studentized conditional residuals checkEquals(round(resid(fit0, term="cond", mode="student"),digits=digits), round(resid(fit1, term="cond", mode="student"),digits=digits), tolerance=tol) # check equality of raw pearson-type transformed conditional residuals checkEquals(round(resid(fit0, term="cond", mode="pearson"),digits=digits), round(resid(fit1,term="cond", mode="pearson"),digits=digits), tolerance=tol) # check equality of raw marginal residuals checkEquals(round(resid(fit0, term="marginal"),digits=digits), round(resid(fit1, term="marginal"),digits=digits), tolerance=tol) # check equality of studentized marginal residuals checkEquals(round(resid(fit0, term="marginal", mode="student"),digits=digits), round(resid(fit1, term="marginal", mode="student"),digits=digits), tolerance=tol) # check equality of pearson-type transformed marginal residuals checkEquals(round(resid(fit0, term="marginal", mode="pearson"),digits=digits), round(resid(fit1,term="marginal", mode="pearson"),digits=digits), tolerance=tol) } TF071.check.Scale.reScale.remlMM <- function() { data(VCAdata1) datS5 <- subset(VCAdata1, sample == 5) datHV <- datS5 # HV = huge value datHV$y <- datHV$y * 1e7 # make response variable really big numbers #checkException(remlMM(y~(device+lot)/(day)/(run), datHV)) # should generate an error fit.remlMM <- Scale("remlMM", y~(device+lot)/(day)/(run), datHV) # should not generate an error fit.remlMM <- reScale(fit.remlMM) # check equivalence of scaling against not scaling # the null model (reference) without scaling fit0 <- remlMM(y~(device+lot)/(day)/(run), datS5) VarCov0 <- vcovVC(fit0) # the model to compare to the null model fit1 <- Scale("remlMM", y~(device+lot)/(day)/(run), datS5) fit1 <- reScale(fit1, VarVC=TRUE) VarCov1 <- vcovVC(fit1) tol <- 1e-5 # lower precision due to differences in the results of the numerical optimization algorithm applied to data on different scales digits <- log10(1/tol) checkEquals(round(VarCov0, 6), round(VarCov1, 6), tolerance=1e-6) # check equality of raw random effects checkEquals(round(ranef(fit0),digits=digits), round(ranef(fit1),digits=digits), tolerance=tol) # check equality of raw studentized random effects checkEquals(round(ranef(fit0, mode="student"),digits=digits), round(ranef(fit1, mode="student"),digits=digits), tolerance=tol) # check equality of raw standardized random effects checkEquals(round(ranef(fit0, mode="standard"),digits=digits), round(ranef(fit1, mode="standard"),digits=digits), tolerance=tol) # check equality of raw conditional residuals checkEquals(round(resid(fit0, term="cond"),digits=digits), round(resid(fit1, term="cond"),digits=digits), tolerance=tol) # check equality of raw studentized conditional residuals checkEquals(round(resid(fit0, term="cond", mode="student"),digits=digits), round(resid(fit1, term="cond", mode="student"),digits=digits), tolerance=tol) # check equality of raw pearson-type transformed conditional residuals checkEquals(round(resid(fit0, term="cond", mode="pearson"),digits=digits), round(resid(fit1,term="cond", mode="pearson"),digits=digits), tolerance=tol) # check equality of raw marginal residuals checkEquals(round(resid(fit0, term="marginal"),digits=digits), round(resid(fit1, term="marginal"),digits=digits), tolerance=tol) # check equality of studentized marginal residuals checkEquals(round(resid(fit0, term="marginal", mode="student"),digits=digits), round(resid(fit1, term="marginal", mode="student"),digits=digits), tolerance=tol) # check equality of pearson-type transformed marginal residuals checkEquals(round(resid(fit0, term="marginal", mode="pearson"),digits=digits), round(resid(fit1,term="marginal", mode="pearson"),digits=digits), tolerance=tol) } # check correct implementation of model with multiple numeric covariables # numeric examples taken from # Jonas Haslbeck # University of Amsterdam # Psychological Methods # Nieuwe Achtergracht 129-B # Postbus 15906 # 1018 WT, Amsterdam # the Netherlands # jonashaslbeck@gmail.com # http://jmbh.github.io TF072.getMM.multi.covariables <- function() { # example provided by Jonas Haslbeck data <- as.data.frame(matrix(rnorm(40), nrow=10)) colnames(data) <- c('V.1', 'V.2', 'V.3', 'V.4') form <- as.formula('V.1 ~ V.2 + V.3 + V.4 + V.2 * V.3 + V.2 * V.4 + V.3 * V.4 + V.2 * V.3 * V.4') X0 <- model.matrix(form, data = data) X1 <- getMM(form, Data = data) checkEquals(c(X0), c(as.matrix(X1))) # check equality to mix data2 <- data data2$subject <- gl(5, 2) form2 <- as.formula('V.1 ~ V.2 + subject*V.2 * V.3') X3 <- model.matrix(form2, data2) X4 <- getMM(form2, data2) # exclude extra-columns in X4 for comparison and re-arrage column-order (content counts) checkEquals(c(X3[,c(1:6, 8:11, 7, 12:20)]), c(as.matrix(X4)[,-c(3,8,14,20)])) } # Check whether CI for intermediary VC are correctly estimated when fitting model # using 'fitVCA' instead of 'anovaVCA' # In version 1.4.1 vcov for intermediary VC were not estimated. TF073.fitVCA_ANOVA_VarVC_TRUE <- function() { data(dataEP05A2_3) fit0 <- anovaVCA(y~day/run, dataEP05A2_3) fit1 <- fitVCA( y~day/run, dataEP05A2_3) inf0 <- VCAinference(fit0, VarVC=TRUE) inf1 <- VCAinference(fit1, VarVC=TRUE) res0 <- unlist(lapply(inf0$ConfInt, function(x) return(c(x$OneSided$LCL, x$OneSided$UCL, x$TwoSided$LCL, x$TwoSided$UCL)))) res1 <- unlist(lapply(inf1$ConfInt, function(x) return(c(x$OneSided$LCL, x$OneSided$UCL, x$TwoSided$LCL, x$TwoSided$UCL)))) checkEquals(res0, res1) } TF074.fitVCA_REML_VarVC_TRUE <- function() { fit0 <- remlVCA(y~day/run, dataEP05A2_3) fit1 <- fitVCA( y~day/run, dataEP05A2_3, "REML") inf0 <- VCAinference(fit0, VarVC=TRUE) inf1 <- VCAinference(fit1, VarVC=TRUE) res0 <- unlist(lapply(inf0$ConfInt, function(x) return(c(x$OneSided$LCL, x$OneSided$UCL, x$TwoSided$LCL, x$TwoSided$UCL)))) res1 <- unlist(lapply(inf1$ConfInt, function(x) return(c(x$OneSided$LCL, x$OneSided$UCL, x$TwoSided$LCL, x$TwoSided$UCL)))) checkEquals(res0, res1, tol=1e-6) } ##################################################################### #### check protectedCAll functionality on all model-fitting functions TF075.protectedCall.anovaVCA <- function() { dat3 <- data.frame( y=rnorm(8,10), day=paste("day",rep(c(1,2),c(4,4)), sep="_"), run=rep(c(2,1), c(4,4))) assign("dat3", dat3, envir=.GlobalEnv) smpl <- try(protectedCall(anovaVCA(y~day/run, dat3), ErrorType="Simple"), silent=TRUE) dtld <- try(protectedCall(anovaVCA(y~day/run, dat3), ErrorType="Detailed"), silent=TRUE) smpl <- attr(smpl, "condition")$message dtld <- attr(dtld, "condition")$message cat("\nTF075.protectedCall.anovaVCA\nsmpl:",smpl) cat("\ndtld:", dtld) checkTrue(nchar(smpl)