# 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/UnitTests/LSMeans_Data.RData"))

fit.vca <- remlMM(y~snp+time+snp:time+sex+(id)+(id):time, dat,VarVC=F)

# 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(fit.fitLMM$aov.tab, 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(fit.fitLMM$aov.tab, 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(fit.fitVCA$aov.tab, 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)
	checkEquals(unlist(inf0$ConfInt), unlist(inf1$ConfInt))
}


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)
	checkEquals(unlist(inf0$ConfInt), unlist(inf1$ConfInt), tol=1e-6)
}


