R/vca.R

Defines functions stepwiseVCA residuals.VCA anovaVCA as.matrix.VCAinference as.matrix.VCA Trace isBalanced print.VCAinference VCAinference print.VCA SattDF getMM getVCvar MPinv getAmatrix getCmatrix getV getMat vcovVC getDDFM getL test.fixef test.lsmeans DfSattHelper vcovFixed getDF vcov.VCA coef.VCA lsmMat lsmeans fixef.VCA fixef ranef.VCA ranef solveMME anovaDF anovaMM getSSQqf getBasis getSSQsweep Sinv Csweep Solve remlMM remlVCA lmerSummary lmerMatrices lmerG getGB check4MKL reScale Scale scaleData load_if_installed fitVCA fitLMM model.matrix.VCA predict.VCA model.frame.VCA orderData onUnload onLoad

Documented in anovaDF anovaMM anovaVCA as.matrix.VCA as.matrix.VCAinference check4MKL coef.VCA Csweep DfSattHelper fitLMM fitVCA fixef fixef.VCA getAmatrix getBasis getCmatrix getDDFM getDF getGB getL getMat getMM getSSQqf getSSQsweep getV getVCvar isBalanced lmerG lmerMatrices lmerSummary load_if_installed lsmeans lsmMat model.frame.VCA model.matrix.VCA MPinv orderData predict.VCA print.VCA print.VCAinference ranef ranef.VCA remlMM remlVCA reScale residuals.VCA SattDF Scale scaleData Sinv Solve solveMME stepwiseVCA test.fixef test.lsmeans Trace VCAinference vcovFixed vcovVC vcov.VCA

# Load/unload C-lib.

.onLoad <- function(lib, pkg)
{
	library.dynam(chname="VCA", package=pkg, lib.loc=lib)
	# create VCA message environment
	msgEnv <<- new.env(parent=emptyenv())
	check4MKL()
}

.onUnload <- function(lib)
{
	library.dynam.unload(chname="VCA", libpath=lib)
}	


#' Re-Order Data.Frame.
#' 
#' Functions attempts to standardize input data for linear mixed model analyses
#' to overcome the problem that analysis results sometimes depend on ordering of
#' the data and definition of factor-levels.
#' 
#' @param Data				(data.frame) with input data intented to put into standard-order
#' @param trms				(formula, terms) object speciying a model to be fitted to \code{Data}
#' @param order.data		(logical) TRUE = variables will be increasingly ordered, FALSE = order of
#' 							the variables remains as is
#' @param exclude.numeric	(logical) TRUE = numeric variables will not be included in the reordering,
#' 							which is required whenever this variable serves as covariate in a LMM, 
#' 							FALSE = numeric variables will also be converted to factors, useful in 
#' 							VCA-analysis, where all variables are interpreted as class-variables	
#' @param quiet				(logical) TRUE = omits any (potentially) informative output regarding
#' 							re-ordering and type-casting of variables
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' @examples 
#' \dontrun{
#' # random ordering
#' data(dataEP05A2_1)
#' dat <- dataEP05A2_1
#' levels(dat$day) <- sample(levels(dat$day))
#' # this has direct impact e.g. on order of estimated effects
#' fit <- anovaVCA(y~day/run, dat, order.data=FALSE)
#' ranef(fit)
#' # to guarantee consistent analysis results
#' # independent of the any data orderings option
#' # 'order.data' is per default set to TRUE:
#' fit <- anovaVCA(y~day/run, dat)
#' ranef(fit)
#' # which is identical to:
#' fit2 <- anovaVCA(y~day/run, orderData(dat, y~day/run), order.data=FALSE)
#' ranef(fit2)
#' }

orderData <- function(Data, trms, order.data=TRUE, exclude.numeric=TRUE, quiet=FALSE)
{
	stopifnot(is.data.frame(Data))
	stopifnot("formula" %in% class(trms))
	if(!"terms" %in% class(trms))
		trms <- terms(trms)
	vars <- rownames(attr(trms, "factors"))[-1]
	if(length(vars) == 0)
		return(Data)
	stopifnot(all(vars %in% colnames(Data)))
	ord.vars <- NULL
	
	for(i in rev(vars))
	{
		if(is.numeric(Data[,i]) && exclude.numeric)
			next
		else
		{
			if(!quiet && is.character(Data[,i]))
				message("Convert variable ", i," from \"character\" to \"factor\"!")
			ord.vars <- c(paste0("Data[,\"",i,"\"]"), ord.vars)
			cha <- as.character(Data[,i])
			num <- suppressWarnings(as.numeric(cha))						# warnings are like to occur here
			if(!any(is.na(num)) && !any(grepl("0[[:digit:]]+", cha)))		# also handle elements like "01", "001", ... those will not give NAs converting them to numeric
			{
				if(order.data)
					Data[,i] <- factor(cha, levels=as.character(sort(unique(num))))
				else
					Data[,i] <- factor(cha, levels=as.character(unique(num)))
			}
			else
			{
				if(order.data)
					Data[,i] <- factor(cha, levels=sort(unique(cha)))
				else
					Data[,i] <- factor(cha, levels=unique(cha))
			}
		}
	}
	if(order.data)
	{
		exprs <- paste0("Data <- Data[order(", paste(ord.vars, collapse=","), "),]")
		eval(parse(text=exprs))
	}
	Data
}


#' Extract the Model Frame from a 'VCA' Object.
#' 
#' Function returns the data-element of 'object' and
#' adds the terms-element as attribute. 
#' 
#' It enables application of functions relying on the existence of 
#' this method, e.g. the functin 'glht' of the 'multcomp'
#' R-package.
#' 
#' @param formula		(VCA) object
#' @param ...			additional arguments
#' @return (data.frame) with attribute 'terms'
#' @method model.frame VCA
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}

model.frame.VCA <- function(formula, ...)
{
	stopifnot(class(formula) == "VCA")
	object <- formula
	rn <- rownames(attr(object$terms, "factors"))
	MF <- object$data[,rn]
	attr(MF, "terms") <- object$terms
	MF
}

#' Predictions from a Model Fitted by \code{fitLMM}.
#' 
#' Model returns fitted values in case \code{newdata} is NULL or evaluates
#' the fitted model employing user-specified data \code{newdata}. The default is that
#' fitted values incorporate fixed effects and random effects, leaving out the (conditional)
#' residuals only. If the interest lies in constraining predictions to the fixed effects only
#' set \code{re=NA} or incorporate just part of the random variability specifying distinct random
#' effects (see \code{re}.
#' 
#' @param object			(VCA) object fitted via function \code{\link{fitLMM}}
#' @param newdata			(data.frame) with all variables required for the specified prediction,
#' 							i.e. the default settings require all variables of the original model, 
#' 							in case of \code{re=NA}, only variables corresponding to fixed effects are
#' 							required.
#' @param re				(character) if NULL (default) all random effects will be included, 
#' 							to restrict predictions to the fixed effects use \code{re=NA}, for
#' 							a subset of random effects included in predictions use any valid
#' 							random effects specification, i.e. \code{object$random}
#' @param allow.new.levels	(logical) if new levels (no part of the original fitted model) in newdata
#' 							are allowed. If FALSE (default), such new values in newdata will trigger
#' 							an error; if TRUE, then the prediction will use the unconditional (population-level)
#' 							values for data with previously unobserved levels (or NAs).
#' @param ...				additional arguments passdo or from other methods
#' 
#' @return (numeric) vector of prediction results
#' @method predict VCA
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples 
#' \dontrun{
#' # fit LMM with fixed lot and device effects and test for lot-differences
#' data(VCAdata1)
#' datS5 <- subset(VCAdata1, sample==5)
#' fitS5 <- fitLMM(y~(lot+device)/(day)/(run), datS5, "anova")
#' fitS5
#' 
#' # fitted values including fixed and random effects
#' pred0 <- predict(fitS5)
#' pred0
#' # sanity check:
#' all(round(pred0 + resid(fitS5) - datS5$y, 12) == 0)
#' # restrict to fixed effects
#' predict(fitS5, re=NA)
#' # restrict to fixed effects and dayly random effects
#' # see required names
#' fitS5$random
#' predict(fitS5, re="lot:device:day")
#' 
#' # check against original 'lmer'-predictions
#' # use version from VCA-package (ordinary data.frame)
#' data(Orthodont, package="VCA")
#' Ortho <- Orthodont
#' Ortho$age2 <- Ortho$age-11
#' # use exactly the same data, same ordering
#' Ortho <- orderData(Ortho, distance ~ Sex * age2 + (Subject) * age2)
#' fit.fitLMM <- fitLMM(distance ~ Sex * age2 + (Subject) * age2, Ortho, "reml")
#' library(lme4)
#' fit.lmer <- lmer(distance ~ Sex + age2 + Sex:age2 + (age2 | Subject), Ortho)
#' # check fitted value first (fixed + random effects)
#' predict(fit.lmer)
#' predict(fit.fitLMM)
#' # restrict to fixed part only
#' predict(fit.lmer, re.form=NA)  
#' predict(fit.fitLMM, re=NA)
#' # user-specified 'newdata'
#' newdata <- Ortho[45:54,]
#' newdata$age2 <- newdata$age2 + 5
#' # include fixed and random effects
#' predict(fit.lmer, newdata)
#' predict(fit.fitLMM, newdata)
#' # generate new data
#' newdata <- Ortho[45:54,]          
#' newdata$age2 <- newdata$age2 + 5
#' # predict on newdata using fixed and random effects
#' predict(fit.lmer, newdata) 
#' predict(fit.fitLMM, newdata)       
#' # restrict prediction to random Subject effects
#' predict(fit.lmer, newdata, re.form=~(1|Subject))        
#' predict(fit.fitLMM, newdata, re="Subject")
#' # restrict prediction to random per-Subject slope
#' predict(fit.lmer, newdata, re.form=~(age2-1|Subject)) 
#' predict(fit.fitLMM, newdata, re="age2:Subject")
#' }

predict.VCA <- function(object, newdata=NULL, re=NULL, allow.new.levels=FALSE, ...)
{
	if(is.null(re))								# all random effects per default
	{
		re <- object$random
		if(is.null(re))							# VCA-model
		{
			if(object$EstMethod == "REML")
				re <- object$VCnames
			else
				re <- object$VCnames[-length(object$VCnames)]
		}
	}
	else
	{
		if(!is.na(re[1]) && !all(re %in% object$random))
			stop("You can only specify random effects in 're' which are available in 'object$random'!")
	}
	
	user.nd <- FALSE							# indicate that no user-defined newdata was specified (default)
	
	if(!is.na(re[1]) && !is.null(newdata))		# all variables found in 're' in 'newdata'?
	{
		re.trms <-  unique(unlist(strsplit(object$random, ":")))
		cnd		<- colnames(newdata)
		if(any(!re.trms %in% cnd))
			stop(paste0("Variable(s): ", paste(re.trms[which(!re.trms %in% cnd)], collapse=", "), " are missing in 'newdata'!"))
	}
	
	if(is.null(attr(object$Type, "fitted.by")))
	{
		warning("Predictions can only be made on models fitted by function 'fitLMM' or 'fitVCA'!")
		return()
	}
	trms  	<- object$fixed.terms
	
	### fixed part first	
	if(is.null(newdata))
	{
		X <- getMat(object, "X")				# get design matrix of fixed effects
		newdata <- object$data	
	}
	else
	{
		user.nd <- TRUE							# indicate user-defined newdata
		Terms	<- delete.response(trms)
		mf 		<- model.frame(Terms, newdata, xlev=object$xlevels)	
		cntr    <- as.list(rep("contr.SAS", length(object$xlevels)))	# use SAS-type contrasts, i.e. last level will be constrained 
		names(cntr) <- names(object$xlevels)
		X		<- model.matrix(Terms, mf, contrasts.arg=cntr)
		asgn	<- attr(X, "assign")
		asgn0	<- object$fe.assign
		tdiff	<- table(asgn0) - table(asgn)	# missing columns in X?
		
		if(any(tdiff == 0))
		{
			X0  <- matrix(0, nrow=nrow(X), ncol=length(asgn0))		# pre-allocate matrix of right dimension
			ind <- match(asgn, asgn0) 								# add all-zero columns to match vector of fixed effects
			for(i in unique(ind))
			{
				tmp <- which(ind == i)
				if(length(tmp) > 1)
					ind[tmp] <- ind[tmp] + 0:(length(tmp)-1)
			}
			X0[,ind] <- X
			X <- X0
		}
	}
	fe 		<- fixef(object, quiet=TRUE)
	preds	<- as.numeric(as.matrix(X) %*% fe[,"Estimate", drop=F])
	
	### add random part to predictions
	if(!is.na(re[1]))
	{
		#Z 	<- getMat(object, "Z")
		ref <- ranef(object, quiet=TRUE)						# at this point all random effects would be used
		
		# generate names for random effects from 'newdata' accounting for user-specified 'newdata' 
		#fac   <- na.omit(names(object$terms.classes[re]))		# only these will get names like VarnameVarlevel 
		rtrms <- re
		#renam <- NULL		
		Z	  <- matrix(nrow=nrow(newdata), ncol=0)
		
		for(i in 1:length(rtrms))						# over all random terms
		{
			tmpZ <- getMM(as.formula(paste0("~", rtrms[i], "-1")), newdata) 	# model.matrix wo intercept
			Z  	 <- cbind(Z, tmpZ) 												# build complete Z-matrix corresponding to newdata
		}	
		
		ref <- ref[colnames(Z),,drop=F]						# align random effects to column-order of Z
		preds <- preds + as.numeric(Z %*% ref)		
	}	
	names(preds) <- rownames(newdata)
	return(preds)
}


#' Model Matrix of a Fitted VCA-Object.
#' 
#' Function returns matrix \code{X} corresponding
#' to the design matrix of fixed effects of the fitted
#' model.
#' 
#' @param object			(VCA) object
#' @param ...				further arguments
#' @method model.matrix VCA

model.matrix.VCA <- function(object, ...)
{
	getMat(object, "X")
}


#' Fit Linear Mixed Model by ANOVA or REML.
#' 
#' Function serves as interface to functions \code{\link{anovaMM}} and \code{\link{remlMM}}
#' for fitting a linear mixed model (LMM) either by ANOVA or REML. All arguments applicable
#' to either one of these functions can be specified (see \code{\link{anovaMM}} or \code{\link{remlMM}} for details).
#' 
#' Besides offering a convenient interface to both functions for fitting a LMM, this function also provides all elements
#' required for standard task of fitted models, e.g. prediction, testing general linear hypotheses via R-package \code{multcomp},
#' etc. (see examples).
#' 
#' @param form		(formula) specifiying the linear mixed model, random effects need to be identified by enclosing
#' 					them in round brackets, i.e. ~a/(b) will model factor 'a' as fixed and 'b' as random
#' @param Data		(data.frame)  containing all variables referenced in 'form', note that variables can only be of type
#'  				"numeric", "factor" or "character". The latter will be automatically converted to "factor"
#' @param method	(character) either "anova" to use ANOVA Type-I estimation of variance components or "reml" to use 
#' 					restricted maximum likelihood (REML) estimation of variance component		
#' @param scale		(logical) TRUE = scale values of the response aiming to avoid numerical problems
#' 					when numbers are either very small or very large, FALSE = use original scale
#' @param VarVC		(logical) TRUE = variance-covariance matrix of variance components will be computed, FALSE = it will not
#' 					be computed
#' @param ...		additional arguments to be passed to function \code{\link{anovaMM}} or function \code{\link{remlMM}}.
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @seealso \code{\link{fitVCA}}, \code{\link{anovaMM}}, \code{\link{remlMM}}
#' 
#' @examples
#' \dontrun{
#' data(dataEP05A2_2)
#' 
#' # assuming 'day' as fixed, 'run' as random
#' # Note: default method is "anova"
#' fitLMM(y~day/(run), dataEP05A2_2)
#' 
#' # explicitly request "reml"
#' fitLMM(y~day/(run), dataEP05A2_2, method="reml")
#' 
#' # assuming both as random leads to same results as
#' # calling anovaVCA (ANOVA is the default)
#' fitLMM(y~(day)/(run), dataEP05A2_2)
#' anovaVCA(y~day/run, dataEP05A2_2)
#' 
#' # now using REML-estimation
#' fitLMM(y~(day)/(run), dataEP05A2_2, "reml")
#' remlVCA(y~day/run, dataEP05A2_2)
#' 
#' # use different approaches to estimating the covariance of 
#' # variance components (covariance parameters)
#' # create unbalanced data
#' dat.ub <- dataEP05A2_2[-c(11,12,23,32,40,41,42),]
#' m1.ub <- fitLMM(y~day/(run), dat.ub, SSQ.method="qf", VarVC.method="scm")
#' # VarVC.method="gb" is an approximation not relying on quadratic forms
#' m2.ub <- fitLMM(y~day/(run), dat.ub, SSQ.method="qf", VarVC.method="gb")		
#' # REML-estimated variance components usually differ from ANOVA-estimates
#' # and so do the variance-covariance matrices
#' m3.ub <- fitLMM(y~day/(run), dat.ub, "reml", VarVC=TRUE)		
#' V1.ub <- round(vcovVC(m1.ub), 12)
#' V2.ub <- round(vcovVC(m2.ub), 12)
#' V3.ub <- round(vcovVC(m3.ub), 12)
#' 
#' # fit a larger random model
#' data(VCAdata1)
#' fitMM1 <- fitLMM(y~((lot)+(device))/(day)/(run), VCAdata1[VCAdata1$sample==1,])
#' fitMM1
#' # now use function tailored for random models
#' fitRM1 <- anovaVCA(y~(lot+device)/day/run, VCAdata1[VCAdata1$sample==1,])
#' fitRM1
#' 
#' # there are only 3 lots, take 'lot' as fixed 
#' fitMM2 <- fitLMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,])
#' # use REML on this (balanced) data
#' fitMM2.2 <- fitLMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,], "reml")
#' 
#' # the following model definition is equivalent to the one above,
#' # since a single random term in an interaction makes the interaction
#' # random (see the 3rd reference for details on this topic)
#' fitMM3 <- fitLMM(y~(lot+(device))/day/run, VCAdata1[VCAdata1$sample==2,])
#' 
#' # fit same model for each sample using by-processing
#' lst <- fitLMM(y~(lot+(device))/day/run, VCAdata1, by="sample")
#' lst
#' 
#' # fit mixed model originally from 'nlme' package
#'  
#' library(nlme)
#' data(Orthodont)
#' fit.lme <- lme(distance~Sex*I(age-11), random=~I(age-11)|Subject, Orthodont) 
#' 
#' # re-organize data
#' Ortho <- Orthodont
#' Ortho$age2 <- Ortho$age - 11
#' Ortho$Subject <- factor(as.character(Ortho$Subject))
#' fit.anovaMM1 <- fitLMM(distance~Sex*age2+(Subject)*age2, Ortho)
#' 
#' # use simplified formula avoiding unnecessary terms
#' fit.anovaMM2 <- fitLMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho)
#' 
#' # and exclude intercept
#' fit.anovaMM3 <- fitLMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho)
#' 
#' # compare results
#' fit.lme
#' fit.anovaMM1
#' fit.anovaMM2
#' fit.anovaMM3
#' 
#' # are there a sex-specific differences?
#' cmat <- getL(fit.anovaMM3, c("SexMale-SexFemale", "SexMale:age2-SexFemale:age2")) 
#' cmat
#' 			 
#' test.fixef(fit.anovaMM3, L=cmat)
#' 
#' # fit LMM with fixed lot and device effects and test for lot-differences
#' data(VCAdata1)
#' fitS5 <- fitLMM(y~(lot+device)/(day)/(run), subset(VCAdata1, sample==5), "reml")
#' fitS5
#' 
#' # apply Tukey-HSD test to screen for lot differences
#' library(multcomp)
#' res.tuk <- glht(fitS5, linfct=mcp(lot="Tukey"))
#' summary(res.tuk)
#' 
#' # compact letter display
#' res.tuk.cld <- cld(res.tuk, col=paste0("gray", c(90,60,75)))
#' plot(res.tuk.cld)
#' }

fitLMM <- function(	form, Data, method=c("anova", "reml"), scale=TRUE, VarVC=TRUE, ...)
{
	call 	<- match.call()
	method	<- match.arg(tolower(method[1]), choices=c("anova", "reml"))
	Args	<- list(...)
	
	Args$form  <- form
	Args$Data  <- Data
	
	if(scale)
	{
		if(method=="anova")
			Args$Fun  <- "anovaMM"			
		else
		{
			Args$Fun  	<- "remlMM"
			Args$VarVC 	<- VarVC
		}
		fit <- do.call("Scale", Args)
	}
	else
	{
		if(method=="anova")
			fit <- do.call("anovaMM", Args)
		else
		{
			Args$VarVC 	<- VarVC
			fit 		<- do.call("remlMM", Args)
		}
	}
	
	wasVCA <- FALSE								# record that result was VCA-object
	
	if(class(fit) == "VCA")						# if by-processing was not used 			
	{
		fit <- list(fit)
		wasVCA <- TRUE
	}
	
	for(i in 1:length(fit))								# over each list-element
	{
		if(scale)
			fit[[i]] <- reScale(fit[[i]], VarVC=VarVC)	# first re-scale
		
		fit[[i]]$call 			<- call											# prevent non-informative object-name, e.g. "form"
		fe 						<- fixef(fit[[i]])
		X  						<- getMat(fit[[i]], "X")		
		fit[[i]]$fitted.values	<- as.numeric(as.matrix(X) %*% fe[,"Estimate", drop=F])
		fixed 					<- fit[[i]]$fixed[!grepl(":", fit[[i]]$fixed)]	# no interactions
		fixed					<- sapply(fit[[i]]$data[,fixed,drop=F], class)
		fixed					<- names(fixed[which(fixed != "numeric")])
		fit[[i]]$xlevels 		<- lapply(fixed, function(x) as.character(unique(fit[[i]]$data[,x])))
		names(fit[[i]]$xlevels) <- fixed
		
		attr(fit[[i]]$Type, "fitted.by") <- "fitLMM"	# only LMM fitted by 'fitLMM' can be handled in standard ways
	}
	
	if(wasVCA)								
		fit <- fit[[1]]
	
	fit
}



#' Fit Variance Component Model by ANOVA or REML.
#' 
#' Function serves as interface to functions \code{\link{anovaVCA}} and \code{\link{remlVCA}}
#' for fitting a variance component models (random models) either by ANOVA or REML. All arguments applicable
#' to either one of these functions can be specified (see \code{\link{anovaVCA}} or \code{\link{remlVCA}} for details).
#' 
#' @param form		(formula) specifiying the variance component model (see \code{\link{anovaVCA}} and/or \code{\link{remlVCA}})
#' @param Data		(data.frame)  containing all variables referenced in 'form'
#' @param method	(character) either "anova" to use ANOVA Type-I estimation of variance components or "reml" to use 
#' 					restricted maximum likelihood (REML) estimation of variance component		
#' @param scale		(logical) TRUE = scale values of the response aiming to avoid numerical problems
#' 					when numbers are either very small or very large, FALSE = use original scale
#' @param VarVC		(logical) TRUE = variance-covariance matrix of variance components will be computed, FALSE = it will not
#' 					be computed
#' @param ...		additional arguments to be passed to function \code{\link{anovaVCA}} or function \code{\link{remlVCA}}.
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @seealso \code{\link{fitLMM}}, \code{\link{anovaVCA}}, \code{\link{remlVCA}}
#' 
#' @examples
#' \dontrun{
#' #load data (CLSI EP05-A2 Within-Lab Precision Experiment) 
#' data(dataEP05A2_2)
#' 
#' # perform ANOVA-estimation of variance components
#' res.anova <- fitVCA(y~day/run, dataEP05A2_2, "anova")
#' # perform REML-estimation of variance components
#' res.reml <- fitVCA(y~day/run, dataEP05A2_2, "reml")
#' 
#' # compare scaling vs. not scaling the response
#' fit0 <- fitVCA(y~day/run, dataEP05A2_2, "anova", scale=TRUE)
#' fit1 <- fitVCA(y~day/run, dataEP05A2_2, "anova", scale=FALSE)
#' }

fitVCA <- function(form, Data, method=c("anova", "reml"), scale=TRUE, VarVC=TRUE, ...)
{
	call <- match.call()
	method <- match.arg(tolower(method[1]), choices=c("anova", "reml"))
	
	if(scale)
	{
		if(method=="anova")
			fit <- Scale("anovaVCA", form=form, Data=Data, ...)
		else
			fit <- Scale(remlVCA, form=form, Data=Data, VarVC=VarVC, ...)
	}
	else
	{
		if(method=="anova")
			fit <- anovaVCA(form=form, Data=Data, ...)
		else
			fit <- remlVCA(form=form, Data=Data, VarVC=VarVC, ...)	
	}
	
	wasVCA <- FALSE								# record that result was VCA-object
	
	if(class(fit) == "VCA")						# if by-processing was not used 			
	{
		fit <- list(fit)
		wasVCA <- TRUE
	}
	
	for(i in 1:length(fit))								# over each list-element
	{
		if(scale)
			fit[[i]] <- reScale(fit[[i]], VarVC=VarVC)				# first re-scale
		
		fit[[i]]$call 			<- call								# prevent non-informative object-name, e.g. "form"
		fe 						<- fixef(fit[[i]])
		X  						<- getMat(fit[[i]], "X")		
		fit[[i]]$fitted.values	<- as.numeric(as.matrix(X) %*% fe[,"Estimate", drop=F])
		fixed 					<- fit[[i]]$fixed[!grepl(":", fit[[i]]$fixed)]	# no interactions
		fit[[i]]$xlevels 		<- lapply(fixed, function(x) as.character(unique(fit[[i]]$data[,x])))
		names(fit[[i]]$xlevels) <- fixed
		
		attr(fit[[i]]$Type, "fitted.by") <- "fitVCA"	# only LMM fitted by 'fitLMM' can be handled in standard ways
	}
	
	if(wasVCA)								
		fit <- fit[[1]]
	
	fit
}


#' Load 'RevoUtilsMath'-package if available.
#' 
#' This function is taken from the Rprofile.site file of Microsoft R Open.
#' It was added to the package namespace to avoid a NOTE during the R CMD check
#' process stating that this function is not gobally defined.
#' 
#' Only change to the original version is a different bracketing scheme to match
#' the one used in the remaining source-code of the package. 
#' 
#' @param package		(character) package name to load, usually this will be package
#' 						'RevoUtilsMath' if available
#' 
#' @author Authors of the Rprofile.site file in Microsoft R Open.

load_if_installed <- function(package) 
{
	if (!identical(system.file(package="RevoUtilsMath"), "")) 
	{
		do.call('library', list(package))
		return(TRUE)
	} 
	else
		return(FALSE)
}

#' Scale Response Variable to Ensure Robust Numerical Calculations.
#' 
#' Function determines scaling factor for transforming the mean of the response to a range
#' between 0.1 and 1, applies scaling of the response and binds the scaling factor to the 
#' data as attribute 'scale'.
#' 
#' @param Data			(data.frame) with the data to be fitted and the response to be scaled
#' @param resp			(character) name of the (numeric) response variable
#' 
#' @return (data.frame) with the response scaled according to the scaling-factor,
#' 			which is recorded in the attribute \code{scale} of the data set
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}

scaleData <- function(Data=NULL, resp=NULL)
{
	Mean 	<- mean(Data[,resp], na.rm=TRUE)								# scale response variable for more robust numerical optimization
	scale 	<- 10^ceiling(log10(abs(mean(Data[,resp], na.rm=TRUE))))	
	Data[,resp] <- Data[, resp] / scale
	attr(Data, "scale") <- scale
	Data
}

#' Automatically Scale Data Calling These Functions: 'anovaVCA', 'anovaMM', 'remlVCA' or 'remlMM'.
#' 
#' This function scales data before fitting a linear mixed model aiming to avoid numerical problems
#' when numbers of the response variable are either very small or very large. It adds attribute "scale" 
#' to the resulting 'VCA'-object, which is used by function \code{\link{reScale}} to transform back the
#' VCA-results of a \code{VCA} or \code{VCAinference} object that was previously scaled.
#' 
#' NOTE: Scaling is applied on the complete data set, without checking whether there are incomplete
#' observations or not!
#' 
#' @param Fun		(expr, function, character) either a complete function call to one of "anovaVCA", "anovaMM", "remlVCA", "remlMM",
#' 					a character string or just the function name without quotes (see example)
#' @param form		(formula) specifying the model to fitted by 'Fun'
#' @param Data		(data.frame) with all variables specified via 'Fun'
#' @param ...		additional arguments applying to one of the four functions \code{\link{anovaVCA}},\code{\link{anovaMM}},
#' 					\code{\link{remlVCA}}, \code{\link{remlMM}}
#' 
#' @return (object) of class 'VCA' which can be used as input for function \code{\link{VCAinference}} 
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @seealso \code{\link{reScale}}
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_3)
#' 
#' # simulate very large numbers of the response
#' dat3   <- dataEP05A2_3
#' dat3$y <- dat3$y * 1e8
#' 
#' # now try to fit 21-day model to this data
#' fit <- anovaVCA(y~day/run, dat3)
#' 
#' # now use 'Scale' function
#' fit1 <- Scale("anovaVCA", y~day/run, dat3)
#' fit2 <- Scale(anovaVCA, y~day/run, dat3)	# also works
#' fit3 <- Scale(anovaVCA(y~day/run, dat3)) # works as well
#' 
#' # back to original scale
#' (fit1 <- reScale(fit1))
#' (fit2 <- reScale(fit2))
#' (fit3 <- reScale(fit3))
#' 
#' # reference values
#' fit0 <- anovaVCA(y~day/run, dataEP05A2_3, MME=TRUE)
#' inf0 <- VCAinference(fit0, VarVC=TRUE)
#' 
#' fit1 <- Scale(anovaVCA(y~day/run, dataEP05A2_3, MME=TRUE))
#' inf1 <- VCAinference(fit1, VarVC=TRUE)
#' inf1 <- reScale(inf1)
#' 
#' # compare to reference
#' print(inf0, what="VC")
#' print(inf1, what="VC")
#' print(inf0, what="SD")
#' print(inf1, what="SD")
#' print(inf0, what="CV")
#' print(inf1, what="CV")
#' 
#' # now use REML-based estimation
#' fit0 <- remlVCA(y~day/run, dataEP05A2_3)
#' inf0 <- VCAinference(fit0)
#' 
#' fit1 <- Scale("remlVCA", y~day/run, dataEP05A2_3)
#' inf1 <- VCAinference(fit1)
#' inf1 <- reScale(inf1)
#' 
#' # compare to reference
#' print(inf0, what="VC")
#' print(inf1, what="VC")
#' print(inf0, what="SD")
#' print(inf1, what="SD")
#' print(inf0, what="CV")
#' print(inf1, what="CV")
#' 
#' # scaling also works with by-processing
#' data(VCAdata1)
#' fit <- Scale(anovaVCA(y~(device+lot)/day/run, VCAdata1, by="sample"))
#' reScale(fit)
#' }

Scale <- function(Fun, form, Data, ...)
{
	call <- as.list(match.call())			# check whether a complete funtion call was the sole argument
	
	if(class(call$Fun) == "call")
	{
		call <- as.list(call$Fun)
		
		Fun  <- as.character(call[[1]])
		form <- as.formula(call[[2]])
		Data <- get(as.character(call[[3]]))
		if(length(call) > 3)
			Args <- call[4:length(call)]
		else
			Args <- NULL
	}
	else	
		Args <- list(...)	
	
	if("by" %in% names(Args))
	{
		stopifnot(is.character(Args$by))
		stopifnot(Args$by %in% colnames(Data))
		stopifnot(is.factor(Data[,Args$by]) || is.character(Data[,Args$by]))
		
		levels  	<- unique(Data[,Args$by])
		tmpArgs 	<- Args
		tmpArgs$by	<- NULL
		if(length(tmpArgs) == 0)
			tmpArgs <- NULL
		res 		<- lapply(levels, function(x) Scale(Fun=Fun, form=form, Data[Data[,Args$by] == x,], tmpArgs))
		names(res) 	<- paste0(Args$by, levels)
		return(res)
	}
	if(is.function(Fun))						# if a function was passed and not the name of it
		Fun <- deparse(substitute(Fun))
	fun  <- match.arg(Fun, choices=c("anovaVCA", "anovaMM", "remlVCA", "remlMM"))
	stopifnot(class(form) == "formula")
	stopifnot(is.data.frame(Data))
	
	tobj <- terms(form)
	stopifnot(attr(tobj, "response")==1)
	resp		<- rownames(attr(tobj, "factors"))[1]
	
	Mean  		<- mean(Data[,resp], na.rm=TRUE)								# scale response variable for more robust numerical optimization
	scale 		<- 10^ceiling(log10(abs(mean(Data[,resp], na.rm=TRUE))))	
	Data[,resp] <- Data[,resp]/(scale)
	
	FunArgs <- list(form=form, Data=Data)
	FunArgs <- c(FunArgs, Args)
	
	obj <- do.call(Fun, FunArgs)
	
	obj$scale <- scale
	obj
}

#' Re-Scale results of 'VCA' or 'VCAinference' objects which were constructed using
#' scaled values of the response variable.
#' 
#' Function adjusts variance components (VC) and standard deviations (SD) and their respective
#' confidence intervals of 'VCAinference' objects, and the 'VCAobj' sub-element. For 'VCA' objects
#' the VC and SD values are adjusted as well as the fixed and random effects and the covariance-matrix
#' of fixed effects.
#' 
#' @param obj		(object) either of class 'VCA' or 'VCAinference'
#' @param VarVC		(logical) TRUE = variance-covariance matrix of the fitted model 'obj'
#' 					will be computed and automatically re-scaled, FALSE = variance-covariance
#' 					matrix will not be computed and re-scaled. This might cause wrong results
#' 					in downstream analyses which require this matrix on the correct scale! Only
#' 					use this option if computation time really matters!
#' 
#' @return (object) either of class 'VCA' or 'VCAinference', where results have been 
#'         transformed back to the original scale of the response variable 
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @seealso \code{\link{Scale}}
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_3)
#' 
#' # reference values
#' fit0 <- anovaVCA(y~day/run, dataEP05A2_3, MME=TRUE)
#' inf0 <- VCAinference(fit0, VarVC=TRUE)
#' 
#' fit1 <- Scale("anovaVCA", y~day/run, dataEP05A2_3, MME=TRUE)
#' inf1 <- VCAinference(fit1, VarVC=TRUE)
#' inf1 <- reScale(inf1)
#' 
#' # compare to reference
#' print(inf0, what="VC")
#' print(inf1, what="VC")
#' print(inf0, what="SD")
#' print(inf1, what="SD")
#' print(inf0, what="CV")
#' print(inf1, what="CV")
#' 
#' # now use REML-based estimation
#' fit0 <- remlVCA(y~day/run, dataEP05A2_3)
#' inf0 <- VCAinference(fit0)
#' 
#' fit1 <- Scale("remlVCA", y~day/run, dataEP05A2_3, MME=TRUE)
#' inf1 <- VCAinference(fit1)
#' inf1 <- reScale(inf1)
#' 
#' # compare to reference
#' print(inf0, what="VC")
#' print(inf1, what="VC")
#' print(inf0, what="SD")
#' print(inf1, what="SD")
#' print(inf0, what="CV")
#' print(inf1, what="CV")
#' }

reScale <- function(obj, VarVC=TRUE)
{
	if(is.list(obj) && !class(obj) %in% c("VCA", "VCAinference"))
	{
		
		if(!all(sapply(obj, class) %in% c("VCA", "VCAinference")))
			stop("Only lists of 'VCA' or 'VCAinference' objects are accepted!")
		
		obj.len <- length(obj)
		
		res <- lapply(obj, FUN=reScale)
		names(res) <- names(obj)
		
		if(obj.len == 1)			# mapply returns a list of length 2 in case that length(obj) was equal to 1
			res <- res[1]
		
		return(res)
	}	
	
	if(!is.null(obj$rescaled))		# re-scaling has already been done on this object
		return(obj)
	
	stopifnot(class(obj) %in% c("VCA", "VCAinference"))
	
	if(class(obj) == "VCAinference")
		VCAobj <- obj$VCAobj
	else
		VCAobj <- obj
	
	scale <- VCAobj$scale	
	
	if(VarVC)						# estimate covariance-matrix of VCs before any re-scaling takes place
	{
		VCAobj <- solveMME(VCAobj)
		VCAobj$VarCov <- vcovVC(VCAobj)
	}	
	
	VCAobj$aov.tab[,"VC"] 		<- VCAobj$aov.tab[,"VC"] * scale^2
	VCAobj$aov.tab[,"SD"] 		<- VCAobj$aov.tab[,"SD"] * scale
	if("Var(VC)" %in% colnames(VCAobj$aov.tab))
		VCAobj$aov.tab[,"Var(VC)"]	<- VCAobj$aov.tab[,"Var(VC)"] * scale^4
	if(VCAobj$EstMethod == "ANOVA")
	{
		VCAobj$aov.tab[,"SS"]	<- VCAobj$aov.tab[,"SS"] * scale^2
		VCAobj$aov.tab[,"MS"]	<- VCAobj$aov.tab[,"MS"] * scale^2
	}
	VCAobj$Mean 				<- VCAobj$Mean * scale
	
	VCAobj$Matrices$y <- VCAobj$Matrices$y * scale
	VCAobj$data[,VCAobj$response] <- VCAobj$data[,VCAobj$response] * scale			# re-scale response variable in the data-element
	
	if(!is.null(VCAobj$FixedEffects))
		VCAobj$FixedEffects 	<- try(VCAobj$FixedEffects * scale,  silent=TRUE)
	if(!is.null(VCAobj$RandomEffects))
		VCAobj$RandomEffects 	<- try(VCAobj$RandomEffects * scale, silent=TRUE)
	if(!is.null(VCAobj$VarFixed))
		VCAobj$VarFixed 		<- try(VCAobj$VarFixed * scale^2, 	 silent=TRUE)
	if(!is.null(VCAobj$VarCov))
		VCAobj$VarCov			<- try(VCAobj$VarCov * scale^4, 	 silent=TRUE)
	
	if(class(obj) == "VCAinference")
	{
		obj$VCAobj <- VCAobj
		
		# variance components
		
		obj$ConfInt$VC$OneSided$LCL <- obj$ConfInt$VC$OneSided$LCL * scale^2
		obj$ConfInt$VC$OneSided$UCL <- obj$ConfInt$VC$OneSided$UCL * scale^2
		
		obj$ConfInt$VC$TwoSided$LCL <- obj$ConfInt$VC$TwoSided$LCL * scale^2
		obj$ConfInt$VC$TwoSided$UCL <- obj$ConfInt$VC$TwoSided$UCL * scale^2
		
		# SD
		
		obj$ConfInt$SD$OneSided$LCL <- obj$ConfInt$SD$OneSided$LCL * scale
		obj$ConfInt$SD$OneSided$UCL <- obj$ConfInt$SD$OneSided$UCL * scale
		
		obj$ConfInt$SD$TwoSided$LCL <- obj$ConfInt$SD$TwoSided$LCL * scale
		obj$ConfInt$SD$TwoSided$UCL <- obj$ConfInt$SD$TwoSided$UCL * scale
	}
	else
		obj <- VCAobj
	
	obj$rescaled <- TRUE			# record that re-scaling was performed
	obj 
}


#' Check for Availability of Intel's Math Kernel Library.
#' 
#' Majority of the code is borrowed from the Microsoft R Open Rprofile.site file.
#' In case MKL can be detected this information will be stored in a separate envrionment, which
#' other function know about. If so, an optimized version of function \code{\link{getGB}}
#' will be used which used ordinary matrix-objects instead of matrices defined by the
#' \code{Matrix}-package. This seems to accelerate computation time for large datasets
#' by up to factor 30.
#' 
#' This function is for internal use only and therefore not exported.
#' 
#' @return variable 'MKL' in envir "msgEnv" will be set to TRUE/FALSE
#' 
#' @author 	Authors of the Rprofile.site file in Microsoft R Open,
#' 			Andre Schuetzenmeister \email{[email protected]@roche.com}

check4MKL <- function()
{
	if("msgEnv" %in% ls(.GlobalEnv) && !is.null(msgEnv$MKL))
		return(msgEnv$MKL)
	else
	{
		msgEnv <<- new.env(parent=emptyenv())				
		
		if(Sys.info()["sysname"] == "Darwin")				# Mac OSx
		{		
			assign("MKL", TRUE, envir=msgEnv)				# set MKL to TRUE, although, it may not be installed --> no good method known to check for MKL under MacOS
		} 
		else 												# other operating systems
		{		
			MRO <- FALSE
			try(MRO <- load_if_installed("RevoUtilsMath"), silent=TRUE)			# function only exists in MRO environment
			if(!inherits(MRO, "try-error") && MRO)
				assign("MKL", TRUE, envir=msgEnv)
			else 
				assign("MKL", FALSE, envir=msgEnv)
		}
		
		return(msgEnv$MKL)
	}
}


#' Giesbrecht & Burns Approximation of the Variance-Covariance Matrix of Variance Components.
#'
#' Compute variance covariance matrix of variance components of a linear mixed model
#' via the method stated in Giesbrecht and Burns (1985). 
#' 
#' This function is not intended to be called by users and therefore not exported.
#'
#' @param obj		(object) with list-type structure, e.g. \code{VCA} object fitted by ANOVA
#' 					or a premature \code{VCA} object fitted by REML
#' @param tol		(numeric) values < 'tol' will be considered being equal to zero
#' 
#' @return 	(matrix) corresponding to the Giesbrecht & Burns approximation
#'			of the variance-covariance matrix of variance components
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com},
#' 		   Florian Dufey \email{[email protected]@contractors.roche.com}
#' 
#' @seealso \code{\link{vcovVC}}, \code{\link{remlVCA}}, \code{\link{remlMM}}
#'
#' @references
#' Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York
#' 
#' Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample
#' Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486 
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_3)
#' fit <- anovaVCA(y~day/run, dataEP05A2_3)
#' fit <- solveMME(fit)		# some additional matrices required
#' getGB(fit)
#' }

getGB <- function (obj, tol = 1e-12)
{
	Nvc <- obj$Nvc
	Z <- obj$Matrices$Zre
	Q <- obj$Matrices$Q
	if (is.null(Q)) 
	{
		X <- getMat(obj, "X")
		Vi <- getMat(obj, "Vi")
		T <- getMat(obj, "T")
		Q <- Vi - Vi %*% X %*% T
	}
	VCvar <- matrix(0, Nvc, Nvc)
	nze <- NULL
	for (i in 1:Nvc) 
	{
		VCi <- obj$VCoriginal[i]
		if (is.null(VCi))
			VCi <- obj$aov.tab[obj$re.assign$terms[i], "VC"]
		if (i < Nvc && abs(VCi) < tol) 
		{
			VCvar[i, ] <- VCvar[, i] <- 0
			next
		}
		else
			nze <- c(nze, i)
		Zi  <- Z[, which(obj$re.assign$ind == i)]
		ZiT <- t(Zi)
		
		for (j in i:Nvc) 
		{			
			Zj    <- Z[, which(obj$re.assign$ind == j)]
			ZjT   <- t(Zj)
			Zpart <- as.matrix(ZjT %*% Q %*% Zi)
			
			if(j == Nvc)
			{
				if(j == i)
					VCvar[i,j] <- sum(sum(Q*Q))
				else
				{
					ZiTQ <- as.matrix(ZiT %*% Q)
					VCvar[i,j] <- VCvar[j,i] <- sum(sum(ZiTQ*ZiTQ))		# trace A*t(B)
				}
			}
			else
				VCvar[j,i] <- VCvar[i,j] <- sum(sum(Zpart*Zpart))		# trace A*t(B)
		}
	}
	nze <- unique(nze)
	VCnames <- obj$VCnames
	if (!"error" %in% VCnames)
		VCnames <- c(VCnames, "error")
	rownames(VCvar) <- colnames(VCvar) <- VCnames
	
	VCvar[nze, nze] <- 2 * solve(VCvar[nze, nze])
	VCvar <- VCvar
	attr(VCvar, "method") <- "gb"
	VCvar
}

#' Construct Variance-Covariance Matrix of Random Effects for Models Fitted by Function 'lmer'.
#' 
#' This function restricts the variance-covariance matrix of random effects \eqn{G} to be either
#' diagonal ('cov=FALSE') or to take any non-zero covariances into account (default, 'cov=TRUE').
#' 
#' This function is not intended to be called directly by users and therefore not exported!
#' 
#' @param obj		(object) inheriting from class 'lmerMod'
#' @param cov		(logical) TRUE = in case of non-zero covariances a block diagonal matrix will be constructed,
#'                  FALSE = a diagonal matrix with all off-diagonal element being equal to zero will be contructed
#' 
#' @return (Matrix) representing the variance-covariance structure of random effects \eqn{G}
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples 
#' \dontrun{
#' data(Orthodont)
#' Ortho <- Orthodont
#' Ortho$age2 <- Ortho$age - 11
#' Ortho$Subject <- factor(as.character(Ortho$Subject))
#' fit <-lmer(distance~Sex+Sex:age2+(age2|Subject), Ortho) 
#' G1 <- VCA:::lmerG(fit, cov=FALSE)
#' G2 <- VCA:::lmerG(fit, cov=TRUE)
#' G1[1:10,1:10]
#' G2[1:10,1:10]
#' }

lmerG <- function(obj, cov=FALSE)
{
	stopifnot(inherits(obj, "lmerMod"))
	
	vc  <- VarCorr(obj)
	Nre <- unlist(lapply(lme4::ranef(obj), nrow))					# number of random effects per variance component 
	
	lst <- list()
	for(i in 1:length(Nre))
	{
		tmp.vc <- vc[[i]]
		
		if(!cov && nrow(tmp.vc) > 1)
		{
			tmp.vc <- diag(diag(tmp.vc))
		}
		
		lst <- c(lst, replicate(Nre[i], tmp.vc, simplify=FALSE))	# remove covariances from off-diagonal
		
	}
	G <- bdiag(lst)
	G
}


#' Derive and Compute Matrices for Objects Fitted by Function 'lmer'.
#' 
#' Function derives and computes all matrices required for down-stream
#' analyses of VCA-objects fitted with REML via function \code{\link{lmer}}.
#' 
#' Mixed Model Equations (MME) are solved for fixed and random effects applying the same
#' constraints as in \code{\link{anovaMM}}. 
#' The most elaborate and therefore time consuming part is to prepare all matrices required for 
#' approximating the variance-covariance matrix of variance components (see \code{\link{getGB}}).
#' To reduce the computational time, this function tries to optimize object-classes depending
#' on whether Intel's (M)ath (K)ernel (L)ibrary could be loaded or not. MKL appears to be more
#' performant with ordinary matrix-objects, whereas all other computations are perfomred using
#' matrix-representations of the \code{Matrix}-package.
#' 
#' This function is not intended to be called directly by users and therefore not exported.
#' 
#' @param obj		(object) inheriting from 'lmerMod'
#' @param tab		(data.frame) representing the basic VCA-table
#' @param terms		(character) vector used for ordering variance components
#' @param cov		(logical) take non-zero covariances among random effects into account (TRUE) or
#' 					not (FALSE), the latter is the default in this package and also implemented in
#' 					\code{\link{remlVCA}}, \code{\link{anovaVCA}}, and \code{\link{anovaMM}}.
#' @param X			(matrix) design matrix of fixed effects as constructed to meet VCA-package requirements
#' 
#' @return (list), a premature 'VCA' object
#' 
#' @seealso \code{\link{remlVCA}}, \code{\link{remlMM}}
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}

lmerMatrices <- function(obj, tab=NULL, terms=NULL, cov=FALSE, X=NULL)
{
	stopifnot(inherits(obj, "lmerMod"))
	
	re.org <- re <- lme4::ranef(obj)	# use lme4's ranef S3-method
	
	VCnam  <- NULL
	reInd  <- list()
	last   <- 0
	count  <- 1
	REnam  <- names(re)
	REnamZ <- NULL
	
	for(i in 1:length(re))				# transform random effects into VCA-compatible structure and derive					
	{									# column-index in random effects design matrix Z for all random effects
		if(ncol(re[[i]]) > 1)			# regression model, i.e. random effects in multi-column matrix
		{
			REi 	<- re[[i]]
			NRi 	<- nrow(REi)
			NCi 	<- ncol(REi)
			trm 	<- names(re)[i]
			tmp.re 	<- ind <- NULL			
			
			for(j in 1:NCi)
			{
				reInd[[count]] <- seq(last+1, last+NRi*NCi, by=NCi)		# columns in Z corresponding to a specific random effect
				
				if(j == 1)
				{
					tmp.nam <- paste(trm, rownames(re[[i]]), sep="")
					names(reInd)[count] <- REnam[count]
				}
				else
				{
					tmp.nam <- c(tmp.nam, paste(colnames(REi)[j], tmp.nam, sep=":"))
					names(reInd)[count] <- paste(colnames(REi)[j], REnam[i], sep=":")
				}
				last  <- last + 1
				count <- count + 1
				
				if(j == 1)
					ind <- eval(parse(text=paste("((",j,"-1)*",NRi,"+1):(", j,"*", NRi, ")", sep="")))
				else
					ind <- paste(ind, eval(parse(text=paste("((",j,"-1)*",NRi,"+1):(", j,"*", NRi, ")", sep=""))), sep=", ")
			}
			ind 	<- paste("c(", paste(ind, collapse=","), ")", sep="")
			ind 	<- eval(parse(text=ind))
			REnamZ 	<- c(REnamZ, tmp.nam[ind])							# column names in Z-matrix			
			cn 		<- colnames(REi)
			cn 		<- paste(names(re)[i], cn, sep=":")
			cn 		<- gsub(":\\(Intercept\\)", "", cn)					# get rif of ()
			VCnam 	<- c(VCnam, cn)										# names of Variance components
		}
		else		# single column random effects matrix
		{
			reInd[[count]] <- (last+1):(last+nrow(re[[i]]))
			names(reInd)[count] <- REnam[count]
			last   <- last + nrow(re[[i]])
			count  <- count + 1
			nami   <- unlist(strsplit(REnam[i], ":"))
			rni    <- rownames(re[[i]])
			VCnam  <- c(VCnam, REnam[i])			
			REnamZ <- c(REnamZ, sapply(rni, function(x) paste(paste(nami, unlist(strsplit(x, ":")), sep=""), collapse=":")))
		}
	}
	REnamZ 		<- as.character(REnamZ)						# names for the random effects and the columns in Z
	reInd  		<- reInd[terms]								# order according to order of variables in user-specified formula
	reInd  		<- reInd[which(!is.na(names(reInd)))]
	Nre	   		<- as.numeric(unlist(lapply(re.org, nrow)))
	sNre   		<- sum(Nre)
	re.assign 	<- list(ind=integer(sNre), terms=names(reInd))
	VC	  	  	<- tab[-c(1, nrow(tab)), "VC"]
	
	for(i in 1:length(reInd))
		re.assign$ind[reInd[[i]]] <- i
	
	Nvc	  <- nrow(tab)-1							# minus total and error								
	Zt 	  <- obj@pp$Zt
	
	if(check4MKL())
		Zt <- as.matrix(Zt)
	
	if(is.null(X) || class(X) != "matrix")
		X <- model.matrix(obj, type="fixed")
	
	Xt <- t(X)
	Z  <- t(Zt)
	colnames(Z) <- REnamZ
	G  <- lmerG(obj, cov=cov)						# construct G-matrix
	
	if(check4MKL())
		R	<- diag(nrow(obj@frame))* tab[nrow(tab),"VC"]
	else
		R	<- Diagonal(nrow(obj@frame))* tab[nrow(tab),"VC"]
	
	V	 <- Z %*% G %*% Zt + R						# variance-covariance matrix of observations
	
	Vi	 <- solve(V)
	ViX  <- Vi %*% X
	XtVi <- Xt %*% Vi
	P 	 <- MPinv(Xt %*% ViX)   					# re-compute since 'vcov(obj)' differs from e.g. SAS PROC MIXED
	Q	 <- Vi - ViX %*% P %*% XtVi
	T    <- P %*% XtVi
	
	res 		 <- list()
	res$Nvc 	 <- Nvc
	res$VCnames  <- c(terms, "error") 
	y  <- obj@resp$y
	fe <- P %*% XtVi %*% y						# solve mixed model equations for fixed effects
	
	colnames(fe) <- "Estimate"
	rownames(fe) <- colnames(X)
	iind <- which(names(fe) == "(Intercept)")
	
	if(length(iind) > 0)
		names(fe)[iind] <- "int"
	
	res$Matrices 		<- list(Zre=Z, G=G, R=R, V=V, Vi=Vi, Q=Q, X=X, T=T, y=y)
	res$FixedEffects  	<- fe	
	re  				<- G %*% Zt %*% Vi %*% (y - X %*% fe) 		# estimate random effects solving Mixed Model Equations
	re 					<- Matrix(re)								# re-order random effects accordingly
	rownames(re) 		<- REnamZ
	res$RandomEffects 	<- re
	res$re.assign 		<- re.assign
	res$VarFixed  		<- P
	res$ColOrderZ 		<- reInd			# keep information about original col-indexing
	res
}

#' Derive VCA-Summary Table from an Object Fitted via Function \code{\link{lmer}}.
#'
#' This function builds a variance components analysis (VCA) table
#' from an object representing a model fitted by function \code{\link{lmer}}
#' of the \code{lme4} R-package. 
#' 
#' It applies the approximation of the variance-covariance
#' matrix of variance components according to Giesbrecht & Burns (1985) and uses this
#' information to approximate the degrees of freedom according to Satterthwaite
#' (see SAS PROC MIXED documentation option 'CL').
#' 
#' This function can be used to create a VCA-results table from almost any fitted 'lmerMod'-object, i.e. one can
#' apply it to a model fitted via function \code{\link{lmer}} of the \code{lme4}-package. The only 
#' additional argument that needs to be used is 'tab.only' (see examples).
#' 
#' @param obj		(lmerMod) object as returned by function \code{\link{lmer}}
#' @param VarVC		(logical) TRUE = the variance-covariance matrix of variance components will be approximated
#' 					following the Giesbrecht & Burns approach, FALSE = it will not be approximated	
#' @param terms		(character) vector, optionally defining the order of variance terms to be used
#' @param Mean		(numeric) mean value used for CV-calculation
#' @param cov		(logical) TRUE = in case of non-zero covariances a block diagonal matrix will be constructed,
#'                  FALSE = a diagonal matrix with all off-diagonal elements being equal to zero will be contructed
#' @param X			(matrix) design matrix of fixed effects as constructed to meet VCA-package requirements
#' @param tab.only	(logical) TRUE = will return only the VCA-results table as 'data.frame', argument 'VarVC' will 
#' 					automatically set to 'FALSE' (see details)
#' 
#' @return (list) still a premature 'VCA'-object but close to a "complete" 'VCA'-object
#' 
#' @seealso \code{\link{remlVCA}}, \code{\link{remlMM}}
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#'
#' @references
#' Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York
#' 
#' Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample
#' Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486 
#' 
#' @examples 
#' \dontrun{
#' # fit a model with a VCA-function first
#' data(VCAdata1)
#' fit0 <- remlVCA(y~(device+lot)/day/run, subset(VCAdata1, sample==5))
#' fit0
#' 
#' # fit the same model with function 'lmer' of the 'lme4'-package
#' library(lme4)
#' fit1 <- lmer(y~(1|device)+(1|lot)+(1|device:lot:day)+(1|device:lot:day:run),
#'              subset(VCAdata1, sample==5))
#' lmerSummary(fit1, tab.only=TRUE)
#' }

lmerSummary <- function(obj, VarVC=TRUE, terms=NULL, Mean=NULL, cov=FALSE, X=NULL, tab.only=FALSE)
{
	stopifnot(inherits(obj, "lmerMod"))	
	if(tab.only)
	{
		VarVC <- FALSE
		
		if(is.null(Mean))
			Mean <- mean(obj@resp$y, na.rm=TRUE)
	}
	
	Sum  <- as.data.frame(summary(obj)$varcor)
	
	Sum[which(Sum[,"var1"] %in% c(NA, "(Intercept)")), "var1"] <- ""
	
	Sum[,"grp"] <- apply(Sum[,c("grp", "var1")], 1, 
			function(x)
			{
				if(x[2] == "")
					return(x[1])
				else
					return(paste(rev(x), collapse=":"))
			}
	)
	re.cor <- NULL
	
	if(any(!is.na(Sum[,"var2"])))
	{
		ind.cor <- which(!is.na(Sum[,"var2"]))
		re.cor  <- Sum[ind.cor,,drop=FALSE]
		Sum <- Sum[-ind.cor,]
	}
	
	if(tab.only)
	{
		if(is.null(terms))
		{
			terms <- Sum[,"grp"]
			terms <- terms[-length(terms)]
		}
	}
	
	Sum  <- Sum[,-c(2,3)]
	rownames(Sum) <- Sum[,"grp"]
	Sum <- Sum[c(terms, "Residual"), ]
	colnames(Sum) <- c("Name", "VC", "SD")
	Sum <- rbind(c(Name="total", VC=sum(Sum[,"VC"]), SD=sqrt(sum(Sum[,"VC"]))), Sum)
	rownames(Sum) <- Sum[,"Name"]
	Sum$VC <- as.numeric(Sum$VC)
	Sum$SD <- as.numeric(Sum$SD)
	Sum$Perc <- 100*Sum$VC/Sum$VC[1]
	Sum$CV	 <- 100*Sum$SD/Mean
	
	if(!tab.only)										# complete VCA-object shall be created
	{
		obj <- lmerMatrices(obj, tab=Sum, terms=terms,	# compute some required matrices
				cov=cov, X=X)							
		obj$aov.tab <- Sum								# required for Giesbrecht & Burns approximation
	}
	
	if(VarVC)
	{
		varVC <- obj$VarCov <- getGB(obj)				# apply Giesbrecht & Burns approximation optimized for MKL
		varVC <- diag(varVC)
		varVC <- c(sum(obj$VarCov), varVC)				# variance of total is sum of all elements of the variance-covariance matrix
		seVC  <- sqrt(varVC)
		Sum$varVC <- varVC
		Sum$Wald <- Sum$VC/seVC							# Wald-statistic
		Sum$DF	 <- 2*Sum$Wald^2						# see SAS-Doc of PROC MIXED
		Sum 	 <- Sum[,c(8,2,4,3,5,6)]
		colnames(Sum)[c(3,5,6)] <- c("%Total", "CV[%]", "Var(VC)")
	}
	else
	{	
		Sum <- Sum[,c(2,4,3,5)]
		colnames(Sum)[c(2,4)] <- c("%Total", "CV[%]")
	}
	
	if(check4MKL())			# remaining part of the package computes with Matrix-package
	{
		obj$Matrices$Zre 	<- Matrix(obj$Matrices$Zre)
		obj$Matrices$G   	<- Matrix(obj$Matrices$G)
		obj$Matrices$R 		<- Matrix(obj$Matrices$R)
		obj$Matrices$V 		<- Matrix(obj$Matrices$V)
		obj$Matrices$Vi 	<- Matrix(obj$Matrices$Vi)
		obj$Matrices$Q 		<- Matrix(obj$Matrices$Q)
		obj$Matrices$X	    <- Matrix(obj$Matrices$X)
		obj$Matrices$T	    <- Matrix(obj$Matrices$T)
	}	
	
	rownames(Sum)[nrow(Sum)] <- "error"
	if(tab.only)							# return just the VCA-table
	{
		Sum <- as.matrix(Sum)
		attr(Sum, "Mean") <- Mean
		return(Sum)
	}
	obj$aov.tab <- Sum
	obj$re.cor <- re.cor					# save correlation among random terms
	obj
}


#' Perform (V)ariance (C)omponent (A)nalysis via REML-Estimation.
#' 
#' Function performs a Variance Component Analysis (VCA) using Restricted Maximum Likelihood (REML)
#' to fit the random model, i.e. a linear mixed model (LMM) where the intercept is the only fixed effect.
#' 
#' Here, a variance component model is fitted by REML using the \code{\link{lmer}} function of the
#' \code{lme4}-package. For all models the Giesbrechnt & Burns (1985) approximation of the variance-covariance
#' matrix of variance components (VC) is applied. A Satterthwaite approximation of the degrees of freedom
#' for all VC and total variance is based on this approximated matrix using \eqn{df=2Z^2}, where
#' \eqn{Z} is the Wald statistic \eqn{Z=\sigma^2/se(\sigma^2)}, and \eqn{\sigma^2} is here used for an
#' estimated variance. The variance of total variability, i.e. the sum of all VC is computed via summing
#' up all elements of the variance-covariance matrix of the VC.
#' Note, that for large datasets approximating the variance-covariance matrix of VC is computationally expensive
#' and may take very long. There is no Fisher-information matrix available for 'merMod' objects, which can
#' serve as approximation. To avoid this time-consuming step, use argument 'VarVC=FALSE' but remember,
#' that no confidence intervals for any VC will be available. If you use Microsoft's R Open, formerly known
#' as Revolution-R, which comes with Intel's Math Kernel Library (MKL), this will be automatically detected
#' and an environment-optimized version will be used, reducing the computational time very much (see examples).
#' 
#' @param form          (formula) specifying the model to be fit, a response variable left of the '~' is mandatory
#' @param Data          (data.frame) containing all variables referenced in 'form'
#' @param by			(factor, character) variable specifying groups for which the analysis should be performed individually,
#' 						i.e. by-processing
#' @param VarVC			(logical) TRUE = the variance-covariance matrix of variance components will be approximated using 
#' 						the method found in Giesbrecht & Burns (1985), which also serves as basis for applying a Satterthwaite
#' 						approximation of the degrees of freedom for each variance component, FALSE = leaves out this step, 
#' 						no confidence intervals for VC will be available
#' @param quiet			(logical) TRUE = will suppress any messages or warnings, which will be issued otherwise
#' @param order.data	(logical) TRUE = class-variables will be ordered increasingly, FALSE = ordering of class-variables
#' 						will remain as is 
#' 
#' @seealso \code{\link{remlMM}}, \code{\link{VCAinference}}, \code{\link{ranef.VCA}}, \code{\link{residuals.VCA}},
#' 			\code{\link{anovaVCA}}, \code{\link{anovaMM}}, \code{\link{plotRandVar}}, \code{\link{lmer}}
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples
#' \dontrun{
#' 
#' # a VCA standard example
#' data(dataEP05A2_3)
#' 
#' # fit it by ANOVA first, then by REML
#' fit0 <- anovaVCA(y~day/run, dataEP05A2_3) 
#' fit1 <- remlVCA(y~day/run, dataEP05A2_3)
#' fit0
#' fit1
#'  
#' # make example unbalanced
#' set.seed(107)
#' dat.ub <- dataEP05A2_3[-sample(1:80, 7),]
#' fit0ub <- anovaVCA(y~day/run, dat.ub) 
#' fit1ub <- remlVCA(y~day/run, dat.ub) 
#' 
#' # not that ANOVA- and REML-results now differ
#' fit0ub
#' fit1ub
#' 
#' ### Use the six sample reproducibility data from CLSI EP5-A3
#' ### and fit per sample reproducibility model
#' data(CA19_9)
#' fit.all <- remlVCA(result~site/day, CA19_9, by="sample")
#' 
#' reproMat <- data.frame(
#'  Sample=c("P1", "P2", "Q3", "Q4", "P5", "Q6"),
#'  Mean= c(fit.all[[1]]$Mean, fit.all[[2]]$Mean, fit.all[[3]]$Mean, 
#' 	        fit.all[[4]]$Mean, fit.all[[5]]$Mean, fit.all[[6]]$Mean),
#' 	Rep_SD=c(fit.all[[1]]$aov.tab["error","SD"], fit.all[[2]]$aov.tab["error","SD"],
#' 	         fit.all[[3]]$aov.tab["error","SD"], fit.all[[4]]$aov.tab["error","SD"],
#'           fit.all[[5]]$aov.tab["error","SD"], fit.all[[6]]$aov.tab["error","SD"]),
#' 	Rep_CV=c(fit.all[[1]]$aov.tab["error","CV[%]"],fit.all[[2]]$aov.tab["error","CV[%]"],
#'           fit.all[[3]]$aov.tab["error","CV[%]"],fit.all[[4]]$aov.tab["error","CV[%]"],
#'           fit.all[[5]]$aov.tab["error","CV[%]"],fit.all[[6]]$aov.tab["error","CV[%]"]),
#'  WLP_SD=c(sqrt(sum(fit.all[[1]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[2]]$aov.tab[3:4, "VC"])),
#'           sqrt(sum(fit.all[[3]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[4]]$aov.tab[3:4, "VC"])),
#'           sqrt(sum(fit.all[[5]]$aov.tab[3:4,"VC"])),sqrt(sum(fit.all[[6]]$aov.tab[3:4, "VC"]))),
#'  WLP_CV=c(sqrt(sum(fit.all[[1]]$aov.tab[3:4,"VC"]))/fit.all[[1]]$Mean*100,
#'           sqrt(sum(fit.all[[2]]$aov.tab[3:4,"VC"]))/fit.all[[2]]$Mean*100,
#'           sqrt(sum(fit.all[[3]]$aov.tab[3:4,"VC"]))/fit.all[[3]]$Mean*100,
#'           sqrt(sum(fit.all[[4]]$aov.tab[3:4,"VC"]))/fit.all[[4]]$Mean*100,
#'           sqrt(sum(fit.all[[5]]$aov.tab[3:4,"VC"]))/fit.all[[5]]$Mean*100,
#'           sqrt(sum(fit.all[[6]]$aov.tab[3:4,"VC"]))/fit.all[[6]]$Mean*100),
#'  Repro_SD=c(fit.all[[1]]$aov.tab["total","SD"],fit.all[[2]]$aov.tab["total","SD"],
#'             fit.all[[3]]$aov.tab["total","SD"],fit.all[[4]]$aov.tab["total","SD"],
#'             fit.all[[5]]$aov.tab["total","SD"],fit.all[[6]]$aov.tab["total","SD"]),
#'  Repro_CV=c(fit.all[[1]]$aov.tab["total","CV[%]"],fit.all[[2]]$aov.tab["total","CV[%]"],
#'             fit.all[[3]]$aov.tab["total","CV[%]"],fit.all[[4]]$aov.tab["total","CV[%]"],
#'             fit.all[[5]]$aov.tab["total","CV[%]"],fit.all[[6]]$aov.tab["total","CV[%]"]))
#'   
#'  for(i in 3:8) reproMat[,i] <- round(reproMat[,i],digits=ifelse(i%%2==0,1,3))
#'  reproMat
#' 
#' # now plot the precision profile over all samples
#' plot(reproMat[,"Mean"], reproMat[,"Rep_CV"], type="l", main="Precision Profile CA19-9",
#' 		xlab="Mean CA19-9 Value", ylab="CV[%]")
#' grid()
#' points(reproMat[,"Mean"], reproMat[,"Rep_CV"], pch=16)
#' 
#' 
#' # REML-estimation not yes optimzed to the same degree as
#' # ANOVA-estimation. Note, that no variance-covariance matrix
#' # for the REML-fit is computed (VarVC=FALSE)!
#' # Note: A correct analysis would be done per-sample, this is just
#' #       for illustration.
#' data(VCAdata1)
#' system.time(fit0 <- anovaVCA(y~sample+(device+lot)/day/run, VCAdata1))
#' system.time(fit1 <- remlVCA(y~sample+(device+lot)/day/run, VCAdata1, VarVC=FALSE))
#' 
#' # The previous example will also be interesting for environments using MKL.
#' # Run it once in a GNU-R environment and once in a MKL-environment
#' # and compare computational time of both. Note, that 'VarVC' is now set to TRUE
#' # and variable "sample" is put into the brackets increasing the number of random
#' # effects by factor 10. On my Intel Xeon E5-2687W 3.1 GHz workstation it takes
#' # ~ 400s with GNU-R and ~25s with MKL support (MRO) both run under Windows.
#' system.time(fit2 <- remlVCA(y~(sample+device+lot)/day/run, VCAdata1, VarVC=TRUE))
#' 
#' # using the SWEEP-Operator is even faster but the variance-covariance matrix of
#' # VC is not automatically approximated as for fitting via REML
#' system.time(fit3 <- anovaVCA(y~(sample+device+lot)/day/run, VCAdata1))
#' fit2
#' fit3
#' }

remlVCA <- function(form, Data, by=NULL, VarVC=TRUE, quiet=FALSE, order.data=TRUE)
{
	Call <- match.call()
	
	if(!is.null(by))
	{
		stopifnot(is.character(by))
		stopifnot(by %in% colnames(Data))
		stopifnot(is.factor(by) || is.character(by))
		
		levels  <- unique(Data[,by])
		res <- lapply(levels, function(x) remlVCA(form=form, Data[Data[,by] == x,], VarVC=VarVC, quiet=quiet))
		names(res) <- paste(by, levels, sep=".")
		return(res)
	}
	
	stopifnot(class(form) == "formula")
	stopifnot(is.logical(VarVC))
	stopifnot(is.logical(quiet))
	stopifnot(is.data.frame(Data))
	stopifnot(nrow(Data) > 2)                                               # at least 2 observations for estimating a variance
	
	if(is.null(.GlobalEnv$msgEnv))											# may removed after loading the package
		msgEnv <<- new.env(parent=emptyenv())
	
	org.form <- form
	trms <- terms(form)														# convert VCA-formula to valid lmer-formula
	stopifnot(attr(trms, "response") == 1)
	lab  <- attr(trms, "term.labels")
	Data <- orderData(Data, trms, quiet=quiet,
			exclude.numeric=FALSE, order.data=order.data)		# convert all variables to factors
	
	allObsEqual <- FALSE
	
	if(length(lab) == 0)
	{
		allObsEqual <- TRUE
		if(!quiet)
			warning("No random effects specified! Call function 'anovaVCA' instead!")
		return(anovaVCA(form, Data))
	}
	
	lab  <- paste("(1|", lab, ")", sep="")
	resp <- rownames(attr(trms, "factors"))[1]
	vars <- rownames(attr(trms, "factors"))[-1]
	
	allObsEqual <- FALSE
	
	if(all(Data[,resp] == Data[1,resp]))										# no variance detectable?
	{
		allObsEqual <- TRUE
		
		if(!quiet)
			warning("All values of response variable ", paste0("'", resp, "'"), " are equal!")
	}
	
	form <- as.formula(paste(resp, "~", paste(lab, collapse="+"), sep=""))
	
	stopifnot(resp %in% colnames(Data))
	stopifnot(is.numeric(Data[,resp]))
	
	rmInd 	<- integer()	
	resp.NA <- is.na(Data[,resp])
	
	if(any(resp.NA))
	{    
		rmInd <- c(rmInd, which(resp.NA))
		if(!quiet)
			message("There are ", length(which(resp.NA))," missing values for the response variable (obs: ", paste(which(resp.NA), collapse=", "), ")!")
	}    
	
	if(!is.null(vars))
	{
		for(i in vars)                                                  # convert all nested factors as factor-objects
		{
			if( any(is.na(Data[,i])))
			{
				NAind <- which(is.na(Data[,i]))
				rmInd <- c(rmInd, NAind)
				if(!quiet)
					message("Variable '", i,"' has ",length(NAind)," missing values (obs: ", paste(NAind, collapse=", "), ")!" )
			}
			
			if(length(levels(Data[,i])) >= 0.75*nrow(Data) && !quiet)
				message("Variable >>> ", i," <<< has at least 0.75 * nrow(Data) levels!")
		}
		rmInd <- unique(rmInd)
	}	
	
	vcol <- rownames(attr(trms, "factors"))
	if(is.null(vcol))
		vcol <- resp
	Ndata <- nrow(Data)
	Data  <- na.omit(Data[,vcol, drop=F])
	Nobs <- nrow(Data)
	
	if(quiet || allObsEqual)
		suppressWarnings(fit <- lmer(form, Data))						# fit via 'lmer'
	else
		fit <- lmer(form, Data)
	
	res <- list(call=Call, Type="Random Model", EstMethod="REML", data=Data, terms=trms,
			intercept=as.logical(attr(trms, "intercept")), response=resp)
	
	res$Mean 	 	 <- mean(Data[,resp], na.rm=TRUE)
	res$formula		 <- org.form												# as specified by the user
	res$form 	 	 <- form
	res$Nvc  	 	 <- length(lab) + 1											# error is one additional VC
	res$VCnames  	 <- c(attr(trms, "term.labels", "error"))
	res$NegVCmsg 	 <- ""														# there will never be anything to report
	res$VarVC.method <- "gb"
	res$balanced <- if(isBalanced(as.formula(trms), Data)) 
				"balanced"  
			else 
				"unbalanced" 
	
	res$terms.classes <- sapply(Data[,rownames(attr(trms, "factors"))[-1]], class)
	
	if(Nobs != Ndata)
		res$Nrm <- Ndata - Nobs                         # save number of observations that were removed due to missing data
	
	res$Nobs <- Nobs
	
	tmp <- lmerSummary(	obj=fit, VarVC=VarVC, 			# construct table similar to aov-table and approximate vcovVC
			terms=attr(trms, "term.labels"),
			Mean=res$Mean)	
	res <- c(res, tmp)
	class(res) <- "VCA"
	
	if(allObsEqual)
	{
		res$aov.tab[,"DF"] <- NA
		res$aov.tab[,"%Total"] <- 0	
	}
	res
}


#' Fit Linear Mixed Models via REML.
#' 
#' Function fits Linear Mixed Models (LMM) using Restricted Maximum Likelihood (REML).
#' 
#' The model is formulated exactly as in function \code{\link{anovaMM}}, i.e. random terms need be enclosed by round brackets.
#' All terms appearing in the model (fixed or random) need to be compliant with the regular expression "^[^[\\.]]?[[:alnum:]_\\.]*$",
#' i.e. they may not start with a dot and may then only consist of alpha-numeric characters, 
#' dot and underscore. Otherwise, an error will be issued.
#' 
#' Here, a LMM is fitted by REML using the \code{\link{lmer}} function of the \code{lme4}-package. 
#' For all models the Giesbrechnt & Burns (1985) approximation of the variance-covariance
#' matrix of variance components (VC) can be applied ('VarVC=TRUE'). A Satterthwaite approximation of the degrees of freedom
#' for all VC and total variance is based on this approximated matrix using \eqn{df=2Z^2}, where
#' \eqn{Z} is the Wald statistic \eqn{Z=\sigma^2/se(\sigma^2)}, and \eqn{\sigma^2} is here used for an
#' estimated variance. The variance of total variability, i.e. the sum of all VC is computed via summing
#' up all elements of the variance-covariance matrix of the VC.
#' One can constrain the variance-covariance matrix of random effects \eqn{G} to be either diagonal ('cov=FALSE'), i.e.
#' all random effects are indpendent of each other (covariance is 0). If 'cov=TRUE' (the default) matrix \eqn{G} will be
#' constructed as implied by the model returned by function \code{\link{lmer}}. 
#' 
#' As for objects returned by function \code{\link{anovaMM}} linear hypotheses of fixed effects or LS Means can be
#' tested with functions \code{\link{test.fixef}} and \code{\link{test.lsmeans}}. Note, that option "contain" does
#' not work for LMM fitted via REML.
#' 
#' Note, that for large datasets approximating the variance-covariance matrix of VC is computationally expensive
#' and may take very long. There is no Fisher-information matrix available for 'merMod' objects, which can
#' serve as approximation. To avoid this time-consuming step, use argument 'VarVC=FALSE' but remember,
#' that no confidence intervals for any VC will be available. If you use Microsoft's R Open, formerly known
#' as Revolution-R, which comes with Intel's Math Kernel Library (MKL), this will be automatically detected
#' and an environment-optimized version will be used, reducing the computational time considerably (see examples).
#' 
#' @param form          (formula) specifying the model to be fit, a response variable left of the '~' is mandatory, random terms
#' 						have to be enclosed in brackets (see details for definition of valid model terms)
#' @param Data          (data.frame) containing all variables referenced in 'form'
#' @param by			(factor, character) variable specifying groups for which the analysis should be performed individually,
#' 						i.e. by-processing
#' @param VarVC			(logical) TRUE = the variance-covariance matrix of variance components will be approximated using 
#' 						the method found in Giesbrecht & Burns (1985), which also serves as basis for applying a Satterthwaite
#' 						approximation of the degrees of freedom for each variance component, FALSE = leaves out this step, 
#' 						no confidence intervals for VC will be available
#' @param cov			(logical) TRUE = in case of non-zero covariances a block diagonal matrix will be constructed,
#'                  	FALSE = a diagonal matrix with all off-diagonal element being equal to zero will be contructed
#' @param quiet			(logical) TRUE = will suppress any messages or warning, which will be issued otherwise 
#' @param order.data	(logical) TRUE = class-variables will be ordered increasingly, FALSE = ordering of class-variables
#' 						will remain as is
#' 
#' @seealso \code{\link{remlVCA}}, \code{\link{VCAinference}}, \code{\link{ranef.VCA}}, \code{\link{residuals.VCA}},
#' 			\code{\link{anovaVCA}}, \code{\link{anovaMM}}, \code{\link{plotRandVar}},  \code{\link{test.fixef}},  
#' 			\code{\link{test.lsmeans}}, \code{\link{lmer}}
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples
#' \dontrun{
#' data(dataEP05A2_2)
#' 
#' # assuming 'day' as fixed, 'run' as random
#' remlMM(y~day/(run), dataEP05A2_2)
#' 
#' # assuming both as random leads to same results as
#' # calling anovaVCA
#' remlMM(y~(day)/(run), dataEP05A2_2)
#' anovaVCA(y~day/run, dataEP05A2_2)
#' remlVCA(y~day/run, dataEP05A2_2)
#' 
#' # fit a larger random model
#' data(VCAdata1)
#' fitMM1 <- remlMM(y~((lot)+(device))/(day)/(run), VCAdata1[VCAdata1$sample==1,])
#' fitMM1
#' # now use function tailored for random models
#' fitRM1 <- anovaVCA(y~(lot+device)/day/run, VCAdata1[VCAdata1$sample==1,])
#' fitRM1
#' 
#' # there are only 3 lots, take 'lot' as fixed 
#' fitMM2 <- remlMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,])
#' 
#' # the following model definition is equivalent to the one above,
#' # since a single random term in an interaction makes the interaction
#' # random (see the 3rd reference for details on this topic)
#' fitMM3 <- remlMM(y~(lot+(device))/day/run, VCAdata1[VCAdata1$sample==2,])
#' 
#' # fit same model for each sample using by-processing
#' lst <- remlMM(y~(lot+(device))/day/run, VCAdata1, by="sample")
#' lst
#' 
#' # fit mixed model originally from 'nlme' package
#'  
#' library(nlme)
#' data(Orthodont)
#' fit.lme <- lme(distance~Sex*I(age-11), random=~I(age-11)|Subject, Orthodont) 
#' 
#' # re-organize data for using 'remlMM'
#' Ortho <- Orthodont
#' Ortho$age2 <- Ortho$age - 11
#' Ortho$Subject <- factor(as.character(Ortho$Subject))
#' fit.remlMM1 <- remlMM(distance~Sex*age2+(Subject)*age2, Ortho)
#' 
#' # use simplified formula avoiding unnecessary terms
#' fit.remlMM2 <- remlMM(distance~Sex+age2+Sex:age2+(Subject)+age2:(Subject), Ortho)
#' 
#' # and exclude intercept
#' fit.remlMM3 <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho)
#' 
#' # now use exclude covariance of per-subject intercept and slope
#' # as for models fitted by function 'anovaMM'
#' fit.remlMM4 <- remlMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho, cov=FALSE)
#' 
#' # compare results
#' fit.lme
#' fit.remlMM1
#' fit.remlMM2
#' fit.remlMM3
#' fit.remlMM4
#' 
#' # are there a sex-specific differences?
#' cmat <- getL(fit.remlMM3, c("SexMale-SexFemale", "SexMale:age2-SexFemale:age2")) 
#' cmat
#' 			 
#' test.fixef(fit.remlMM3, L=cmat)
#' }

remlMM <- function(form, Data, by=NULL, VarVC=TRUE, cov=TRUE, quiet=FALSE, order.data=TRUE)
{
	Call <- match.call()
	
	if(!is.null(by))
	{
		stopifnot(is.character(by))
		stopifnot(by %in% colnames(Data))
		stopifnot(is.factor(by) || is.character(by))
		
		levels  <- unique(Data[,by])
		res <- lapply(levels, function(x) remlMM(form=form, Data[Data[,by] == x,], VarVC=VarVC, cov=cov, quiet=quiet))
		names(res) <- paste(by, levels, sep=".")
		return(res)
	}
	
	stopifnot(class(form) == "formula")
	stopifnot(is.logical(VarVC))
	stopifnot(is.logical(quiet))
	stopifnot(is.data.frame(Data))
	stopifnot(nrow(Data) > 2)                                               # at least 2 observations for estimating a variance
	
	if(is.null(.GlobalEnv$msgEnv))												# may removed after loading the package
		msgEnv <<- new.env(parent=emptyenv())
	
	org.form <- form
	trms <- terms(form, simplify=TRUE, keep.order=TRUE)						# convert VCA-formula to valid lmer-formula
	stopifnot(attr(trms, "response") == 1)
	resp <- rownames(attr(trms, "factors"))[1]
	org.form <- form
	Data <- orderData(Data, trms, quiet=quiet,
			order.data=order.data)
	
	if(any(!grepl("^[^[\\.]]?[[:alnum:]_\\.]*$", rownames(attr(trms, "factors")))))
		stop("There are terms in the model formula where regular expression '^[^[\\.]]?[[:alnum:]_\\.]*$' does not fit!")
	
	stopifnot(resp %in% colnames(Data))
	stopifnot(is.numeric(Data[,resp]))
	
	res <- list(call=Call,  EstMethod="REML", data=Data, terms=trms,
			response=resp)
	
	allObsEqual <- FALSE
	
	if(all(Data[,resp] == Data[1,resp]))										# no variance detectable?
	{
		allObsEqual <- TRUE
		if(!quiet)
			warning("All values of response variable ", paste0("'", resp, "'"), " are equal!")
	}
	
	int <- res$intercept <- attr(trms, "intercept") == 1						# has intercept	
	rf <- gregexpr("\\([[:alnum:]_\\.]*\\)", as.character(org.form)[3])			# check for random effects
	
	if(rf[[1]][1] != -1)														# identify random variables
	{
		len <- attr(rf[[1]], "match.length")
		pos <- rf[[1]]
		tmp <- NULL
		for(i in 1:length(len))
		{
			tmp <- c(tmp, substr(as.character(org.form)[3], pos[i]+1, pos[i]+len[i]-2))				# remember random factors, exclude brackets
		}
		rf <- tmp
	}
	
	Ndata <- nrow(Data)	
	rmInd <- integer()	
	resp.NA <- is.na(Data[,resp])
	
	if(any(resp.NA))								# remove missing data and warn
	{    
		rmInd <- c(rmInd, which(resp.NA))
		if(!quiet)
			message("There are ", length(which(resp.NA))," missing values for the response variable (obs: ", paste(which(resp.NA), collapse=", "), ")!")
		res$resp.NA <- rmInd
	}    
	
	fac  	<- attr(trms, "term.labels")
	if(length(fac) > 1)
	{
		rf.ind  <- which(apply(sapply(rf, function(x) regexpr(x, fac)), 1, function(x) any(x>0)))		
	}
	else
	{
		if(length(rf) > 0)
		{
			if(rf == fac)
				rf.ind <- 1
			else
				rf.ind <- numeric(0)
		}
	}
	vars    <- rownames(attr(trms, "factors"))[-1]	                        # remove response
	Nvc     <- length(fac) + 1 
	Order   <- NULL
	
	for(i in rev(vars))														# check Data for consistency
	{
		if( any(is.na(Data[,i])))
		{
			NAind <- which(is.na(Data[,i]))
			rmInd <- c(rmInd, NAind)
			if(!quiet)
				message("Variable '", i,"' has ",length(NAind)," missing values (obs: ", paste(NAind, collapse=", "), ")!" )
		}
	}
	
	if(length(rf.ind) > 0)													# at least one random term in 'form'
	{
		res$random <- fac[rf.ind]
		res$fixed  <- fac[-rf.ind]
	}
	else
	{
		if(!quiet)
			warning("No random terms in the model! Call 'anovaMM' insead!")
		
		return(anovaMM(form, Data))
	}
	
	res$terms.classes <- sapply(Data[,rownames(attr(trms, "factors"))[-1]], class)
	
	res$Type <- if(length(res$fixed) == 0)
				"Random Model"
			else
				"Mixed Model"			
	
	lab       <- attr(trms, "term.labels")	
	fixed     <- paste(res$fixed, collapse="+")	
	var.num   <- names(res$terms.classes[which(res$terms.classes == "numeric")])
	fac 	  <- attr(trms, "factors")[var.num,res$fixed,drop=FALSE]		# restrict to columns of fixed effects variables and rows with numeric variables
	fe.num 	  <- character()
	
	fac 	  <- fac[which(apply(fac, 1, any)),,drop=FALSE]
	fac 	  <- fac[,which(apply(fac, 2, any)),drop=FALSE]
	fe.num 	  <- colnames(fac)
	iacs      <- which(grepl(":", res$random))										# interaction terms in the formula
	num.rand  <- sapply(var.num, function(x) grepl(x, res$random)) 
	num.rand  <- as.matrix(num.rand)
	num.fixed <- sapply(var.num, function(x) grepl(x, res$fixed))
	random 	  <- ""
	trmMat    <- matrix(nrow=length(res$random), ncol=2, dimnames=list(NULL, c("form", "subject")))
	
	for(i in 1:length(res$random))			# process each random term
	{
		if(nchar(random) == 0)
			sep <- ""
		else
			sep <- "+"
		
		if(length(num.rand) > 0 && any(num.rand[i,]))				# random term with numeric variable found
		{
			if(i %in% iacs)
			{
				splt 	<- unlist(strsplit(res$random[i], ":"))
				tmp.num <- which(splt %in% var.num)
				numVar 	<- splt[ tmp.num]
				others  <- paste(splt[-tmp.num], collapse=":")				
				ind 	<- which(trmMat[,2] == others)				# subject terms may only occur once
				
				if(length(ind) > 0)
					trmMat[ind,2] <- NA
				
				trmMat[i,] <- c(numVar, others)	
			}
			else
				stop("Numeric variables may only occur in random interaction terms!")
		}
		else
			trmMat[i,] <- c("1", res$random[i])
	}
	trmMat 	<- na.omit(trmMat)
	random 	<- paste( apply(trmMat, 1, function(x) paste("(", x[1], "|", x[2], ")", sep="")), collapse="+")	
	form 	<- paste(  resp, "~", fixed, "+", random, sep="")
	form 	<- as.formula(form)
	vcol 	<- rownames(attr(trms, "factors"))
	
	if(is.null(vcol))
		vcol <- resp
	
	Ndata <- nrow(Data)
	Data  <- na.omit(Data[,vcol, drop=F])
	Nobs  <- nrow(Data)
	
	if(quiet || allObsEqual)
		suppressWarnings(fit <- lmer(form, Data))		# fit via 'lmer'
	else
		fit <- lmer(form, Data)
	
	
	res$Mean 		 <- mean(Data[,resp], na.rm=TRUE)
	res$formula		 <- org.form						# as.specified by the user
	res$form 	 	 <- form							# as used in the lmer-call
	res$NegVCmsg 	 <- ""
	res$VarVC.method <- "gb"
	
	if(int)
	{
		X <- matrix(1, ncol=1, nrow=nrow(Data))			# design matrix of fixed effects: include intercept --> needs a restriction
		colnames(X) <- "int"
		fe.assign <- 0									# '0' indicates intercept
	}
	else
	{
		fe.assign <- NULL
		X <- matrix(1, ncol=0, nrow=nrow(Data))	
	}
	
	fixed.form <- as.formula(paste(	resp, "~", 
					paste(	ifelse(int, "1", "-1"), 
							fixed, sep=ifelse(fixed=="", "", "+")),
					sep=""))
	
	old.opt <- options(contrasts=c("contr.SAS", "contr.poly"))
	
	suppressWarnings({
				X1 	<- model.matrix(fixed.form, data=Data)						# with SAS contrasts but without the columns where restrictions apply
				X10	<- apply(X1, 2, function(x) all(x==0))
				X2 <- model.matrix(	fixed.form, data=Data,						# full model matrix with re-setted contrasts
						contrasts.arg=lapply(Data[,sapply(Data, is.factor), drop=FALSE],
								contrasts, contrasts=FALSE))
				X20	<- apply(X2, 2, function(x) all(x==0))
				
				X2.asgn 	<- attr(X2, "assign")									# keep info		
				
				if(any(X10))
					X1 <- X1[,-which(X10),drop=FALSE]
				
				if(any(X20))
				{
					X2 <- X2[,-which(X20),drop=FALSE]
					X2.asgn 	<- X2.asgn[-which(X20)]
				}				
			})
	
	options(old.opt)													# reset contrasts option
	
	fixed.terms <- terms(fixed.form)
	fe.terms 	<- attr(fixed.terms, "term.labels")
	
	X2[,!colnames(X2) %in% colnames(X1)] <- 0							# transfer restriction from X1 to X2
	if(int)
	{
		colnames(X2)[1] <- "int"
		fe.terms <- c("int", fe.terms)
	}
	
	fe.assign <- X2.asgn
	attr(fe.assign, "terms") <- fe.terms
	
	res$fe.assign 	<- fe.assign										# mapping columns of X to fixed terms in the model formula
	res$fixed.terms <- fixed.terms
	res$balanced <- if(isBalanced(as.formula(trms), Data)) 
				"balanced"  
			else 
				"unbalanced" 
	
	if(Nobs != Ndata)
		res$Nrm <- Ndata - Nobs                         # save number of observations that were removed due to missing data
	
	res$Nobs <- Nobs
	
	tmp <- lmerSummary(	obj=fit, VarVC=VarVC, 			# construct table similar to aov-table and approximate vcovVC
			terms=res$random,
			Mean=res$Mean,
			cov=cov, X=X2)				
	
	tmp$Matrices$y <- Matrix(Data[,resp], ncol=1)
	res <- c(res, tmp)
	class(res) <- "VCA"
	
	if(allObsEqual)
	{
		res$aov.tab[,"DF"] <- NA
		res$aov.tab[,"%Total"] <- 0	
	}
	res
}


#' Solve System of Linear Equations using Inverse of Cholesky-Root.
#' 
#' Function solves a system of linear equations, respectively, inverts a matrix
#' by means of the inverse Cholesky-root.
#' 
#' This function is intended to reduce the computational time in function
#' \code{\link{solveMME}} which computes the inverse of the square variance-
#' covariance Matrix of observations. It is considerably faster than function
#' \code{\link{solve}} (see example).
#' Whenever an error occurs, which is the case for non positive definite matrices
#' 'X', function \code{\link{MPinv}} is called automatically yielding a generalized
#' inverse (Moore-Penrose inverse) of 'X'.
#' 
#' @param X			(matrix, Matrix) object to be inverted
#' @param quiet		(logical) TRUE = will suppress any warning, which will be issued otherwise 
#' 
#' @return (matrix, Matrix) corresponding to the inverse of X
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples
#' \dontrun{
#' # following complex (nonsense) model takes pretty long to fit
#' system.time(res.sw <- anovaVCA(y~(sample+lot+device)/day/run, VCAdata1))
#' # solve mixed model equations (not automatically done to be more efficient)
#' system.time(res.sw <- solveMME(res.sw))
#' # extract covariance matrix of observations V
#' V1 <- getMat(res.sw, "V")
#' V2 <- as.matrix(V1)
#' system.time(V2i <- solve(V2))
#' system.time(V1i <- VCA:::Solve(V1))
#' V1i <- as.matrix(V1i)
#' dimnames(V1i) <- NULL
#' dimnames(V2i) <- NULL
#' all.equal(V1i, V2i)
#' } 

Solve <- function(X, quiet=FALSE)
{
	stopifnot(ncol(X) == nrow(X))
	
	clsMatrix <- inherits(X, "Matrix")
	if(clsMatrix)
	{
		cls <- class(X)
		X <- as.matrix(X)
	}
	Xi <- try(chol2inv(chol(X)), silent=TRUE)
	
	if(class(Xi) == "try-error")			# use Moore-Penrose inverse instead in case of an error
	{										# using the Cholesky-decomposition approach
		if(!quiet)
			warning("\tMatrix inversion via 'chol2inv' failed!\n\tUse generalized (Moore-Penrose) inverse (MPinv)!", sep="\n")
		Xi <- MPinv(X)
	}
	
	if(clsMatrix)
	{
		Xi <- as(Xi, "dgCMatrix")			# Solve used here for inverting V-matrix
	}
	return(Xi)
}



#' Calling C-implementation of the SWEEP-Operator
#' 
#' Function calls a fast C-implementation of the SWEEP operator using the
#' transpose of the original augmented matrix \eqn{X'X} (see \code{\link{getSSQsweep}}).
#' 
#' Transposing prior to applying the SWEEP-operator speeds up things since the 
#' complete matrix is stored in memory in consecutive manner. 
#' 
#' This is an utility-function not intended to be called directly.
#' 
#' @param M			(matrix) matrix, representing the augmented matrix \eqn{X'X}
#' @param asgn		(integer) vector, identifying columns in \eqn{M} corresponding to variables, 
#' 					respectively, to their coefficients
#' @param thresh	(numeric) value used to check whether the influence of the a coefficient
#' 					to reducing the error sum of squares is small enough to conclude that the
#' 					corresponding column in \eqn{X'X} is a linear combination of preceding 
#' 					columns
#' @param tol		(numeric) value used to check numerical equivalence to zero
#' @param Ncpu		(integer) number of cores to be used for parallel processing
#'                  (not yet used)
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @return (list) with two elements:\cr
#' 			\item{SSQ}{(numeric) vector of ANOVA sum of squares}
#' 			\item{LC}{(integer) vector indicating linear dependence of each column}
#' 
#' @references 
#' Goodnight, J.H. (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, 149-158
#' 
#' @examples 
#' \dontrun{
#' # use example from 'lm' Rdoc
#' ctl    <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
#' trt    <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
#' group  <- gl(2, 10, 20, labels = c("Ctl","Trt"))
#' weight <- c(ctl, trt)
#' lm.D9  <- lm(weight ~ group)
#' anova(lm.D9)
#' 
#' # create augmented matrix
#' X  <- model.matrix(lm.D9)
#' Xt <- t(X)
#' y  <- matrix(weight, ncol=1)
#' yt <- t(y)
#' M  <- rbind(	cbind(as.matrix(Xt%*%X), as.matrix(Xt%*%y)), 
#' 		        cbind(as.matrix(yt%*%X), as.matrix(yt%*%y)))
#' swept <- VCA:::Csweep(M, asgn=c(0,1))					
#' LC    <- swept$LC
#' SS    <- swept$SSQ
#' SSQ   <- abs(diff(SS))
#' SSQ   <- c(SSQ, tail(SS,1))
#' SSQ		# compare to column "Sum Sq" in 'anova(lmD9)' output
#' }

Csweep <- function(M, asgn, thresh=1e-12, tol=1e-12, Ncpu=1)
{
	nr <- nrow(M)
	stopifnot(nr == ncol(M))
	LC <- ZeroK <- rep(0, length(asgn))
	
	ind <- is.nan(M) | is.na(M) | M == Inf | M == -Inf;
	if(any(ind))
		M[which(ind)] <- 0
	
	swept <- .C("Tsweep", M=as.double(t(M)), k=as.integer(asgn), thresh=as.double(thresh), 
				NumK=as.integer(length(asgn)), nr=as.integer(nr), LC=as.integer(LC), 
				tol=as.double(tol), SSQ=as.double(rep(0, length(unique(asgn)))), 
				PACKAGE="VCA")
	
	res <- list(SSQ=swept$SSQ, 
			LC=tapply(swept$LC, asgn, function(x) length(which(x==1))))
	
	return(res)
}


#' Calling C-implementation of the SWEEP-Operator for Matrix-Inversion
#' 
#' Function calls a fast C-implementation of the SWEEP operator using the
#' transpose of the matrix to be swept for generating a generalized inverse of 
#' a matrix for which no regular matrix invers exists.
#' 
#' Transposing prior to applying the SWEEP-operator speeds up things since the 
#' complete matrix is stored in memory in consecutive manner. 
#' This version of the SWEEP-operator is intended for matrix inversion only, thus,
#' not computing ANOVA sum of squares and number of linear dependencies (see function
#' \code{\link{Csweep}}).
#' 
#' This is an utility-function not intended to be called directly.
#' 
#' @param M			(matrix) matrix, representing the augmented matrix \eqn{X'X}
#' @param tol		(numeric) value used to check numerical equivalence to zero
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @return (Matrix) object corresponding to the inverted matrix
#' 
#' @examples
#' \dontrun{
#' M <- matrix(c(4,-6,6,-9),2)
#' solve(m)			# regular inverse does not exist
#' Mi1 <- MPinv(M)			# MASS-implementation 'ginv'
#' Mi1
#' M %*% Mi1 %*% M			# should be M
#' Mi2 <- VCA:::Sinv(M)
#' Mi2
#' M %*% Mi2 %*% M			# should be M
#' }

Sinv <- function(M, tol=.Machine$double.eps)
{
	nr <- nrow(M)
	stopifnot(nr == ncol(M))
	
#	ind <- is.nan(M) | is.na(M) | M == Inf | M == -Inf;
#	if(any(ind))
#		M[which(ind)] <- 0
#	
	swept <- .C("TsweepFull", M=as.double(t(M)), nr=as.integer(nr), 
				tol=as.double(tol*max(abs(diag(M)))), PACKAGE="VCA")
	
	if(class(M) == "matrix")
		return(matrix(round(swept$M, abs(log10(tol))), nrow=nr, ncol=nr, byrow=TRUE))
	else
		return(Matrix(round(swept$M, abs(log10(tol))), nrow=nr, ncol=nr, byrow=TRUE))
}



#' ANOVA Sum of Squares via Sweeping.
#' 
#' Compute ANOVA Type-1 sum of squares for linear models.
#' 
#' This function performs estimation of ANOVA Type-1 sum of squares
#' using the SWEEP-operator (see reference), operating on the augmented
#' matrix \eqn{X'X}, where \eqn{X} represents the design matrix not differentiating
#' between fixed and random factors. See the numerical example in \code{\link{Csweep}}
#' exemplifying the type of augmentation of \eqn{X'X} on which sweeping is carried out.
#' 
#' This is an utility function not intended to be called directly.
#' For each term in the formula the design-matrix \eqn{Z} is constructed.
#' Matrix \eqn{X} corresponds to binding all these \eqn{Z}-matrices together column-wise.
#' 
#' Degrees of freedom for each term are determined by subtracting the number of
#' linearly dependent columns from the total number of column in X asigned to a
#' specific term.
#' 
#' @param Data			(data.frame) with the data
#' @param tobj			(terms) object derived from original formula object
#' @param random		(character) vector, optionally containing information about each
#' 						model term, whether it is random or fixed (only used in mixed models)
#' 
#' @return (list) representing the  with variables:\cr
#' 			\item{aov.tab}{basic ANOVA-table with degrees of freedom (DF), SS and MS}
#' 			\item{Lmat}{(list) with components 'Z' and 'A'}
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @seealso \code{\link{Csweep}}
#' 
#' @references 
#' 
#' Goodnight, J.H., (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, p.149-158
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_1)
#' res <- VCA:::getSSQsweep(dataEP05A2_1, terms(y~day/run))
#' str(res)
#' }
getSSQsweep <- function(Data, tobj, random=NULL)
{	
	form     <- formula(tobj)
	resp     <- as.character(form)[2]
	fac      <- attr(tobj, "term.labels")
	int      <- attr(tobj, "intercept") == 1 
	
	N        <- nrow(Data)															
	SS       <- numeric()
	DF       <- numeric()
	Lmat     <- list()                                                      	# compute A-matrices, used in constructing covariance-matrix of VCs
	Lmat$Z   <- list()
	Lmat$Zre <- Matrix(nrow=N, ncol=0)
	y        <- Matrix(Data[,resp], ncol=1)
	
	ord <- attr(tobj, "order")
	vars <- rownames(attr(tobj, "factors"))
	fac <- c(fac, "error")
	DF  <- rep(0, length(fac))
	done <- rep(FALSE, length(fac))					# keep track which DFs have been computed
	names(DF) <- fac
	vars <- vars[-1]
	Nvc  <- length(fac) 							# error not included
	N <- nrow(Data)
	asgn <- NULL
	Lmat <- list()
	Lmat$Z <- list()
	
	if(int)
	{
		Lmat$Zre <- Matrix(1, nrow=N, ncol=1)
		colnames(Lmat$Zre) <- "int"
		asgn     <- 0
	}	
	else
		Lmat$Zre <- Matrix(0, nrow=N, ncol=0)
	
	NumK <- c(rep(0, Nvc-1), N)
	
	for(i in 1:Nvc)                 				# construct design-matrices Z for each term, and matrix X (here Zre)                       
	{
		if(i < Nvc)
		{
			tmpMM <- model.matrix(as.formula(paste(resp, "~", fac[i], "-1", sep="")), Data)
			Lmat$Z[[i]] <- Matrix(tmpMM)
			all0 <- apply(Lmat$Z[[i]], 2, function(x) all(x==0))
			if(any(all0))
			{
				ZeroCol <- which(all0)
				Lmat$Z[[i]] <- Lmat$Z[[i]][,-ZeroCol]
			}	
			if(!is.null(random))
				attr(Lmat$Z[[i]], "type") <- ifelse(fac[i] %in% random, "random", "fixed")
			
			attr(Lmat$Z[[i]], "term") <- fac[i]
			
			NumK[i] <- ncol(Lmat$Z[[i]])
			Lmat$Zre    <- Matrix(cbind(as.matrix(Lmat$Zre), as.matrix(Lmat$Z[[i]])))
			asgn <- c(asgn, rep(i, ncol(Lmat$Z[[i]])))
		}
		else
		{
			Lmat$Z[[Nvc]] <- Diagonal(N)							# include error design matrix
			attr(Lmat$Z[[Nvc]], "type") <- "random"
			attr(Lmat$Z[[Nvc]], "term") <- "error"
		}
	}
	
	attr(Lmat$Zre, "assign") <- asgn
	X <- Lmat$Zre
	Xt <- t(X)
	y  <- Matrix(Data[,resp], ncol=1)
	yt <- t(y)
	
	M <- rbind(	cbind(as.matrix(Xt%*%X), as.matrix(Xt%*%y)), 
				cbind(as.matrix(yt%*%X), as.matrix(yt%*%y)))	
	
	uind <- unique(asgn)							# all factors
	SS <- LC <- NULL
	nr <- nrow(M)
	tmpM <- M
	
	Int <- ifelse(int, 1, 0)
	
	V <- rep(1, ncol(M))
	swept <- Csweep(M, asgn=asgn)					# sweep matrix M
	LC    <- swept$LC
	SS	  <- swept$SSQ
	
	if(int)											# no DF-adjustment for intercept
		LC <- LC[-1]
	
	if(Nvc > 1)
	{
		DF[1:(Nvc-1)] <- NumK[1:(Nvc-1)]-LC				# adjust for linearly dependent variables --> resulting in degrees of freedom
		DF["error"] <- tail(NumK,1)-sum(DF[1:(length(DF)-1)]) - Int
	}
	else
		DF["error"] <- NumK - Int
	
	if(!int)
		SS <- c(M[nr, nr], SS)
	
	SSQ <- abs(diff(SS))
	SSQ <- c(SSQ, tail(SS,1))
	
	aov.tab <- data.frame(DF=DF, SS=SSQ, MS=SSQ/DF)	# basic ANOVA-table
	
	return(list(aov.tab=aov.tab, Lmat=Lmat))
}


#' Get Orthogonal Basis of a Matrix.
#' 
#' Calculate an orthogonal Basis for the span of the column vectors of the input matrix.
#' 
#' @param X			(object) two-dimensional, for which a Basis inverse
#'                  has to be computed
#' @param tol		(numeric) tolerance value to be used in comparisons
#' 
#' @return (object) A matrix whose columns form an orthogonal basis for the span of X
#' 
#' @author Florian Dufey \email{[email protected]@contractors.roche.com}

getBasis <- function (X, tol = sqrt(.Machine$double.eps)) 
{
	if (length(dim(X)) > 2L) 
		stop("'X' must be two-dimensional")
	Xsvd 		<- svd(X,,nv=0)
	Positive 	<- Xsvd$d > max(tol * Xsvd$d[1L], 0)
	if (all(Positive)) 
		Xsvd$u
	else if (!any(Positive)) 
		array(0, dim(X)[2L:1L])
	else 
		Xsvd$u[, Positive, drop = FALSE] 
}


#' ANOVA Sum of Squares via Quadratic Forms
#' 
#' Compute ANOVA Type-1 sum of squares for linear models.
#' 
#' This function performs estimation of ANOVA Type-1 sum of squares
#' using an approach of expressing them as quadratic forms in \code{y},
#' the column vector of observations. This is an utility function not
#' intended to be called directly.
#' For each term in the formula the design-matrix \code{Z} and the corresponding
#' \code{A}-matrix 
#' Degrees of freedom for each term are determined calling function \code{\link{anovaDF}}.
#' 
#' @param Data			(data.frame) with the data
#' @param tobj			(terms) object derived from original formula object
#' @param random		(character) vector, optionally containing information about each
#' 						model term, whether it is random or fixed (only used in mixed models)
#' 
#' @return (list) representing the  with variables:\cr
#' 			\item{aov.tab}{basic ANOVA-table with degrees of freedom (DF), SS and MS}
#' 			\item{Lmat}{(list) with components 'Z' and 'A'}
#' 
#' @author 	Andre Schuetzenmeister \email{[email protected]@roche.com},
#' 			Florian Dufey \email{[email protected]@contractors.roche.com}
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_1)
#' res <- VCA:::getSSQqf(dataEP05A2_1, terms(y~day/run))
#' str(res)
#' }

getSSQqf<- function(Data, tobj, random=NULL)
{
	form     <- formula(tobj)
	resp     <- as.character(form)[2]
	fac      <- attr(tobj, "term.labels")
	int      <- attr(tobj, "intercept") == 1 
	
	N        <- nrow(Data)															
	SS       <- numeric()
	DF       <- c()
	Lmat     <- list()                # compute Baiss-matrices, used in constructing covariance-matrix of VCs
	Lmat$Z   <- list()
	asgn     <- NULL
	
	if(int)
	{
		Lmat$Zre <- Matrix(1, nrow=N, ncol=1)
		colnames(Lmat$Zre) <- "int"
		asgn     <- 0
	}
	else
		Lmat$Zre <- Matrix(0, nrow=N, ncol=0)
	
	y        <- Matrix(Data[,resp], ncol=1)
	Lmat$B   <- list()
	if (int)
		Bre <- Matrix(1/sqrt(N),ncol=1,nrow=N) #cumulative vector of Basis-matrices
	else
		Bre <- Matrix(,nrow=N, ncol=0)
	
	sumDF  	<- ifelse(int,1,0)
	Nvc 	<- length(fac) + 1
	
	X <- NULL							# X-matrix for building B-matrices 
	
	for(i in 1:Nvc)                     # over variance components (including error)                        
	{
		if(i < Nvc)
		{
			Lmat$Z[[i]] <- Matrix(model.matrix(as.formula(paste(resp, "~", fac[i], "-1", sep="")), Data))
			
			all0 <- apply(Lmat$Z[[i]], 2, function(x) all(x==0))
			if(any(all0))
				Lmat$Z[[i]] <- Lmat$Z[[i]][,-which(all0)]
			
			asgn 		<- c(asgn, rep(i, ncol(Lmat$Z[[i]])))
			Lmat$B[[i]] <- as.matrix(getBasis(X=as.matrix(Lmat$Z[[i]]-Bre %*% t(Bre) %*% Lmat$Z[[i]])))
			DFtemp    	<- ncol(Lmat$B[[i]])
			sumDF     	<-sumDF+DFtemp
			DF        	<- c(DF,DFtemp)
			Bre       	<- cbind(Bre, Lmat$B[[i]])
			
			if(!is.null(random))
			{
				attr(Lmat$Z[[i]], "type") <- ifelse(fac[i] %in% random, "random", "fixed")
				attr(Lmat$B[[i]], "type") <- ifelse(fac[i] %in% random, "random", "fixed")
			}			
			attr(Lmat$Z[[i]], "term") <- fac[i]
			attr(Lmat$B[[i]], "term") <- fac[i]
			
			Lmat$Zre    <- Matrix(cbind(as.matrix(Lmat$Zre), as.matrix(Lmat$Z[[i]])))
			By 			<- t(y) %*% Lmat$B[[i]]
			SS 			<- c(SS, as.numeric(By %*% t(By)))				# calculate ANOVA sum of squares
		}
		else
		{
			Lmat$Z[[i]] <- Diagonal(N)                                		# error
			Lmat$B[[i]] <- matrix(1,nrow=N,ncol=N-sumDF) ## dummy B as only DF of this matrix is calculated in getVCvar
			DF        	<- c(DF,N-sumDF)
			attr(Lmat$Z[[i]], "type") <- "random"
			attr(Lmat$B[[i]], "type") <- "random"
			attr(Lmat$Z[[i]], "term") <- "error"
			attr(Lmat$B[[i]], "term") <- "error"
			By <- t(y)-t(y) %*% Bre %*% t(Bre)
			SS <- c(SS, as.numeric(By %*% t(By)))				# calculate ANOVA sum of squares
		}
	}
	
	Nlvl <- sapply(Lmat$Z, ncol)									# number of unique factor levels
	names(Nlvl) <- c(fac,"error")
	names(DF) <- c(fac, "error")
	attr(DF, "Nlevel") <- Nlvl
	attr(Lmat$Zre, "assign") <- asgn
	aov.tab <- data.frame(DF=DF, SS=SS, MS=SS/DF)
	rownames(aov.tab) <- c(attr(tobj, "term.labels"), "error")
	return(list(aov.tab=aov.tab, Lmat=Lmat))
}


#' ANOVA-Type Estimation of Mixed Models.
#' 
#' Estimate/Predict random effects employing ANOVA-type estimation and obtain generalized least squares estimates
#' of fixed effects for any linear mixed model including random models and linear models.
#' 
#' A Linear Mixed Model, noted in standard matrix notation, can be written as {y = Xb + Zg + e}, where
#' \eqn{y} is the column vector of observations, \eqn{X} and \eqn{Z}{Z} are design matrices assigning fixed (\eqn{b}),
#' respectively, random (\eqn{g}) effects to observations, and \eqn{e} is the column vector of residual errors.
#' Whenever there is an intercept in the model, i.e. the substring "-1" is not part of the model formula, the same
#' restriction as in SAS PROC MIXED is introduced setting the last fixed effect equal to zero. Note, that the results
#' of an linear contrasts are not affected by using an intercept or not, except that constrained fixed effects cannot
#' be part of such contrasts (one could use the intercept estimated instead).
#' 
#' Here, no further restrictions on the type of model are made. One can fit mixed models as well as random models, which
#' constitute a sub-set of mixed models (intercept being the only fixed effect). Variables must be either of type "numeric"
#' or "factor". "character" variables are automatically converted to factors and the response variable has to be numeric, of course. 
#' In case that 'class(Data[,i])' is neither one of these three options, an error is issued. 
#' Even simple linear models can be fitted, i.e. models without a random part (without \eqn{Zg}{Zg}) besides the
#' residual errors. In this case, an Analysis of Variance (ANOVA) table is computed in the same way as done by function 'anova.lm'.
#' 
#' One drawback of using ANOVA-type estimation of random effects is, that random effects are independent, i.e they have
#' zero covariance by definition \eqn{cov(g_i,g_j) = 0}. Another one is that estimated variance components may become negative
#' under certain conditions. The latter situation is addressed by setting negative variance estimates equal to zero and adapting
#' ANOVA mean squares (MS) as \eqn{MS = C * VC}, where \eqn{C} is a coefficient matrix and a function of the design matrix \eqn{[X Z]}
#' and \eqn{VC} is the column-vector of adapted variance components. The Satterthwaite approximation of total degrees of freedom 
#' (DF for total variance) will use adapted \eqn{MS}-values. 
#' 
#' Note, that setting negative VCs equal to zero results in a conservative estimate of the total variance, i.e. it will be larger than
#' the estimate including the negative VC(s). Use parameter 'NegVC=TRUE' to explicitly allow negative variance estimates. 
#' 
#' For further details on ANOVA Type-I estimation methods see \code{\link{anovaVCA}}.
#'
#' @param form				(formula) specifying the linear mixed model (fixed and random part of the model),
#' 							all random terms need to be enclosed by round brackets. Any variable not being bracketed
#'                      	will be considered as fixed. Interaction terms containing at least one random factor
#'                      	will automatically be random (Piepho et al. 2003). All terms appearing in the model 
#' 							(fixed or random) need to be compliant with the regular expression "^[^[\\.]]?[[:alnum:]_\\.]*$",
#' 							i.e. they may not start with a dot and may then only consist of alpha-numeric characters, 
#'							dot and underscore. Otherwise, an error will be issued.
#' @param Data				(data.frame) containing all variables referenced in 'form', note that variables can only be
#'                          of type "numeric", "factor" or "character". The latter will be automatically converted to "factor".
#' @param by				(factor, character) variable specifying groups for which the analysis should be performed individually,
#' 							i.e. by-processing
#' @param VarVC.method		(character) string specifying whether to use the algorithm given in Searle et al. (1992) which corresponds to \code{VarVC.method="scm"} or in
#' 							Giesbrecht and Burns (1985) which can be specified via "gb". Method "scm" (Searle, Casella, McCulloch)
#'                      	is the exact algorithm but slower, "gb" (Giesbrecht, Burns) is termed "rough approximation"
#' 							by the authors, but sufficiently exact compared to e.g. SAS PROC MIXED (method=type1) which
#' 							uses the inverse of the Fisher-Information matrix as approximation. For balanced designs all
#'                      	methods give identical results, only in unbalanced designs differences occur. 
#' @param SSQ.method		(character) string specifying the method used for computing ANOVA Type-1 sum of squares and respective degrees of freedom.
#' 							In case of "sweep" funtion \code{\link{getSSQsweep}} will be called, otherwise, function \code{\link{getSSQqf}}
#' @param NegVC         	(logical) FALSE = negative variance component estimates (VC) will be set to 0 and they will not 
#' 							contribute to the total variance (as done e.g. in SAS PROC NESTED, conservative estimate of total variance). 
#' 							The original ANOVA estimates can be found in element 'VCoriginal'. 
#'                      	The degrees of freedom of the total variance are based on adapted mean squares (MS) (see details).
#' 							TRUE = negative variance component estimates will not be set to 0 and they will contribute to the total 
#' 							variance (original definition of the total variance).
#' @param quiet				(logical) TRUE = will suppress any warning, which will be issued otherwise 
#' @param order.data		(logical) TRUE = class-variables will be ordered increasingly, FALSE = ordering of class-variables
#' 							will remain as is
#' @return (VCA) object
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @references	
#' 
#' Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York	
#' 
#' Goodnight, J.H. (1979), A Tutorial on the SWEEP Operator, The American Statistician, 33:3, 149-158
#' 
#' Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample
#' Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486
#' 
#' H.P.Piepho, A.Buechse and K.Emrich (2003), A Hitchhiker's Guide to Mixed Models for Randomized Experiments,
#' J.Agronomy & Crop Science 189, p. 310-322
#' 
#' Gaylor,D.W., Lucas,H.L., Anderson,R.L. (1970), Calculation of Expected Mean Squares by the Abbreviated Doolittle and Square Root Methods., 
#' Biometrics 26 (4): 641-655 
#' 
#' SAS Help and Documentation PROC MIXED, SAS Institute Inc., Cary, NC, USA
#' 
#' @seealso \code{\link{anovaVCA}}
#' 
#' @examples 
#' \dontrun{
#' 
#' data(dataEP05A2_2)
#' 
#' # assuming 'day' as fixed, 'run' as random
#' anovaMM(y~day/(run), dataEP05A2_2)
#' 
#' # assuming both as random leads to same results as
#' # calling anovaVCA
#' anovaMM(y~(day)/(run), dataEP05A2_2)
#' anovaVCA(y~day/run, dataEP05A2_2)
#' 
#' # use different approaches to estimating the covariance of 
#' # variance components (covariance parameters)
#' dat.ub <- dataEP05A2_2[-c(11,12,23,32,40,41,42),]			# get unbalanced data
#' m1.ub <- anovaMM(y~day/(run), dat.ub, SSQ.method="qf", VarVC.method="scm")
#' m2.ub <- anovaMM(y~day/(run), dat.ub, SSQ.method="qf", VarVC.method="gb")		# is faster
#' V1.ub <- round(vcovVC(m1.ub), 12)
#' V2.ub <- round(vcovVC(m2.ub), 12)
#' all(V1.ub == V2.ub)
#' 
#' # make it explicit that "gb" is faster than "scm"
#' # compute variance-covariance matrix of VCs 10-times
#' 
#' system.time(for(i in 1:500) vcovVC(m1.ub))	# "scm"
#' system.time(for(i in 1:500) vcovVC(m2.ub))	# "gb"
#' 
#' 
#' # fit a larger random model
#' data(VCAdata1)
#' fitMM1 <- anovaMM(y~((lot)+(device))/(day)/(run), VCAdata1[VCAdata1$sample==1,])
#' fitMM1
#' # now use function tailored for random models
#' fitRM1 <- anovaVCA(y~(lot+device)/day/run, VCAdata1[VCAdata1$sample==1,])
#' fitRM1
#' 
#' # there are only 3 lots, take 'lot' as fixed 
#' fitMM2 <- anovaMM(y~(lot+(device))/(day)/(run), VCAdata1[VCAdata1$sample==2,])
#' 
#' # the following model definition is equivalent to the one above,
#' # since a single random term in an interaction makes the interaction
#' # random (see the 3rd reference for details on this topic)
#' fitMM3 <- anovaMM(y~(lot+(device))/day/run, VCAdata1[VCAdata1$sample==2,])
#' 
#' # fit same model for each sample using by-processing
#' lst <- anovaMM(y~(lot+(device))/day/run, VCAdata1, by="sample")
#' lst
#' 
#' # fit mixed model originally from 'nlme' package
#'  
#' library(nlme)
#' data(Orthodont)
#' fit.lme <- lme(distance~Sex*I(age-11), random=~I(age-11)|Subject, Orthodont) 
#' 
#' # re-organize data for using 'anovaMM'
#' Ortho <- Orthodont
#' Ortho$age2 <- Ortho$age - 11
#' Ortho$Subject <- factor(as.character(Ortho$Subject))
#' fit.anovaMM1 <- anovaMM(distance~Sex*age2+(Subject)*age2, Ortho)
#' 
#' # use simplified formula avoiding unnecessary terms
#' fit.anovaMM2 <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2, Ortho)
#' 
#' # and exclude intercept
#' fit.anovaMM3 <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho)
#' 
#' # compare results
#' fit.lme
#' fit.anovaMM1
#' fit.anovaMM2
#' fit.anovaMM3
#' 
#' # are there a sex-specific differences?
#' cmat <- getL(fit.anovaMM3, c("SexMale-SexFemale", "SexMale:age2-SexFemale:age2")) 
#' cmat
#' 			 
#' test.fixef(fit.anovaMM3, L=cmat)
#' 
#' # former versions of the package used R-function 'lm' and 'anova',
#' # which is significantly slower for sufficiently large/complex models
#' data(realData)
#' datP1 <- realData[realData$PID==1,]
#' system.time(anova.lm.Tab <- anova(lm(y~lot/calibration/day/run, datP1)))
#' # using the quadratic forms approach for estimating ANOVA Type-1 sums of squares
#' system.time(anovaMM.Tab1  <- anovaMM(y~lot/calibration/day/run, datP1, SSQ.method="qf"))
#' # using the sweeping approach for estimating ANOVA Type-1 sums of squares
#' # this is now the default setting (Note: only "gb" method works as VarVC.method)
#' # Also see ?anovaVCA for a comparison of the computational efficiency of "qf" and "sweep".
#' system.time(anovaMM.Tab2  <- anovaMM(y~lot/calibration/day/run, datP1, SSQ.method="sweep"))
#' 
#' # compare results, note that the latter corresponds to a linear model,
#' # i.e. without random effects. Various matrices have already been computed,
#' # e.g. "R", "V" (which are identical in this case).
#' anova.lm.Tab
#' anovaMM.Tab1
#' anovaMM.Tab2
#' }
#' 
#' @seealso \code{\link{anovaVCA}}, \code{\link{VCAinference}}, \code{\link{remlVCA}}, \code{\link{remlMM}}
#' 			\code{\link{ranef}}, \code{\link{fixef}}, \code{\link{vcov}}, \code{\link{vcovVC}}, 
#' 			\code{\link{test.fixef}}, \code{\link{test.lsmeans}}, \code{\link{plotRandVar}}

anovaMM <- function(form, Data, by=NULL, VarVC.method=c("gb", "scm"), SSQ.method=c("sweep", "qf"), 
					NegVC=FALSE, quiet=FALSE, order.data=TRUE)
{
	if(!is.null(by))
	{
		stopifnot(is.character(by))
		stopifnot(by %in% colnames(Data))
		stopifnot(is.factor(by) || is.character(by))
		
		levels  <- unique(Data[,by])
		res <- lapply(levels, function(x) anovaMM(form=form, Data[Data[,by] == x,], NegVC=NegVC, SSQ.method=SSQ.method, VarVC.method=VarVC.method, quiet=quiet))
		names(res) <- paste(by, levels, sep=".")
		return(res)
	}
	
	stopifnot(class(form) == "formula")
	stopifnot(is.data.frame(Data))
	stopifnot(nrow(Data) > 2)                                               	# at least 2 observations for estimating a variance
	
	if(is.null(.GlobalEnv$msgEnv))												# may removed after loading the package
		msgEnv <<- new.env(parent=emptyenv())
	
	VarVC.method <- match.arg(VarVC.method)
	SSQ.method   <- match.arg(SSQ.method)
	VarVC.method <- ifelse(SSQ.method == "sweep", "gb", VarVC.method)			# always use "gb", since A-matrices will not be computed	
	
	res <- list()
	res$call <- match.call()
	if(is.null(rownames(Data)))													# assign rownames if missing to record results of re-ordering taking place before the analysis
		rownames(Data) <- 1:nrow(Data)
	tobj <- terms(form, simplify=TRUE, keep.order=TRUE)                    		# expand nested factors if necessary retain order of terms in the formula
	Data <- orderData(Data, tobj, quiet=quiet, order.data=order.data)
	res$data <- Data															# as provided 
	org.form <- form
	
	if(length(attr(tobj, "term.labels")) == 0)									# handle pure-error models with 'anovaVCA'
		return(anovaVCA(form, Data))
	
	int   <- res$intercept <- attr(tobj, "intercept") == 1						# has intercept
	form  <- formula(tobj)
	res$terms <- tobj
	
	if(any(!grepl("^[^[\\.]]?[[:alnum:]_\\.]*$", rownames(attr(tobj, "factors")))))
		stop("There are terms in the model formula where regular expression '^[^[\\.]]?[[:alnum:]_\\.]*$' does not fit!")
	
	if(!attr(tobj, "response"))
		stop("You need to include a response variable in the fixed effects formula!")
	resp <- as.character(form)[2]
	res$response <- resp
	
	stopifnot(resp %in% colnames(Data))
	stopifnot(is.numeric(Data[,resp]))
	
	rf <- gregexpr("\\([[:alnum:]_\\.]*\\)", as.character(org.form)[3])			# check for random effects
	
	if(rf[[1]][1] != -1)														# identify random variables
	{
		len <- attr(rf[[1]], "match.length")
		pos <- rf[[1]]
		tmp <- NULL
		for(i in 1:length(len))
		{
			tmp <- c(tmp, substr(as.character(org.form)[3], pos[i]+1, pos[i]+len[i]-2))				# remember random factors, exclude brackets
		}
		rf <- tmp
	}
	
	Ndata <- nrow(Data)	
	rmInd <- integer()	
	resp.NA <- is.na(Data[,resp])
	
	if(any(resp.NA))
	{    
		rmInd <- c(rmInd, which(resp.NA))
		if(!quiet)
			message("There are ", length(which(resp.NA))," missing values for the response variable (obs: ", paste(which(resp.NA), collapse=", "), ")!")
		res$resp.NA <- rmInd
	}    
	
	fac  	<- attr(tobj, "term.labels")
	
	if(length(fac) > 1)
	{
		rf.ind  <- which(apply(sapply(rf, function(x) regexpr(x, fac)), 1, function(x) any(x>0)))		
	}
	else
	{
		if(length(rf) > 0)
		{
			if(rf == fac)
				rf.ind <- 1
			else
				rf.ind <- numeric(0)
		}
	}
	vars    <- rownames(attr(tobj, "factors"))[-1]	                        # remove response
	Nvc     <- length(fac) + 1 
	
	if(length(rf.ind) > 0)													# at least one random term in 'form'
	{
		res$random <- fac[rf.ind]
		res$fixed  <- fac[-rf.ind]
	}
	else
	{
		res$random <- character(0)											# only fixed effects
		res$fixed  <- fac
	}
	
	res$VCnames <- c(res$random, "error")
	res$Nvc  	<- length(res$VCnames)											# error is one additional VC
	res$Type 	<- if(length(res$fixed) == 0)
				"Random Model"
			else
				"Mixed Model"	
	
	for(i in rev(vars))														# check Data for consistency
	{
		if( any(is.na(Data[,i])))
		{
			NAind <- which(is.na(Data[,i]))
			rmInd <- c(rmInd, NAind)
			if(!quiet)
				message("Variable '", i,"' has ",length(NAind)," missing values (obs: ", paste(NAind, collapse=", "), ")!" )
		}
	}
	attr(res$data, "analysis.order") <- rownames(Data)						# actually record the ordering of the data 
	
	rmInd <- unique(rmInd)
	
	if(length(rmInd) > 0)
		Data <- Data[-rmInd,]												
	
	Data <- na.omit(Data[,rownames(attr(tobj, "factors"))])					# get rid of incomplete observations
	Mean <- mean(Data[,resp], na.rm=TRUE)                                   # mean computed after removing incomplete observations
	Nobs <- N <- nrow(Data)
	y    <- matrix(Data[,resp], ncol=1)										# vector of observations
	
	allObsEqual <- FALSE
	if(all(Data[,resp] == Data[1,resp]))									# no variance detectable?
	{
		allObsEqual <- TRUE
		Data.org <- Data
		Data[,resp] <- Data[,resp] + rnorm(nrow(Data))
		
		if(!quiet)
			warning("All values of response variable ", paste0("'", resp, "'"), " are equal!")
	}
	
	if(SSQ.method == "qf")
		tmp.res <- getSSQqf(Data, tobj, res$random)
	else
		tmp.res <- getSSQsweep(Data, tobj, res$random)
	
	Lmat    <- tmp.res$Lmat
	aov.tab <- tmp.res$aov.tab
	DF		<- aov.tab[,"DF"]
	SS		<- aov.tab[,"SS"]
	
	if(allObsEqual)
	{
		aov.tab[,c("SS", "MS")] <- 0
		Data <- Data.org
		SS[1:length(SS)] <- 0
	}
	
	rownames(aov.tab) <- c(attr(tobj, "term.labels"), "error")
	
	res$Mean <- Mean
	res$formula <- org.form
	res$Nobs <- Nobs
	res$aov.org <- aov.tab
	
	rf.ind  <- c(rf.ind, nrow(aov.tab))
	
	C <- getCmatrix(form, Data, aov.tab[,"DF"], "SS", MM=Lmat$Zre)          # compute coefficient matrix C in ss = C * s
	
	# at this point Zre comprises fixed and random effects
	Ci  <- solve(C)
	C2  <- apply(C, 2, function(x) x/DF)                                    # coefficient matrix for mean squares (MS)
	Ci2 <- solve(C2)
	VC  <- VCorg <- as.matrix( Ci %*% SS)                                   # solve for VCs (p.173)
	VCnam <- rownames(aov.tab)
	VCnam[length(VCnam)] <- "error"
	rownames(VC) <- rownames(VCorg) <- VCnam
	colnames(VC) <- colnames(VCorg) <- "Estimate"
	
	VC  <- VC[rf.ind] 														# only use VC for random terms, inter-fixed effects variation not in scope
	aov.tab <- aov.tab[rf.ind,,drop=F]
	
	rownames(aov.tab)[nrow(aov.tab)] <- "error"
	aov.tab$VC <- VC
	
	res$NegVCmsg <- ""
	res$VCoriginal <- aov.tab[, "VC"]
	
	if(!NegVC)
	{
		IndNegVC <- which(aov.tab[,"VC"] < 0)  
		if(length(IndNegVC) > 0)                                            # there are negative VC
		{
			aov.tab[IndNegVC, "VC"] <- 0                                    # set negative VC to 0
			res$NegVCmsg <- "* VC set to 0"
		}
	} 
	
	totVC  <- sum(aov.tab$VC)
	
	aov.tab <- rbind(total=c(NA, NA, NA, totVC), aov.tab) 	
	aov.tab["total", "DF"] <- SattDF(c(C2[rf.ind, rf.ind] %*% aov.tab[-1, "VC"]), 	# will automatically adapt ANOVA-MS if any VCs were set to 0 
			Ci=Ci2[rf.ind, rf.ind, drop=F], DF=DF[rf.ind])  
	
	suppressWarnings(aov.tab <- cbind(aov.tab, SD=sqrt(aov.tab[,"VC"])))    		# warnings suppressed because sqrt of negative numbers doese not exists
	aov.tab <- cbind(aov.tab, "CV[%]"=aov.tab[,"SD"]*100/Mean)
	aov.tab <- cbind(aov.tab, "%Total"=aov.tab[,"VC"]*100/totVC)
	aov.tab <- aov.tab[,c("DF", "SS", "MS", "VC", "%Total", "SD", "CV[%]")]    
	aov.tab <- as.matrix(aov.tab)
	aov.tab <- apply(aov.tab, 1:2, function(x) ifelse(is.nan(x), NA, x))
	
	res$EstMethod <- "ANOVA"
	
	if(Nobs != Ndata)
		res$Nrm <- Ndata - Nobs                            				# save number of observations that were removed due to missing data
	
	Lmat$y <- y															# column vector of observations
	if(int)
	{
		Lmat$X    <- matrix(1, ncol=1, nrow=nrow(Data))					# design matrix of fixed effects: include intercept --> needs a restriction
		colnames(Lmat$X) <- "int"
		fe.assign <- 0													# '0' indicates intercept
	}
	else
		fe.assign <- NULL
	
	fe <- fac[-rf.ind]
	
	INT <- ifelse(int, "1", "-1")
	if(length(fe) > 0)
		INT <- paste(INT, "+", sep="")
	
	fixed.form <- paste(resp, "~", INT, paste(fe, collapse="+"), sep="")
	fixed.form <- as.formula(fixed.form)
	
	old.opt <- options(contrasts=c("contr.SAS", "contr.poly"))
	
	suppressWarnings({
				X1 	<- model.matrix(fixed.form, data=Data)						# with SAS contrasts but without the columns where restrictions apply
				X10	<- apply(X1, 2, function(x) all(x==0))
				X2 <- model.matrix(	fixed.form, data=Data,						# full model matrix with re-setted contrasts
						contrasts.arg=lapply(Data[,sapply(Data, is.factor), drop=FALSE],
								contrasts, contrasts=FALSE))
				X20	<- apply(X2, 2, function(x) all(x==0))
				
				X2.asgn 	<- attr(X2, "assign")									# keep info		
				
				if(any(X10))
					X1 <- X1[,-which(X10),drop=FALSE]
				
				if(any(X20))
				{
					X2 <- X2[,-which(X20),drop=FALSE]
					X2.asgn 	<- X2.asgn[-which(X20)]
				}				
			})
	
	options(old.opt)													# reset contrasts option
	
	fixed.terms <- terms(fixed.form)
	fe.terms 	<- attr(fixed.terms, "term.labels")
	
	X2[,!colnames(X2) %in% colnames(X1)] <- 0							# transfer restriction from X1 to X2
	if(int)
	{
		colnames(X2)[1] <- "int"
		fe.terms <- c("int", fe.terms)
	}
	
	fe.assign <- X2.asgn
	attr(fe.assign, "terms") <- fe.terms
	
	Lmat$X <- X2
	
	res$fe.assign 	<- fe.assign										# mapping columns of X to fixed terms in the model formula
	res$fixed.terms <- fixed.terms
	
	Lmat$rf.ind <- rf.ind												# indices of random effects
	Lmat$VCall  <- VCorg												# VC-estimates as if all factors were random
	Lmat$C.SS   <- C
	Lmat$C.MS   <- C2 
	Lmat$Ci.SS  <- Ci
	Lmat$Ci.MS  <- Ci2
	
	res$SSQ.method   <- SSQ.method
	res$VarVC.method <- VarVC.method
	res$aov.tab  <- aov.tab
	res$Matrices <- Lmat
	res$balanced <- if(isBalanced(form, Data)) 
				"balanced"  
			else 
				"unbalanced"
	
	class(res)     <- "VCA"
	
	if(length(Lmat$rf.ind) == 1)										# contains only the residual error
	{
		res$Type <- "Linear Model"
		tmp      <- res$aov.org
		tmp 	 <- cbind(tmp, "F value" = tmp[,"MS"]/tmp[nrow(tmp),"MS"])
		tmp		 <- cbind(tmp, "Pr(>F)" = pf(tmp[,"F value"]/tmp[nrow(tmp), "F value"], tmp[, "DF"], tmp[nrow(tmp), "DF"], lower.tail=FALSE))
		tmp[nrow(tmp), c("F value", "Pr(>F)")] <- NA	
		res$aov.org <- tmp
	}
	
	res <- solveMME(res)
	
	if(allObsEqual)
		res$aov.tab[,"%Total"] <- 0	
	
	gc(verbose=FALSE)													# trigger garbage collection
	return(res)
}

#' ANOVA Type-I Degrees of Freedom.
#' 
#' Depending on the type of model, e.g. fully-nested, crossed-nested, etc. algorithms
#' are applied which are believed to be reasonably fast for the respective type of model, whenever
#' ANOVA sums of squares are constructed via quadradic forms in y (SSQ.method="qf").
#' 
#' This function is not meant to be called directly. It is invoked by functions \code{\link{anovaVCA}}
#' and \code{\link{anovaMM}}.
#' 
#' Computing the rank of a corresponding \code{A}-matrix, which generates ANOVA sum of squares as
#' quadratic form in \eqn{y} is a general method applicable to all types of models. This usually
#' envokes a singular value decomposition of \eqn{A} which makes it rather slow. 
#' Here, we try to speed up things by differentiating three classes of models, 1) fully-nested 
#' models where DFs are computed as the number of factor-levels minus the sum of higher order terms
#' minus 1, 2) models with only main factors (in this case \code{\link{anova.lm}} is used), 3) models
#' with main factors and interaction terms. Main factors DFs are computed employing function \code{\link{anova.lm}}.
#' DFs for interaction terms are computed by determining the rank of \eqn{A}-matrices. There are certain 
#' designs for which \code{anova.lm} becomes very fast (see examples section of \code{\link{anovaMM}}).
#' 
#' @param form			(formula) object specifying the model for which ANOVA DFs are requested
#' @param Data			(data.frame) with all variables appearing in 'form'
#' @param Zmat			(list) of Z-matrices, representing the design matrices of all model-terms, interpreted 
#' 						as fixed effects (number of columns represent number of unique levels)
#' @param Amat			(list) of A-matrices, generating ANOVA sums of squares as quadratic forms in \eqn{y},
#'                      see \code{\link{anovaVCA}} for details
#' @param tol			(numeric) constant, representing a numeric tolerance used in computing the rank of
#'                      \eqn{A}-matrices
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples
#' \dontrun{
#' # fully-nested design
#' data(realData)
#' datP1 <- realData[realData$PID==1,]
#' system.time(anova.lm.Tab1 <- anova(lm(y~lot/calibration/day/run, datP1)))
#' system.time(anovaMM.Tab1  <- anovaMM(y~lot/calibration/day/run, datP1))
#' anova.lm.Tab1
#' anovaMM.Tab1
#' 
#' # use SSQ.method="qf" (based on quadratic forms)
#' system.time(anovaMM.Tab1.qf  <- anovaMM(y~lot/calibration/day/run, datP1, SSQ.method="qf"))
#' 
#' # compute degrees of freedom
#' VCA:::anovaDF( y~lot/calibration/day/run, datP1,
#' 				  Zmat=anovaMM.Tab1.qf$Matrices$Z,
#' 				  Amat=anovaMM.Tab1.qf$Matrices$A)
#' 
#' # design with only main-factors
#' system.time(anova.lm.Tab2 <- anova(lm(y~lot+calibration+day+run, datP1)))
#' system.time(anovaMM.Tab2  <- anovaMM(y~lot+calibration+day+run, datP1))
#' anova.lm.Tab2
#' anovaMM.Tab2
#' 
#' # use SSQ.method="qf" (based on quadratic forms)
#' system.time(anovaMM.Tab2.qf  <- anovaMM(y~lot+calibration+day+run, datP1, SSQ.method="qf"))
#' 
#' # compute degrees of freedom
#' VCA:::anovaDF( y~lot+calibration+day+run, datP1,
#' 				  Zmat=anovaMM.Tab2.qf$Matrices$Z,
#' 				  Amat=anovaMM.Tab2.qf$Matrices$A)
#' 
#' # design with main-factors and interactions
#' system.time(anova.lm.Tab3 <- anova(lm(y~(lot+calibration)/day/run, datP1)))
#' system.time(anovaMM.Tab3  <- anovaMM( y~(lot+calibration)/day/run, datP1))
#' anova.lm.Tab3
#' anovaMM.Tab3
#' 
#' # use SSQ.method="qf" (based on quadratic forms)
#' system.time(anovaMM.Tab3.qf  <- anovaMM(y~(lot+calibration)/day/run, datP1, SSQ.method="qf"))
#' 
#' # compute degrees of freedom
#' VCA:::anovaDF( y~(lot+calibration)/day/run, datP1,
#' 				  Zmat=anovaMM.Tab3.qf$Matrices$Z,
#' 				  Amat=anovaMM.Tab3.qf$Matrices$A)
#' }

anovaDF <- function(form, Data, Zmat, Amat, tol=1e-8)
{
	form <- terms(form, simplify=TRUE, keep.order=TRUE)
	
	fac  <- attr(form, "term.labels")
	
	if(length(fac) == 0)										# handle pure intercept models
		return(nrow(Data)-1)
	fmat <- attr(form, "factors")[-1,,drop=FALSE]				# remove row corresponding to the response
	csum <- apply(fmat, 2, function(x) length(which(x>0)))
	fn   <- all(diff(csum) == 1)								# fully-nested model? ( i-th term part of (i+1)-th term --> csum[i]+1 == csum[i+1] )
	
	DF <- numeric(length(fac))
	Split <- list()
	
	Nlvl <- sapply(Zmat, ncol)									# number of unique factor levels
	names(Nlvl) <- fac
	
	if(fn)														# very fast DF-algorithm applicable
	{
		for(i in 1:length(fac))
		{
			if(i == 1)
				DF[i] <- Nlvl[i] - ifelse(attr(form, "intercept"), 1, 0)
			else
				DF[i] <- Nlvl[i] - sum(DF[1:(i-1)]) - ifelse(attr(form, "intercept"), 1, 0)
		}
	}
	else														# not fully-nested
	{
		ia <- csum > 1											# interactions? 
		mf <- !ia												# main factors
		
		if(!any(ia))											# only main-factors in the model
		{
			DF <- anova(lm(form, Data))
			DF <- DF[-nrow(DF),"Df"]
		}
		else													# main factors and interactions
		{
			if(any(mf))											# get DFs for all main factors
			{
				resp <- as.character(form)[2]					# response variable
				
				tmp <- anova(lm(paste(	resp, "~", paste(fac[which(mf)], collapse="+"), 
										ifelse(attr(form, "intercept"),"", "-1"), sep=""), Data))
				
				DF[which(mf)]  <- tmp[-nrow(tmp), "Df"]
			}
			
			if(any(ia))											# apply computationally expensive DF-algorithm on interactions
			{
				for(i in which(ia))
					DF[i] <- rankMatrix(Amat[[i]], tol=tol)
			}			
		}
	}
	
	DF <- c(DF, nrow(Data) - sum(DF) - ifelse(attr(form, "intercept"), 1, 0))									# residual degrees of freedom
	names(DF) <- c(fac, "error")
	attr(DF, "Nlevel") <- Nlvl
	return(DF)
}


#' Solve Mixed Model Equations.
#' 
#' Function solves the Mixed Model Equations (MME) to estimate fixed and random effects.
#' 
#' This function is for internal use only, thus, not exported.
#' 
#' @param obj			... (VCA) object
#' 
#' @return 	(VCA) object, which has additional elements "RandomEffects" corresponding to the column vector 
#' 		   	of estimated random effects, "FixedEffects" being the column vector of estimated fixed effects. 
#' 			Element "Matrices" has additional elements referring to the elements of the MMEs and element
#' 			"VarFixed" corresponds to the variance-covariance matrix of fixed effects.
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples
#' \dontrun{
#' data(dataEP05A2_1)
#' fit <- anovaVCA(y~day/run, dataEP05A2_1, NegVC=TRUE)
#' fit <- solveMME(fit)
#' ranef(fit)
#' }

solveMME <- function(obj)
{
	stopifnot(class(obj) == "VCA")
	mats   	<- obj$Matrices
	V      	<- mats[["V"]]
	
	if(is.null(V))						# will be the case for ANOVA-type fitted models
	{
		obj  <- getV(obj)
		mats <- obj$Matrices
		V	 <- mats[["V"]]
	}
	Z      	<- mats$Zre
	R 	   	<- mats$R
	G      	<- mats$G
	X      	<- Matrix(mats$X)
	y      	<- mats$y
	
	if(is.null(mats$Vi))
		Vi  <- Solve(V, quiet=TRUE)
	else
		Vi	<- mats$Vi
	
	K 		<- Solve(t(X) %*% Vi %*% X, quiet=TRUE)		# variance-covariance matrix of fixed effects
	T	   	<- K %*% t(X) %*% Vi
	fixed  	<- T %*% y
	
	mats$Vi <- Vi
	mats$T  <- T
	rownames(fixed) <- colnames(X)
	colnames(fixed) <- "Estimate"
	re		<- obj$RandomEffects
	
	if(is.null(Z))
		re <- NULL
	else
	{
		if(is.null(re))
		{
			re  <- G %*% t(Z) %*% Vi %*% (y - X %*% fixed) 
			rownames(re) <- colnames(Z)
			colnames(re) <- "Estimate"
		}
	}
	obj$RandomEffects <- re
	obj$FixedEffects  <- fixed
	obj$Matrices	  <- mats
	obj$VarFixed      <- K
	return(obj)
}


#' Generic Method for Extracting Random Effects from a Fitted Model.
#' @param object		(object)
#' @param ...			additional parameters
#' @seealso \code{\link{ranef.VCA}}

ranef <- function(object, ...)
	UseMethod("ranef")

#' Extract Random Effects from 'VCA' Object.
#' 
#' Extract random effects and possibly apply a transformation to them (standardization,
#' studentization).
#' 
#' Extracting the 'RandomEffects' element of an 'VCA' object if this exists and applying
#' standardization (mean 0, sd 1) or studentization. For studentized random effects 
#' the i-th random effects is divided by the i-th main diagonal element of matrix \eqn{O = GZ^{T}QZG}{O = GZ'QZG},
#' where \eqn{G} is the covariance-matrix of random effects, \eqn{Z} is a design matrix assigning 
#' random effects to observations and matrix \eqn{Q = V^{-1}(I - H)}{Q = V"(I - H)} (see \code{\link{residuals.VCA}} for further details). 
#' 
#' @param object		(VCA) object from which random effects shall be extracted
#' @param term			(character) string specifying a term (factor) for which random effects 
#'                      should be extracted, one can also specify an integer which is interpreted
#'                      as i-th element of 'obj$res.assign$terms'
#' @param mode			(character) string or abbreviation specifying whether "raw" residuals
#'                      should be returned or a transformed version c("student" or "standard")
#' @param quiet			(logical) TRUE = will suppress any warning, which will be issued otherwise 
#' @param ...			additional parameters
#' 
#' @method ranef VCA
#' 
#' @references 
#' 
#' Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York	
#' 
#' Laird, N.M., Ware, J.H., 1982. Random effects models for longitudinal data. Biometrics 38, 963-974.
#' 
#' Schuetzenmeister, A. and Piepho, H.P. (2012). Residual analysis of linear mixed models using a simulation approach.
#' Computational Statistics and Data Analysis, 56, 1405-1416
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_1)
#' fit <- anovaVCA(y~day/run, dataEP05A2_1)
#' ranef(fit)
#' 
#' # get variable-specific random effects (REs)
#' # both extract the same REs
#' ranef(fit, "day")
#' ranef(fit, 1)
#' 
#' # get standardized REs
#' ranef(fit, "day:run", "standard")
#' 
#' # or studentized REs
#' ranef(fit, 2, "stu")
#' }

ranef.VCA <- function(object, term=NULL, mode=c("raw", "student", "standard"), quiet=FALSE, ...)
{
	Call <- match.call()
	
	obj <- object
	
	if(is.list(obj) && class(obj) != "VCA")
	{
		if(!all(sapply(obj, class) == "VCA"))
			stop("Only lists of 'VCA' object are accepted!")
		
		obj.len <- length(obj)
		
		if(!"msgEnv" %in% ls(.GlobalEnv))
			msgEnv <<- new.env(parent=emptyenv())
		
		assign("VCAinference.obj.is.list", TRUE, envir=msgEnv)			# indicate that a list-type object was passed intially
		
		if(is.null(term))
		{
			res <- mapply(	FUN=ranef.VCA, object=obj,
					mode=mode[1], quiet=quiet, SIMPLIFY=FALSE)
		}
		else
		{
			res <- mapply(	FUN=ranef.VCA, object=obj, term=term,
					mode=mode[1], quiet=quiet, SIMPLIFY=FALSE)
		}
		names(res) <- names(obj)
		
		if(obj.len == 1)			# mapply returns a list of length 2 in case that length(obj) was equal to 1
			res <- res[1]
		
		rm("VCAinference.obj.is.list", envir=msgEnv)
		
		return(res)
	}	
	
	stopifnot(class(obj) == "VCA")
	mode <- match.arg(mode)
	
	if(!is.null(obj$scale) && is.null(obj$rescaled))
		warning("The fitted model has not been re-scaled yet! Results are likely to differ from correct results!")
	
	if(!is.null(obj$scale) && !mode %in% c("raw", "standard"))
		scale <- obj$scale
	else
		scale <- 1
	
	ObjNam  <- deparse(Call$object)
	ObjNam2 <- sub("\\[.*", "", ObjNam)
	
	if(is.null(obj$RandomEffects))
	{
		obj  <- solveMME(obj)
		
		if(length(ObjNam2) == 1 && ObjNam2 %in% names(as.list(.GlobalEnv)))
		{
			expr <- paste(ObjNam, "<<- obj")		# update object missing MME results
			eval(parse(text=expr))
		}
		else
		{
			if( !"VCAinference.obj.is.list" %in% names(as.list(msgEnv)) && !quiet)
				message("Mixed model equations were solved but results could not be assigned to 'VCA' object!")
		}
	}
	
	if(mode == "student" && is.null(obj$Matrices$Q))
	{
		mats <- obj$Matrices
		X  <- mats$X
		T  <- mats$T
		Vi <- mats$Vi
		mats$H <- H  <- X %*% T
		mats$Q <- Q  <- Vi %*% (diag(nrow(H))-H)
		obj$Matrices <- mats
		
		if(length(ObjNam2) == 1 && ObjNam2 %in% names(as.list(.GlobalEnv)))
		{
			expr <- paste(ObjNam, "<<- obj")		# update object missing MME results
			eval(parse(text=expr))
		}
		else
		{
			if( !"VCAinference.obj.is.list" %in% names(as.list(msgEnv)) && !quiet)
				message("Matrices 'H' and 'Q' were comuted but could not be assigned to 'VCA' object!")
		}
	}
	
	re <- obj$RandomEffects
	nam <- rownames(re)
	if(!is.null(term) && !is.character(term))
	{
		term <- obj$re.assign$terms[as.integer(term)]
	}
	
	if(mode == "standard")
	{
		tmp <- tapply(re, obj$re.assign$ind, scale)
		re  <- matrix(nrow=nrow(re))
		for(i in 1:length(obj$re.assign$terms))
			re[which(obj$re.assign$ind == i)] <- tmp[[i]]
		
		rownames(re) <- nam
	}
	else if(mode == "student")
	{
		G <- getMat(obj, "G")
		Q <- getMat(obj, "Q")
		Z <- getMat(obj, "Z")
		O <- G %*% t(Z) %*% Q %*% Z %*% G
		re <- re / sqrt(diag(O))
		re <- as.matrix(re)
		rownames(re) <- nam
	}
	
	re <- re/scale
	
	if(is.null(term) || !term %in% obj$re.assign$terms)
	{
		if(!is.null(term) && !term %in% obj$re.assign$terms && !quiet)
		{
			warning("There is no term in the random part of the formula corresponding to specified 'term'!")
		}
		
		ind <- NULL
		for(i in 1:length(obj$re.assign$terms))
			ind <- c(ind, which(obj$re.assign$ind == i))
		
		re <- re[ind,,drop=FALSE]
		
		attr(re, "mode") <- mode
		attr(re, "term") <- "all"
		
		return(re)
	}
	else
	{
		ind <- which(obj$re.assign$ind == which(obj$re.assign$terms == term)) 
		re  <- re[ind,,drop=F]
		attr(re, "mode") <- mode
		attr(re, "term") <- term
		
		return(re)
	}
}



#' Generic Method for Extracting Fixed Effects from a Fitted Model.
#' @param object		(object)
#' @param ...			additional parameters
#' @seealso \code{\link{fixef.VCA}}

fixef <- function(object, ...)
	UseMethod("fixef")

#' Extract Fixed Effects from 'VCA' Object.
#' 
#' Conveniently extracting the 'FixedEffects' element of an 'VCA' object. 
#' 
#' The default is to return the fixed effects estimates together with their standard errors.
#' If setting 'type="complex"' or to an abbreviation (e.g. "c") additional inferential statistics
#' on these estimates will be returned, i.e. "t Value", "DF" and respective p-value "Pr > |t|". 
#' One can choose one of three denominator degrees of freedom ('ddfm')-methods. The implementation
#' of these methods are an attempt to align with the results of SAS PROC MIXED. See the respective
#' SAS-documentation for details.
#' 
#' @param object		(VCA) object where fixed effects shall be extracted
#' @param type			(character) string or partial string, specifying whether
#'                      to return "simple" (reduced) or a rather "complex" (more detailed) 
#'                      information about fixed effects
#' @param ddfm			(character) string specifying the method used for computing the 
#'                      degrees of freedom of the t-statistic. Only used when type="complex".
#' 						Available methods are "contain", "residual", and "satterthwaite".
#' @param tol			(numeric) value representing the numeric tolerance use in comparisons, values
#' 						smaller than 'tol' will be considered equal to 0
#' @param quiet			(logical) TRUE = suppress warning messages, e.g. for non-estimable contrasts
#' @param ...			additional parameters
#' 
#' @method fixef VCA
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_1)
#' fit <- anovaVCA(y~day/(run), dataEP05A2_1)
#' fixef(fit)
#' 
#' # for complex models it might take some time computing complex output
#' data(VCAdata1)
#' fit <- anovaMM(y~(lot+device)/(day)/(run), VCAdata1[VCAdata1$sample==2,])
#' fixef(fit, "c")
#' }

fixef.VCA <- function(	object, type=c("simple", "complex"), ddfm=c("contain", "residual", "satterthwaite"), 
						tol=1e-12, quiet=FALSE, ...)
{
	Call <- match.call()
	
	obj <- object
	
	if(is.list(obj) && class(obj) != "VCA")
	{
		if(!all(sapply(obj, class) == "VCA"))
			stop("Only lists of 'VCA' object are accepted!")
		
		obj.len <- length(obj)
		
		if(!"msgEnv" %in% ls(.GlobalEnv))
			msgEnv <<- new.env(parent=emptyenv())
		
		assign("VCAinference.obj.is.list", TRUE, envir=msgEnv)			# indicate that a list-type object was passed intially
		
		res <- mapply(	FUN=fixef.VCA, obj=obj, type=type[1],
						ddfm=ddfm[1], tol=tol, quiet=quiet,
						SIMPLIFY=FALSE)
		
		names(res) <- names(obj)
		
		if(obj.len == 1)			# mapply returns a list of length 2 in case that length(obj) was equal to 1
			res <- res[1]
		
		rm("VCAinference.obj.is.list", envir=msgEnv)
		
		return(res)
	}	
	
	stopifnot(class(obj) == "VCA")
	type <- match.arg(type)
	if(length(ddfm) > 1 && type == "complex")
	{
		ddfm <- "satterthwaite"
		if(!quiet)
			message("Note: 'ddfm' not specified, option \"satterthwaite\" was used!")
	}
	ddfm <- match.arg(ddfm)
	
	vc <- vcov(obj, quiet=quiet)
	se <- suppressWarnings(sqrt(diag(vc)))
	fe <- obj$FixedEffects
	
	if(is.null(fe))								# solve mixed model equations first
	{
		obj  <- solveMME(obj)
		fe   <- obj$FixedEffects
		nam0 <- deparse(Call$object)
		
		nam1 <- sub("\\[.*", "", nam0)			# remove any index-operators 
		if(length(nam1) == 1 && nam1 %in% names(as.list(.GlobalEnv)))
		{
			expr <- paste(nam0, "<<- obj")		# update object missing MME results
			eval(parse(text=expr))
		}
		else			# message only if not called on list of VCA-objects
		{
			if( !"VCAinference.obj.is.list" %in% names(as.list(msgEnv)) && !quiet)
				message("Some required information missing! Usually solving mixed model equations has to be done as a prerequisite!")
		}
	}
	nam <- rownames(fe)
	
	if(type == "simple")
	{		
		fe <- matrix(fe, ncol=1)
		colnames(fe) <- "Estimate"
		fe <- cbind(fe, SE=se)		
		rownames(fe) <- nam
	}
	else
	{
		if(is.null(obj$VarCov))
			obj$VarCov <- vcovVC(obj)
		if(is.null(obj$VarFixed))
			obj$VarFixed <- vcov(obj)
		LC  <- diag(nrow(fe))	
		rownames(LC) <- rownames(fe)
		colnames(LC) <- rownames(fe)
		tst <- test.fixef(obj, LC, ddfm=ddfm, quiet=TRUE)
		tst <- cbind(tst, SE=se)
		fe  <- tst[,c(1,5,2,3,4),drop=F]
		NAs <- is.na(fe[,"Estimate"])
		if(any(NAs))
		{
			fe[which(NAs),"Estimate"] <- 0
			fe[which(NAs), 2:5] <- NA
		}
	}
	
	return(fe)
}




#' Least Squares Means of Fixed Effects.
#' 
#' Computes Least Squares Means (LS Means) of fixed effects for fitted mixed models of class 'VCA'.
#' 
#' Function computes LS Means of fixed effects and their corresponding 
#' standard errors. In case of setting argument 'type' equal to "complex" (or any
#' abbreviation) a \eqn{t}-test is performed on each LS Mean, returning degrees
#' of freedom, t-statistic and corresponding p-values. One can choose from one of three
#' denominator degrees of freedom ('ddfm')-methods.
#' 
#' Actually, function \code{\link{test.fixef}} is called with the "no intercept" 
#' version of the fitted model. The "complex" option is significantly slower for unbalanced
#' designs (see \code{\link{test.fixef}} for details). In case that the 'VarCov' element of
#' the 'VCA' object already exists (calling \code{\link{vcovVC}}), which is the most time 
#' consuming part, results can be obtained in less amount of time.
#' 
#' Standard Errors of LS Means are computed as \eqn{TPT^{T}}{T * P * T'}, where \eqn{T} is the
#' LS Means generating contrast matrix and \eqn{P} is the variance-covariance matrix of
#' fixed effects.
#' 
#' Argument \code{at} can be used to modify the values of covariables when computing LS Means and/or
#' to apply different weighting schemes for (fixed) factor varialbes in the model, e.g. when the prevelance
#' of factor-levels differs from a uniform distribution. Usually, if the weighting scheme is not modified,
#' each factor-level will contribute \eqn{1/N} to the LS Mean, where \eqn{N} corresponds to the number of factor-levels. 
#' 
#' Covariables have to be specified as 'name=value', where value can be a vector of length > 1. 
#' Each value will be evaluated for each row of the original LS Means contrast matrix. 
#' If multiple covariables are specified, the i-th element of covariable 1 will be matched with
#' the i-th element of covariable(s) 2...M, where \eqn{M} is the number of covariables in the model.
#' 
#' To apply a different weighting scheme for factor-variables one has to specify 'factor-name=c(level-name_1=value_1,
#' level-name_2=value_2, ..., level-name_N=value_N)'. The sum of all 'value_i' elements must be equal to 1, otherwise,
#' this factor-variable will be skipped issuing a warning. If any levels 'level-name_i' cannot be found for 
#' factor-variable 'factor-name', this variable will also be skipped and a warning will be issued.
#' See the examples section to get an impression of how this works.
#' 
#' @param obj			(VCA) object having at least one fixed effect
#' @param var			(character) string specifying a fixed effects variable for which
#'                      LS Means should be computed, defaults to all fixed effects, i.e. for
#'         				each level of a fixed effects variable ls means will be computed
#' @param type			(character) "simple" = fast version of computing LS means
#' @param ddfm			(character) string specifying the method used for computing the 
#'                      degrees of freedom of the t-statistic. Only used when type="complex".
#' 						Available methods are "contain", "residual", and "satterthwaite".
#' @param at			(list) where each element corresponds either to a (numeric) covariable or
#' 						to a factor-variable for which the weighting scheme should be adjusted.
#' 						See details section for a thorough description of how argument 'at' works
#' 						and also see the examples.
#' @param contr.mat		(logical) TRUE = the LS Means generating contrast-matrix will be added to the
#' 						result as attribute \code{contrasts}
#' @param quiet			(logical) TRUE = suppress warning messages, e.g. for non-estimable contrasts
#' 
#' @return (matrix) with LS Means of fixed effects and respective standard errors,
#'         in case of 'type="complex"'
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples #
#' \dontrun{
#' data(dataEP05A2_2)
#' fit1 <- anovaMM(y~day/(run), dataEP05A2_2)
#' lsmeans(fit1)
#' lsmeans(fit1,, "complex")
#' 
#' # a more complex model
#' data(VCAdata1)
#' fit2 <- anovaMM(y~(lot+device)/(day)/(run), VCAdata1[VCAdata1$sample==2,])
#' lsmeans(fit2, "lot")
#' lsmeans(fit2, "device", "complex")
#' 
#' # pre-computed 'VarCov' element saves time
#' system.time(lsm1 <- lsmeans(fit2, "device", "complex"))
#' fit2$VarCov <- vcovVC(fit2)
#' system.time(lsm2 <- lsmeans(fit2, "device", "complex"))
#' lsm1
#' lsm2 
#' 
#' # simulate some random data 
#' set.seed(212)
#' id <- rep(1:10,10)
#' x <- rnorm(200)
#' time <- sample(1:5,200,replace=T)
#' y <- rnorm(200)+time
#' snp <- sample(0:1,200,replace=T)
#' dat <- data.frame(id=id,x=x,y=y,time=time,snp=snp)
#' dat$snp <- as.factor(dat$snp)
#' dat$id <- as.factor(dat$id)
#' dat$time <- as.numeric(dat$time)
#' dat$sex <- gl(2, 100, labels=c("Male", "Female"))
#' dat$y <- dat$y + rep(rnorm(2, 5, 1), c(100, 100))
#' 
#' fit3 <- remlMM(y~snp+time+snp:time+sex+(id)+(id):time, dat)
#' 
#' # comute standard LS Means for variable "snp"
#' lsmeans(fit3, var="snp")
#' lsmeans(fit3, var="snp", type="c")    # comprehensive output
#' 
#' # compute LS Means at timepoints 1, 2, 3, 4
#' # Note: original LS Means are always part of the output
#' lsmeans(fit3, var="snp", at=list(time=1:4))
#' 
#' # compute LS Means with different weighting scheme
#' # for factor-variable 'sex'
#' lsmeans(fit3, var="snp", at=list(sex=c(Male=.3, Female=.7)))
#' 
#' # combine covariables at some value and altering the
#' # weighting scheme
#' lsmeans(fit3, var="snp", at=list(time=1:4, sex=c(Male=.3, Female=.7)))
#' 
#' # now with comprehensive output and requesting the
#' # LS Means generating contrast matrix
#' lsmeans(fit3, var="snp", type="complex", contr.mat=TRUE,
#'         at=list(time=1:4, sex=c(Male=.3, Female=.7)))
#' }

lsmeans <- function(obj, var=NULL, type=c("simple", "complex"), ddfm=c("contain", "residual", "satterthwaite"), 
					at=NULL, contr.mat=FALSE, quiet=FALSE)
{
	Call <- match.call()
	
	if(is.list(obj) && class(obj) != "VCA")
	{
		if(!all(sapply(obj, class) == "VCA"))
			stop("Only lists of 'VCA' object are accepted!")
		
		obj.len <- length(obj)
		
		if(is.null(var))
		{
			res <- mapply(	FUN=lsmeans, obj=obj, 
					type=type[1], ddfm=ddfm[1], quiet=quiet,
					SIMPLIFY=FALSE)
		}
		else
		{
			res <- mapply(	FUN=lsmeans, obj=obj, var=var, 
					type=type, ddfm=ddfm, quiet=quiet,
					SIMPLIFY=FALSE)
		}
		
		names(res) <- names(obj)
		
		if(obj.len == 1)			# mapply returns a list of length 2 in case that length(obj) was equal to 1
			res <- res[1]
		
		return(res)
	}	
	
	stopifnot(class(obj) == "VCA")
	stopifnot(obj$Type %in% c("Linear Model", "Mixed Model"))		# won't work for random models
	
	type <- match.arg(type)
	
	if(!is.null(var))												# does this fixed effects variable exist
		stopifnot(var %in% obj$fixed)
	
	if(length(ddfm) > 1 && type == "complex")
	{		
		ddfm <- "satterthwaite"
		if(!quiet)
			message("Note: 'ddfm' not specified, option \"satterthwaite\" was used!")
	}
	ddfm <- match.arg(ddfm)
	
	if(obj$EstMethod == "REML" && ddfm == "contain" && type=="complex")
	{
		ddfm <- "satterthwaite"
		if(!quiet)
			message("Note: 'ddfm' set to \"satterthwaite\", currently option \"contain\" does not work for REML-estimation!")
	}
	
	if(obj$Type != "Mixed Model" && !obj$intercept)
	{
		if(!quiet)
			warning("There are no fixed effects for which LS Means could be calculated!")
		return(NA)
	}
	
	T <- lsmMat(obj, var=var, quiet=quiet)					# LS Means generating contrast matrix
	
	id.mat <- NULL											# being non-NULL means that something was done via 'at'
	
	if(!is.null(at))										# LS Means at some fixed values of covariates and/or with different weighting scheme for fixed factor-variables
	{
		Means <- attr(T, "means")
		
		dat.class <- attr(T, "var.class")
		
		nam <- names(at)
		
		T2  <- T
		cnT <- colnames(T)
		fe.terms <- attr(obj$fe.assign, "terms")
		
		num.at <- nam[dat.class[nam] == "numeric"]
		fac.at <- nam[dat.class[nam] == "factor"]
		
		if(all(is.na(c(num.at, fac.at))) && !quiet)
			warning("Argument 'at' was not correctly specified! Neither covariables nor factor variables could be matched!")
		
		if(length(num.at) == 1 && is.na(num.at))								
			num.at <- character()							# ensure feasibility tests below
		if(length(fac.at) == 1 && is.na(fac.at))
			fac.at <- character()
		
		if(length(num.at) > 0)								# if there are any covariables
		{
			at.mat <- matrix(nrow=max(sapply(at[num.at], length)), ncol=length(num.at))
			colnames(at.mat) <- num.at
			
			for(i in 1:length(num.at))						# fill matrix, could also be user-defined
			{
				if(!num.at[i] %in% cnT)
				{
					if(!quiet)
						warning("There is no fixed term ", paste0("'", nam[i],"'!"),"Possible terms are:", paste(cnT, collapse=", "))
					next
				}
				
				if(length(at[[num.at[i]]]) < nrow(at.mat))
				{
					if(!quiet)
						warning("at[[",num.at[i],"]] does not match the max-length element of 'at', it will be replicated as necessary!")
					
					at.mat[,i] <- rep(at[[num.at[i]]], ceiling(nrow(at.mat)/length(at[[num.at[i]]])))[1:nrow(at.mat)]	# replicate
				}
				else
					at.mat[,i] <- at[[num.at[i]]]
			}
		}
		else
			at.mat <- matrix(nrow=length(fac.at), ncol=0)					# empty matrix
		
		if(length(fac.at) > 0)										# treat factor-variables, i.e. different weighting scheme
		{
			fac.lst <- vector("list", length(fac.at))
			names(fac.lst) <- fac.at
			
			for(i in 1:length(fac.at))								# over all specified factor variables
			{	
				if(sum(at[[fac.at[i]]]) != 1)
				{
					if(!quiet)
						warning("Sum of all coefficients of factor-variable ", paste0("'", fac.at[i],"'"), " is not equal to 1! It will be skipped!")
					next
				}
				
				skip <- FALSE
				
				tmp.mat <- matrix(nrow=nrow(at.mat), ncol=0)
				
				for(j in 1:length(at[[fac.at[i]]]))					# over all levels of the i-th factor variable
				{
					tmp <- matrix(at[[fac.at[i]]][j], ncol=1, nrow=nrow(at.mat))
					colnames(tmp) <- paste0(fac.at[i], names(at[[fac.at[i]]])[j])
					
					if(!colnames(tmp) %in% colnames(T))
					{
						if(!quiet)
							warning("Factor-level ", paste0("'",colnames(tmp), "'"),
									" of variable ",paste0("'", fac.at[i], "'"),
									" does not correspond to a fixed effect!\n  This element of 'at' will be skipped!")
						skip <- TRUE
					}
					tmp.mat <- cbind(tmp.mat, tmp)
				}
				
				if(skip)
					next
				else
				{
					at.mat <- cbind(at.mat, tmp.mat)
					fac.lst[[i]] <- tmp.mat
				}
			}
		}
		at.mat <- na.omit(at.mat)
		
		id.mat <- matrix(nrow=nrow(T), ncol=ncol(at.mat))
		cn.id  <- colnames(at.mat)
		colnames(id.mat) <- cn.id
		
		if(length(num.at) > 0)
		{
			for(i in 1:length(num.at))								# identifies combination of covar-levels
				id.mat[,num.at[i]] <- Means[[num.at[i]]]["mean"]
			
			tmp.ind <- cn.id[!cn.id %in% num.at]					# those columns not corresponding to covariables
		}
		else
			tmp.ind <- cn.id										# only columns corresponding to factor variable weights
		
		if(length(fac.at) > 0)
			id.mat[,tmp.ind] <- T[,tmp.ind]
		
		if(ncol(at.mat) > 0)								# only if there is anything to evaluate
		{
			for(i in 1:nrow(at.mat))						# for each combination specified
			{
				T3  <- T
				
				if(length(num.at) > 0)						# if there any covariables
				{
					for(j in 1:length(num.at))				# over covariables
					{
						cnTi <- cnT[grepl(num.at[j], cnT)]					# all relevant columns in matrix T					
						
						for(k in 1:length(cnTi))							# over all columns in which the current covariable appears (main-effect or interaction)
						{
							if(grepl(":", cnTi[k]))							# interaction
							{
								splt <- unlist(strsplit(cnTi[k], ":"))		# split interaction into atomic terms
								
								tmp.val <- 1
								
								for(l in 1:length(splt))							# over all terms in the interaction
								{
									if(splt[l] == num.at[j])						# user-specified level	
										tmp.val <- tmp.val * at.mat[i,num.at[j]]
									else
									{
										if(splt[l] %in% names(Means))						# another numeric covariable -> use mean of this value
											tmp.val <- tmp.val * Means[[splt[l]]]["mean"]
										else												# check with rowname
										{
											if(splt[l] %in% rownames(T3))
											{
												tmp.fac <- rep(0, nrow(T3))
												tmp.fac[grepl(splt[l], rownames(T3))] <- 1		# only that LS Mean affected, which is part of the interaction in column cnTi[k]
												tmp.val <- tmp.val * tmp.fac
											}
											else										# not part of the current interaction
												tmp.val <- tmp.val * rep(0, nrow(T3))
										}
									}
								}
								T3[,cnTi[k]] <- tmp.val
							}
							else
								T3[,cnTi[k]] <- at.mat[i,num.at[j]]					# main effects		
						}					
					}
				}
				
				if(length(fac.at) > 0)
				{
					for(j in 1:length(fac.at))									# over factor variables
					{
						tmp.row <- fac.lst[[fac.at[j]]][i,]
						T3[, tmp.ind] <- rep(tmp.row, rep(nrow(T3), length(tmp.row)))  
					}
				}
				
				T2 <- rbind(T2, T3)
				
				tmp.id.mat <- matrix(rep(at.mat[i,], nrow(T)), ncol=ncol(at.mat), byrow=TRUE)
				colnames(tmp.id.mat) <- colnames(at.mat)
				
				id.mat <- rbind(id.mat, tmp.id.mat)
			}
		}		
		T <- T2										# T2 might be identical to T
	}
	
	attr(T, "var.class") <- NULL
	attr(T, "means") 	 <- NULL
	
	x   <- obj
	if(x$intercept)
	{
		x$intercept <- FALSE		
	}
	
	if(type == "complex")						# complex output --> slower
	{	
		if(is.null(obj$VarCov))
			obj$VarCov <- vcovVC(obj)			# determine vcov of VCs if missing
		lsm <- test.fixef(obj, T, ddfm=ddfm, quiet=TRUE, lsmeans=TRUE)
		rownames(lsm) <- rownames(T)
	}
	else										# simple output --> faster
	{
		vc  <- vcov(obj)	
		vc  <- T %*% vc %*% t(T)
		se  <- sqrt(diag(vc))
		lsm <- T %*% obj$FixedEffects
		lsm <- cbind(as.matrix(lsm), SE=se)
	}
	
	if(!is.null(id.mat) && ncol(id.mat) > 0)
		lsm <- cbind(id.mat, lsm)
	
	if(contr.mat)
		attr(lsm, "contrasts") <- T
	
	return(lsm)
}


#' Contrast Matrix for LS Means.
#' 
#' Function determines appropriate contrast matrix for computing the LS Means of
#' each factor level of one or multiple fixed effects variables. 
#' 
#' This functions implements the 5 rules given in the documentation of SAS PROC GLM for computing the LS Means.#' 
#' The LS Means correspond to marginal means adjusted for bias introduced by unbalancedness.
#' 
#' @param obj			(VCA) object
#' @param var			(character) string specifyig the fixed effects variable for which
#'                      the LS Means generating matrices should be computed
#' @param quiet			(logical) TRUE = will suppress any warning, which will be issued otherwise 
#' 
#' @return	(matrix) where each row corresponds to a LS Means generating contrast
#'          for each factor level of one or multiple fixed effects variable(s)
#' 
#' @author Andre Schutzenmeister \email{[email protected]@roche.com}
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_1)
#' fit1 <- anovaMM(y~day/run, dataEP05A2_1)
#' 
#' VCA:::lsmMat(fit1, "day")	# function not exported
#' VCA:::lsmMat(fit1, "run")
#' VCA:::lsmMat(fit1)			# is equal to listing all fixed terms
#' 
#' # a more complex and unbalanced model
#' data(VCAdata1)
#' datS1 <- VCAdata1[VCAdata1$sample == 1, ]
#' set.seed(42)
#' datS1ub <- datS1[-sample(1:nrow(datS1))[1:25],]
#' fit2 <- anovaMM(y~(lot+device)/day/(run), datS1ub)
#' VCA:::lsmMat(fit2, c("lot", "device"))
#' }

lsmMat <- function(obj, var=NULL, quiet=FALSE)
{
	stopifnot(class(obj) == "VCA")
	if(!is.null(var))
		stopifnot(var %in% obj$fixed)
	X <- getMat(obj, "X")
	
	if(is.null(var))
		var <- obj$fixed
	terms 	  <- attr(obj$terms, "factors")
	dat.cls   <- sapply(obj$data[,rownames(terms)[-1]], class)	
	variables <- names(dat.cls)
	fe.assign <- obj$fe.assign						# assignment by integers
	fe.terms  <- attr(fe.assign, "terms")	
	
	fe.class  <- list()
	
	for(i in 1:length(fe.terms))
	{
		if(fe.assign[i] == 0)
		{
			fe.class[[i]] <- NULL
		}
		tmp <- unlist(strsplit(fe.terms[i], ":"))
		fe.class[[i]] <- dat.cls[tmp]
	}
	names(fe.class) <- fe.terms
	
	if(!all(var %in% fe.terms))
		stop("At least one element of 'var' is not a fixed term!")
	
	fe.tvec   <- fe.terms[fe.assign+ifelse(obj$intercept, 1, 0)]				# assignment by names
	
	terms.cls <- rep("factor", length(fe.terms))
	names(terms.cls) <- fe.terms
	fe.tab    <- table(fe.tvec)
	
	lsmm <- matrix(nrow=0, ncol=ncol(X))				
	cn  <- colnames(lsmm) <- colnames(X)
	rn  <- Means <- NULL 
	
	if(any(dat.cls == "numeric"))					# find non-dummy variable columns in X
	{
		num.var <- names(dat.cls[which(dat.cls == "numeric")])
		for(i in 1:length(terms.cls))
		{
			if(any(sapply(num.var, grepl, fe.terms[i])))
			{
				terms.cls[i] <- "numeric"	
			}
		}		
	}
	
	fe.cls 		<- terms.cls[fe.assign+ifelse(obj$intercept, 1, 0)]
	num.terms 	<- terms.cls == "numeric"
	
	if(any(num.terms))														# are there any numeric terms
	{
		ind <- which(num.terms)
		Means <- list()
		
		for(i in 1:length(ind))
		{
			tmp		 <- numeric()
			tmp      <- c(tmp, Nlev=as.numeric(fe.tab[fe.terms[ind[i]]]))	# number of columns in X for the current numeric term
			splt     <- unlist(strsplit(fe.terms[ind[i]], ":"))
			tmp.cls  <- dat.cls[splt]
			num.var  <- which(tmp.cls == "numeric")
			fac.var  <- which(tmp.cls == "factor")
			tmp.mean <- apply(obj$data[,names(num.var), drop=FALSE], 1, function(x) eval(parse(text=paste(x, collapse="*"))))	# multiply each value of multiple numeric variables for each row		
			tmp 	 <- c(tmp, mean=mean(tmp.mean))							# mean of current combination of numeric variables
			
			if(length(fac.var) > 0)											# at least one factor variable
			{
				for(j in 1:length(fac.var))									# determine number of levels for each factor variable
				{
					exprs <- paste("c(tmp,",names(fac.var[j]),"=length(unique(obj$data[,\"",names(fac.var[j]),"\"])))", sep="") 
					tmp   <- eval(parse(text=exprs))
				}
			}
			eval(parse(text=paste("Means[[\"", fe.terms[ind[i]], "\"]] <- tmp", sep="")))	# add to list 'Means' and use name of the numeric variable as name of the list element
		}
	}
	
	for(i in 1:length(var))									# over all fixed terms for which LS Means shall be computed
	{
		tsplt <- unlist(strsplit(var[i], ":"))				# LS Means can only be generated for pure factor variables or combinations of such, no numeric variable may interact
		
		if(any(dat.cls[tsplt] == "numeric"))
		{
			if(!quiet)
				warning("'",var[i],"' is \"numeric\", LS Means cannot be estimated,'",var[i],"' will be skipped!")
			next
		}
		lvl.ind  <- which(fe.tvec == var[i])				# indices of all columns representing levels of var[i]
		lvl.nam  <- cn[lvl.ind]								# use column names of X instead of indices
		lvl.asgn <- unique(fe.assign[lvl.ind])				# value of the fixed effects assignment indicator for the current term var[i]
		tmp.lsm  <- matrix(0, nrow=0, ncol=ncol(X))
		colnames(tmp.lsm) <- cn
		rem.ind.init  <- which(fe.tvec != var[i])			# (init)ial (rem)aining columns which need to be treated differently
		rem.nam.init  <- cn[rem.ind.init]
		rem.asgn.init <- fe.assign[rem.ind.init] 
		
		for(j in 1:length(lvl.ind))							# over levels (columns) of the current factor
		{
			rem.ind  <- rem.ind.init						# re-set for each level of the current (i-th) term
			rem.nam  <- rem.nam.init
			rem.asgn <- rem.asgn.init
			
			con.mat <- matrix(0, nrow=1, ncol=ncol(X))
			colnames(con.mat) <- cn
			splt <- unlist(strsplit(lvl.nam[j], ":"))		# get all sub-terms
			
			if(obj$intercept)
			{
				con.mat[1,"int"] <- 1
				rem.ind  <- rem.ind[ -which(cn == "int")]
				rem.nam  <- rem.nam[ -which(cn == "int")]
				rem.asgn <- rem.asgn[-which(cn == "int")]
			}
			covar <- fe.cls[rem.ind] == "numeric"			# covariates involved?
			
			if(any(covar))									# [rule 1] (SAS PROC GLM documentation contrast matrix for LS Means)
			{
				cov.ind   <- rem.ind[ which(covar)]			# column-index of current covariate in contrast matrix
				rem.ind   <- rem.ind[-which(covar)]
				rem.nam   <- rem.nam[-which(covar)]
				rem.asgn  <- rem.asgn[-which(covar)]
				cov.asgn  <- fe.assign[cov.ind]				# covariate-indices
				ucov.asgn <- unique(cov.asgn)				
				
				for(k in 1:length(ucov.asgn))				# over all terms representing covariates
				{
					tmp.term <- fe.terms[ucov.asgn[k] + ifelse(obj$intercept, 1, 0)]
					tmp.info <- Means[[tmp.term]]
					tmp.ind  <- cov.ind[which(cov.asgn == ucov.asgn[k])]
					cn.splt  <- t(sapply(cn[tmp.ind], function(x) unlist(strsplit(x, ":"))))
					
					if(nrow(cn.splt) == 1)													# atomic covariate use mean value
					{
						con.mat[,which(fe.assign == ucov.asgn[k])] <- tmp.info["mean"]
						next
					}
					
					cov.cont <- apply(cn.splt, 1, function(x) any(splt %in% x))				# any sub-terms of the current term found in the levels of the current covariate
					
					if(any(cov.cont))														# covariate is distributed according to a term that is also a sub-term  of the current factor-level 
					{
						con.mat[1,tmp.ind[which(cov.cont)]] <- tmp.info["mean"]/length(which(cov.cont))
					}
					else																	# covariate independent of the current factor-level
					{
						con.mat[1,tmp.ind] <- tmp.info["mean"]/tmp.info["Nlev"]
					}
				}
			}
			
			tmp.splt  <- unlist(strsplit(lvl.nam[j], ":"))									# [rule 2] handling all effects that are contained by the current effect
			contained <- sapply(rem.nam, function(x) 
								all(unlist(strsplit(x, ":")) %in% tmp.splt))				
			
			if(any(contained))
			{
				tmp.asgn  <- rem.asgn[which(contained)]										# this might be multiple elements and this might be different values, e.g. 
				contained <- rem.nam[which(contained)]
				rem.ind   <- rem.ind[-which(rem.asgn %in% tmp.asgn)]						# column belongig to the same term must not be hanled again --> remove them from remaining elements
				rem.nam   <- rem.nam[-which(rem.asgn %in% tmp.asgn)]
				rem.asgn  <- rem.asgn[-which(rem.asgn %in% tmp.asgn)] 
				con.mat[,contained] <- 1
			}
			
			con.mat[1,lvl.ind[j]] <- 1														# [rule 3] setting the columns corresponding to the current effect to 1, all others remain 0
			
			contain <- sapply(	rem.nam, function(x) 											# [rule 4] consider effects that contain the current effect
								all(tmp.splt %in% unlist(strsplit(x, ":"))))  
			
			if(any(contain))
			{
				tmp.asgn  <- rem.asgn[which(contain)]
				utmp.asgn <- unique(tmp.asgn)
				
				for(k in 1:length(utmp.asgn))	
				{
					tmp.ind <- rem.ind[which(contain)]
					tmp.ind <- tmp.ind[which(tmp.asgn == utmp.asgn[k])]						# only those indices that belong to the k-th term (assignment value)
					con.mat[1, tmp.ind] <- 1 / length(tmp.ind)								
				}
				
				rem.ind  <- rem.ind[ -which(rem.asgn %in% utmp.asgn)]
				rem.nam  <- rem.nam[ -which(rem.asgn %in% utmp.asgn)]
				rem.asgn <- rem.asgn[-which(rem.asgn %in% utmp.asgn)]
			}
			
			if(length(rem.ind) > 0)															# [rule 5] there are still remaining, not yet treated effects
			{																				#          set these to 1/number of levels
				ufac.asgn <- unique(rem.asgn)
				
				for(k in 1:length(ufac.asgn))
				{
					con.mat[1, which(fe.assign == ufac.asgn[k])] <- 1/fe.tab[fe.terms[ufac.asgn[k] + ifelse(obj$intercept, 1, 0)] ]
				}
			}
			
			tmp.lsm <- rbind(tmp.lsm, con.mat)		
		}
		rownames(tmp.lsm) <- lvl.nam
		lsmm <- rbind(lsmm, tmp.lsm)
	}
	
	attr(lsmm, "var.class") <- dat.cls
	
	if(!is.null(Means))
		attr(lsmm, "means")	<- Means
	
	return(lsmm)
}


#' Extract Fixed Effects from 'VCA' Object.
#' 
#' For conveniently using objects of class 'VCA' with other packages expecting this
#' function, e.g. the 'multcomp' package for general linear hypotheses for parametric
#' models (currently not fully implemented).
#' 
#' @param object		(VCA) object where fixed effects shall be extracted
#' @param quiet			(logical) TRUE = will suppress any warning, which will be issued otherwise 
#' @param ...			additional parameters
#'
#' @method coef VCA 
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_1)
#' fit1 <- anovaMM(y~day/(run), dataEP05A2_1)
#' coef(fit1)
#' fit2 <- anovaVCA(y~day/run, dataEP05A2_1)
#' coef(fit2)
#' }

coef.VCA <- function(object, quiet=FALSE, ...)
{
	Call <- match.call()
	obj <- object
	stopifnot(class(obj) == "VCA")
	fe  <- fixef(obj)
	if(is.null(fe))								# solve mixed model equations first
	{
		obj  <- solveMME(obj)
		fe   <- fixef(obj)
		nam  <- as.character(as.list(Call)$object)
		if(length(nam) == 1 && nam %in% names(as.list(.GlobalEnv)))
		{
			expr <- paste(nam, "<<- obj")		# update object missing MME results
			eval(parse(text=expr))
		}
		else
		{
			if(!quiet)
				message("Some required information missing! Usually solving mixed model equations has to be done as a prerequisite!")
		}
	}
	fe  <- fe[,"Estimate", drop=F]
	nam <- rownames(fe)
	fe  <- c(fe)
	names(fe) <- nam
	return(fe)
}


#' Calculate Variance-Covariance Matrix of Fixed Effects for an 'VCA' Object.
#' 
#' Return the variance-covariance matrix of fixed effects for a linear mixed model
#' applicable for objects of class 'VCA'.
#' 
#' Actually this function only extracts this matrix or, if not available, calls function
#' \code{\link{vcovFixed}} which performs calculations. It exists for compatibility reasons,
#' i.e. for coneniently using objects of class 'VCA' with other packages expecting this
#' function, e.g. the 'multcomp' package for general linear hypotheses for parametric
#' models.  
#' 
#' @param object		 	(VCA) object for which the variance-covariance matrix of
#'                          fixed effects shall be calculated
#' @param quiet				(logical) TRUE = will suppress any warning, which will be issued otherwise 
#' @param ...				additional parameters
#' 
#' @return (matrix) corresponding to the variance-covariance matrix of fixed effects
#' 
#' @method vcov VCA
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_1)
#' fit1 <- anovaMM(y~day/(run), dataEP05A2_1)
#' vcov(fit1)
#' 
#' fit2 <- anovaVCAy~day/run, dataEP05A2_1)
#' vcov(fit2)
#' }

vcov.VCA <- function(object, quiet=FALSE, ...)
{
	obj <- object
	stopifnot(class(obj) == "VCA")
	if(!obj$intercept && length(obj$fixed) == 0)
	{
		if(!quiet)
			warning("There is no variance-convariance matrix of fixed effects for this object!")
		return(NA)
	}
	if(is.null(obj$VarFixed))
		return(vcovFixed(obj, quiet=quiet))
	else
		return(obj$VarFixed)
}


#' Extract Degrees of Freedom from Linear Hypotheses of Fixed Effects or LS Means.
#' 
#' Determine degrees of freedom for custom linear hypotheses of fixed effects or LS Means
#' using one of three possible approximation methods. 
#' 
#' This is a convenience function to determine the DFs for linear hypotheses in the same way
#' as function \code{\link{test.fixef}}. Only the "DF" part is returned here which can be passed
#' to other functions expecting DFs as input.
#' 
#' @param obj			(VCA) object
#' @param L				(matrix) specifying one or multiple linear hypothese, as returned by function
#'                      \code{\link{getL}} 
#' @param method		(character) the method to be used to determine the degrees of freedom for a 
#'                      linear hypothesis
#' @param ...			additional parameters
#' 
#' @return (numeric) vector with the DFs for each row of 'L'
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples 
#' \dontrun{
#' data(VCAdata1)
#' tmpDat <- VCAdata1[VCAdata1$sample==1,]
#' tmpDat <- tmpDat[-c(11,51,73:76),]
#' fitMM <- anovaMM(y~(lot+device)/(day)/(run), tmpDat)
#' fitMM
#' L <- getL(fitMM, c("lot1-lot2", "device1-device2"))
#' getDF(fitMM, L)						# method="contain" is Default
#' getDF(fitMM, L, method="res")
#' 
#' getDF(fitMM, L, method="satt")		# takes quite long for this model
#' }

getDF <- function(obj, L, method=c("contain", "residual", "satterthwaite"), ...)
{
	stopifnot(class(obj) == "VCA")
	method <- match.arg(method)
	return(test.fixef(obj, L=L, ddfm=method, onlyDF=TRUE))
}


#' Calculate Variance-Covariance Matrix and Standard Errors of Fixed Effects for an 'VCA' Object.
#' 
#' The variance-covariance matrix of fixed effects for the linear mixed model in 'obj' is calculated.
#' 
#' The variance-covariance matrix of fixed effects for a linear mixed model corresponds to matrix
#' \eqn{(X^{T}V^{-1}X)^{-}}{(X'V"X)`}, where >\eqn{^{T}}{'}< denotes the transpose operator, >\eqn{^{-1}}{"}< 
#' the regular matrix inverse, and >\eqn{^{-}}{`}< the generalized (Moore-Penrose) inverse of a matrix.
#' 
#' @param obj			(VCA) object for which the variance-covariance matrix of
#'                      fixed effects shall be calculated
#' @param quiet			(logical) TRUE = will suppress any warning, which will be issued otherwise 
#' 
#' @return (matrix) corresponding to the variance-covariance matrix of fixed effects
#' 
#' @examples 
#' \dontrun{
#' data(dataEP05A2_1)
#' fit1 <- anovaMM(y~day/(run), dataEP05A2_1)
#' vcov(fit1)
#' 
#' fit2 <- anovaVCA(y~day/run, dataEP05A2_1)
#' vcov(fit2)
#' }

vcovFixed <- function(obj, quiet=FALSE)
{
	Call <- match.call()
	if(!obj$intercept && length(obj$fixed) == 0)
	{
		if(!quiet)
			warning("There is no variance-convariance matrix of fixed effects for this object!")
		return(NA)
	}
	
	X  <- getMat(obj, "X")
	if(is.null(obj$Matrices$Vi))			# MME not yet been solved
	{
		obj  <- solveMME(obj)
		nam  <- as.character(as.list(Call)$obj)
		if(length(nam) == 1 && nam %in% names(as.list(.GlobalEnv)))
		{
			expr <- paste(nam, "<<- obj")		# update object missing MME results
			eval(parse(text=expr))
		}
		else
		{
			if(!quiet)
				message("Some required information missing! Usually solving mixed model equations has to be done as a prerequisite!")
		}
	}
	Vi <- getMat(obj, "Vi")
	VCov <- obj$VarFixed
	if(is.null(VCov))
		VCov <- Solve(t(X) %*% Vi %*% X, quiet=TRUE)
	#VCov <- MPinv(t(X) %*% Vi %*% X)
	rownames(VCov) <- colnames(VCov) <- rownames(obj$FixedEffects)
	
	return(VCov)
}


#' Variance-Covariance Matrix of Fixed Effects as Function of Covariance Parameter Estimates.
#' 
#' This is a helper function for function \code{\link{test.fixef}} approximating degrees of freedom for 
#' linear contrasts of fixed effect parameter estimates.
#' 
#' @param obj			(VCA) object
#' @param x				(numeric) vector of covariance parameter estimates
#' 
#' @return (matrix) corresponding to the variance-covariance matrix of fixed effects

DfSattHelper <- function(obj, x)
{
	stopifnot(class(obj) == "VCA")
	
	rf.ind <- obj$Matrices$rf.ind
	Zi <- obj$Matrices$Zre
	Z <- Vs <- nam <- NULL 
	
	ura <- unique(obj$re.assign[[1]])
	
	for(i in 1:length(x))
	{
		if(i != length(x))
		{
			ind <- which(obj$re.assign[[1]] == i)
			Z <- cbind(Z, as.matrix(Zi[,ind]))
			Vs <- c(Vs, rep(x[i], length(ind)))
			nam <- c(nam, colnames(Zi)[ind])
		}
		else
		{
			Z  <- cbind(Z, diag(obj$Nobs))
			Vs <- c(Vs, rep(x[i], obj$Nobs)) 
		}
	}
	
	G  <- diag(Vs)												# estimated variance-covariance matrix of random effects
	Z  <- Matrix(Z)
	R  <- getMat(obj, "R")
	X  <- getMat(obj, "X")
	V  <- Z %*% G %*% t(Z) + R
	Vi <- Solve(V)
	P  <- Solve(t(X) %*% Vi %*% X, quiet=TRUE)
	P <- as.matrix(P)	
	return(P)
}


#' Perform t-Tests for Linear Contrasts on LS Means.
#' 
#' Perform custom hypothesis tests on Least Squares Means (LS Means) of fixed effect.
#' 
#' This function is similar to function \code{\link{test.fixef}} and represents a convenient way of specifying
#' linear contrasts of LS Means. 
#' 
#' @param obj			(VCA) object
#' @param L				(matrix) specifying one or multiple custom hypothesis tests as linear contrasts of LS Means.
#'                      Which LS Means have to be used is inferred from the column names of matrix \eqn{L}, which need to 
#'                      be in line with the naming of LS Means in function \code{\link{lsmeans}}.
#' @param ddfm			(character) string specifying the method used for computing the denominator
#'                      degrees of freedom of t-tests of LS Means. Available methods are "contain", 
#' 						"residual", and "satterthwaite".
#' @param quiet			(logical) TRUE = will suppress any warning, which will be issued otherwise 
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @seealso \code{\link{test.fixef}}, \code{\link{lsmeans}}
#' 
#' @examples
#' \dontrun{
#' data(dataEP05A2_2)
#' ub.dat <- dataEP05A2_2[-c(11,12,23,32,40,41,42),]
#' fit1 <- anovaMM(y~day/(run), ub.dat)
#' fit2 <- remlMM(y~day/(run), ub.dat)
#' lsm1 <- lsmeans(fit1)
#' lsm2 <- lsmeans(fit2)
#' lsm1
#' lsm2
#' 
#' lc.mat <- getL(fit1, c("day1-day2", "day3-day6"), "lsm")
#' lc.mat[1,c(1,2)] <- c(1,-1)
#' lc.mat[2,c(3,6)] <- c(1,-1)
#' lc.mat
#' test.lsmeans(fit1, lc.mat) 
#' test.lsmeans(fit2, lc.mat)
#' 
#' # fit mixed model from the 'nlme' package
#'  
#' library(nlme)
#' data(Orthodont)
#' fit.lme <- lme(distance~Sex*I(age-11), random=~I(age-11)|Subject, Orthodont) 
#' 
#' # re-organize data for using 'anovaMM'
#' Ortho <- Orthodont
#' Ortho$age2 <- Ortho$age - 11
#' Ortho$Subject <- factor(as.character(Ortho$Subject))
#' 
#' # model without intercept
#' fit.anovaMM <- anovaMM(distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho)
#' fit.remlMM1 <- remlMM( distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho)
#' fit.remlMM2 <- remlMM( distance~Sex+Sex:age2+(Subject)+(Subject):age2-1, Ortho, cov=FALSE)
#' lsm0 <- lsmeans(fit.anovaMM)
#' lsm1 <- lsmeans(fit.remlMM1)
#' lsm2 <- lsmeans(fit.remlMM2)
#' lsm0
#' lsm1
#' lsm2
#' 
#' lc.mat <- matrix(c(1,-1), nrow=1, dimnames=list("int.Male-int.Female", c("SexMale", "SexFemale")))
#' lc.mat
#' test.lsmeans(fit.anovaMM, lc.mat)	
#' test.lsmeans(fit.remlMM1, lc.mat)
#' test.lsmeans(fit.remlMM2, lc.mat)
#' }	

test.lsmeans <- function(obj, L, ddfm=c("contain", "residual", "satterthwaite"), quiet=FALSE)
{
	stopifnot(class(obj) == "VCA")
	lsmm <- lsmMat(obj, NULL, quiet=quiet)
	stopifnot(all(colnames(L) %in% rownames(lsmm)))			# only existing LS Means can be used
	
	lsmm <- lsmm[which(colnames(L) %in% rownames(lsmm)),]
	lsmm <- as.matrix(lsmm)
	
	U <- matrix(nrow=0, ncol=ncol(lsmm))
	colnames(U) <- colnames(lsmm)
	
	for(i in 1:nrow(L))
	{
		U <- rbind(U, matrix(L[i,], nrow=1) %*% lsmm)
	}
	
	if(is.null(obj$VarCov))
		obj$VarCov 	<- vcovVC(obj)
	
	res <- test.fixef(obj, U, ddfm=ddfm, quiet=quiet)
	rownames(res) <- rownames(L)
	return(res)
}


#' Perform t-Tests for Linear Contrasts on Fixed Effects.
#' 
#' This function performs t-Tests for one or multiple linear combinations (contrasts) of estimated 
#' fixed effects.
#' 
#' Here, the same procedure as in \code{SAS PROC MIXED ddfm=satterthwaite} (sat) is implemented. 
#' This implementation was inspired by the code of function 'calcSatterth' of R-package 'lmerTest'. 
#' Thanks to the authors for this nice implementation. \cr
#' Note, that approximated Satterthwaite degrees of freedom might differ from 'lmerTest' and SAS PROC MIXED.
#' Both use the inverse Fisher-information matrix as approximation of the variance-covariance matrix
#' of variance components (covariance parameters). Here, either the exact algorithm for ANOVA-estimators of
#' variance components, described in Searle et. al (1992) p. 176, or the approximation presented in Giesbrecht and 
#' Burns (19985) are used. For balanced designs their will be no differences, usually. 
#' In case of balanced designs, the Satterthwaite approximation is equal to the degrees of freedom of the highest
#' order random term in the model (see examples).
#' 
#' @param obj			(VCA) object
#' @param L				(numeric) vector or matrix, specifying linear combinations of the fixed effects, in the latter case,
#'                      each line represents a disctinct linear contrast
#' @param ddfm			(character) string specifying the method used for computing the denominator
#'                      degrees of freedom for tests of fixed effects or LS Means. Available methods are
#' 						"contain", "residual", and "satterthwaite".
#' @param method.grad	(character) string specifying the method to be used for approximating the gradient
#'                      of the variance-covariance matrix of fixed effects at the estimated covariance parameter
#'                      estimates (see function 'grad' (numDeriv) for details)
#' @param tol			(numeric) value specifying the numeric tolerance for testing equality to zero
#' @param quiet			(logical) TRUE = suppress warning messages, e.g. for non-estimable contrasts
#' @param opt			(logical) TRUE = tries to optimize computation time by avoiding unnecessary computations
#'                      for balanced datasets (see details). 
#' @param onlyDF		(logical) TRUE = only the specified type of degrees of freedom are determined without carrying out
#' 						the actual hypothesis test(s)
#' @param ...			further parameters (for internal use actually)
#'  
#' @return (numeric) vector or matrix with 4 elements/columns corresponding to "Estimate", "t Value", "DF", and
#'         "Pr > |t|".
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com} inspired by authors of R-package 'lmerTest'
#' 
#' @references 
#' 
#' Searle, S.R, Casella, G., McCulloch, C.E. (1992), Variance Components, Wiley New York	
#' 
#' Giesbrecht, F.G. and Burns, J.C. (1985), Two-Stage Analysis Based on a Mixed Model: Large-Sample
#' Asymptotic Theory and Small-Sample Simulation Results, Biometrics 41, p. 477-486
#' 
#' SAS Help and Documentation PROC MIXED (MODEL-statement, Option 'ddfm'), SAS Institute Inc., Cary, NC, USA
#' 
#' @seealso \code{\link{test.lsmeans}}, \code{\link{getL}}
#' 
#' @examples
#' \dontrun{
#' data(dataEP05A2_2)
#' ub.dat <- dataEP05A2_2[-c(11,12,23,32,40,41,42),]
#' fit1 <- anovaMM(y~day/(run), ub.dat)
#' fit2 <- remlMM(y~day/(run), ub.dat)
#' fe1 <- fixef(fit1)
#' fe1
#' fe2 <- fixef(fit2)
#' fe2
#' lc.mat <- getL( fit1, c("day1-day2", "day3-day6"))
#' lc.mat
#' test.fixef(fit1, lc.mat, ddfm="satt") 
#' test.fixef(fit2, lc.mat, ddfm="satt")
#' 
#' # some inferential statistics about fixed effects estimates
#' L <- diag(nrow(fe1))
#' rownames(L) <- colnames(L) <- rownames(fe1)
#' test.fixef(fit1, L)
#' test.fixef(fit2, L)
#' 
#' # using different "residual" method determining DFs
#' test.fixef(fit1, L, ddfm="res")
#' test.fixef(fit2, L, ddfm="res")  
#' 
#' # having 'opt=TRUE' is a good idea to save time 
#' # (in case of balanced designs)
#' data(VCAdata1)
#' datS3 <- VCAdata1[VCAdata1$sample==3,]
#' fit3 <- anovaMM(y~(lot+device)/(day)/run, datS3)
#' fit4 <- remlMM(y~(lot+device)/(day)/run, datS3)  
#' fit3$VarCov <- vcovVC(fit3)
#' fe3 <- fixef(fit3)
#' fe4 <- fixef(fit4)
#' L <- diag(nrow(fe3))
#' rownames(L) <- colnames(L) <- rownames(fe3)
#' system.time(tst1 <- test.fixef(fit3, L))
#' system.time(tst2 <- test.fixef(fit3, L, opt=FALSE))
#' system.time(tst3 <- test.fixef(fit4, L, opt=FALSE))
#' tst1
#' tst2
#' tst3
#' }

test.fixef <- function(	obj, L, ddfm=c("contain", "residual", "satterthwaite"),
		method.grad="simple", tol=1e-12, quiet=FALSE, opt=TRUE,
		onlyDF=FALSE, ...)
{
	stopifnot(class(obj) == "VCA")
	ddfm <- match.arg(ddfm)
	
	if(obj$EstMethod == "REML" && ddfm == "contain")
	{
		ddfm <- "satterthwaite"
		if(!quiet)
			message("Note: 'ddfm' set to \"satterthwaite\", currently option \"contain\" does not work for REML-estimation!")
	}
	
	args <- list(...)
	lsmeans <- args$lsmeans						# this function is called when LS Means with complex output is requested
	if(is.null(lsmeans))
		lsmeans <- FALSE
	
	b <- fixef(obj)[,"Estimate", drop=F]
	
	if(is.null(dim(L)))
		L <- matrix(L, nrow=1)
	
	stopifnot(ncol(L) == nrow(b))
	
	if(nrow(L) > 1)
	{
		res <- matrix(ncol=5, nrow=nrow(L))
		colnames(res) <- c("Estimate", "DF", "SE", "t Value", "Pr > |t|")
		
		for(i in 1:nrow(L))
			res[i,] <- test.fixef(obj=obj, L=L[i,,drop=F], ddfm=ddfm, method.grad=method.grad, quiet=quiet, opt=opt, onlyDF=onlyDF, lsmeans=lsmeans)
		
		rownames(res) <- rownames(L)
		if(onlyDF)
			return(res[,"DF"])
		attr(res, "ddfm") <- ddfm
		return(res)
	}	
	else
	{
		ind <- which(L != 0)						# test whether fixef(obj, "complex") was called
		
		if(length(ind) == 1 && any(abs(b[ind,]) < tol) && !lsmeans)
		{
			if(!quiet)
				warning("Linear contrast'", L ,"' not estimable!")
			res <- rep(NA, 5)
			names(res) <- c("Estimate", "DF", "SE", "t Value", "Pr > |t|")
			return(res)
		}		
		est		<- as.numeric(L %*% b)
		sgn		<- sign(est)
		X		<- getMat(obj, "X")
		P    	<- vcov(obj)
		lPl  	<- L %*% P %*% t(L)
		lPli 	<- solve(lPl)
		r       <- as.numeric(rankMatrix(lPli))
		SVD 	<- eigen(lPl)		
		items   <- list(SVD=SVD, r=r, b=b)	
		
		se <- L %*% P %*% t(L)						# compute standard error of linear contrast
		se <- sqrt(diag(se))
		
		DF		<- getDDFM(obj, L, ddfm, tol=tol, method.grad=method.grad, opt=opt, items=items)
		if(onlyDF)
		{
			return(c(NA, DF, NA, NA, NA))
		}
		t.stat 	<- as.numeric(sqrt((t(L %*% b) %*% lPli %*% (L %*% b))/r))
		
		res <- matrix(c(est, DF, se, sgn*t.stat, 2*pt(abs(t.stat), df=DF, lower.tail=FALSE)), nrow=1,
						dimnames=list(NULL, c("Estimate", "DF", "SE", "t Value", "Pr > |t|"))) 
		
		attr(res, "ddfm") <- ddfm
		
		return( res )
	}	
}


#' Construct Linear Contrast Matrix for Hypothesis Tests.
#' 
#' Function constructs coefficient/contrast matrices from a string-representation of linear hypotheses.
#' 
#' Function constructs matrices expressing custom linear hypotheses of fixed effects or
#' LS Means. The user has to specify a string denoting this contrast which is then 
#' transformed into a coefficient/contrast matrix. This string may contain names of fixed effects
#' belonging to same same fixed term, numeric coefficients and mathematical operators "+"
#' and "-" (see examples).
#' 
#' @param obj			(VCA) object
#' @param s				(character) string or vector of strings, denoting one or multiple
#'                      linear contrasts
#' @param what			(character) string specifying whether to construct contrast matrices
#'                      of fixed effects ("fixed") or LS Means ("lsmeans"), abbreviations are allowed.
#' 
#' @return (matrix) representing one linear hypothesis of fixed effects or LS Means per row 
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @examples
#' \dontrun{
#' data(dataEP05A2_2)
#' fit <- anovaMM(y~day/(run), dataEP05A2_2)
#' L <- getL(fit, c("day1-day2", "day5-day10"), what="fixef")
#' L
#' test.fixef(fit, L=L)
#' 
#' # another custom hypothesis
#' L2 <- getL(fit, "0.25*day1+0.25*day2+0.5*day3-0.5*day4-0.5*day5")
#' L2
#' 
#' # more complex model
#' data(VCAdata1)
#' dataS2 <- VCAdata1[VCAdata1$sample==2,]
#' fit.S2 <- anovaMM(y~(lot+device)/day/(run), dataS2)
#' L3 <- getL(fit.S2, c("lot1-lot2", "lot1:device3:day19-lot1:device3:day20", 
#' 						"lot1:device1:day1-lot1:device1:day2"))
#' L3
#' test.fixef(fit.S2, L3)
#' }

getL <- function(obj, s, what=c("fixef", "lsmeans"))
{
	what <- match.arg(what)
	nwhat <- c(fixef="fixed effects", lsmeans="LS Means")
	
	if(what == "fixef")
		b <- fixef(obj, quiet=TRUE)
	else
		b <- lsmeans(obj, quiet=TRUE)
	n <- rownames(b)
	
	if(length(s) > 1)
	{
		L <- try(t(sapply(s, function(x) getL(obj, x, what=what))), silent=TRUE)
		if(class(L) == "try-error")
			stop(L[1])
		colnames(L) <- n
		return(L)
	}	
	contr <- rep(0, nrow(b))
	cname <- names(s)
	s <- gsub("\\*", "", s)
	
	splt <- sapply(unlist(strsplit(s, "\\+")), strsplit, "\\-")
	
	fac <- NULL
	sgn <- NULL
	
	for(i in 1:length(splt))
	{
		if(i == 1)
		{
			if(splt[[1]][1] == "")
				sgn <- -1
			else
				sgn <- 1 
		}
		else
			sgn <- 1									# 1st element always "+" since it was firstly used as split-char
		
		sgn <- c(sgn, rep(-1, length(splt[[i]])-1))
		
		for(j in 1:length(splt[[i]]))
		{
			tmp <- regexpr("^[[:digit:]]*\\.?[[:digit:]]*", splt[[i]][j])
			
			if(tmp > -1 && attr(tmp, "match.length") > 0)				# there is a factor at the beginning of the string
			{
				tmp.fac <- substr(splt[[i]][j], 1, attr(tmp, "match.length"))
				fac <- c(fac, as.numeric(tmp.fac) * sgn[j])	
				splt[[i]][j] <- sub(tmp.fac, "", splt[[i]][j])			# remove factor sub-string
			}
			else
				fac <- c(fac, 1 * sgn[j])
		}
	}
	
	splt <- unlist(splt)
	names(splt) <- NULL
	
	if(!all(splt %in% n))
		stop("\nError: There are terms which do not belong to the '",nwhat[what],"'!")
	
	res <- sapply(n, regexpr, splt)
	rownames(res) <- splt
	cnl <- nchar(colnames(res))	
	
	tms <- apply(res, 1, function(x){
							ind <- which(x > 0)
							len <- cnl[ind]
							return(ind[which(len == max(len))])
						 })
	contr[tms] <- fac	
	return(matrix(contr, nrow=1, dimnames=list(cname, n)))
}

#' Degrees of Freedom for Testing Linear Contrasts of Fixed Effects and Least Square Means.
#' 
#' There are three methods implemented, which are all available in SAS PROC MIXED, namely 
#' "contain", "residual", and "satterthwaite" approximations. See the documentation of SAS
#' PROC MIXED for details on this topic.
#' 
#' The implementation of the Satterthwaite approximation was inspired by the code of function 
#' 'calcSatterth' of R-package 'lmerTest'.
#' 
#' @param obj			(VCA) object
#' @param L				(numeric) vector specifying the linear combination of the fixed effect or
#'                      LS Means
#' @param ddfm			(character) string specifying the method used for computing the denominator
#'                      degrees of freedom for tests of fixed effects or LS Means. Available methods are
#' 						"contain", "residual", and "satterthwaite".
#' @param tol			(numeric) value specifying the numeric tolerance for testing equality to zero
#' @param method.grad	(character) string specifying the method to be used for approximating the gradient
#'                      of the variance-covariance matrix of fixed effects at the estimated covariance parameter
#'                      estimates (see function 'grad' (numDeriv) for details)
#' @param opt			(logical) TRUE = tries to optimize computation time by avoiding unnecessary computations
#'                      for balanced datasets (see \code{\link{test.fixef}}). 
#' @param items			(list) of pre-computed values
#' 
#' @return (numeric) vector with the specified type of degrees of freedom
#' 
#' @author Andre Schuetzenmeister \email{[email protected]@roche.com}
#' 
#' @seealso \code{\link{test.fixef}}

getDDFM <- function(obj, L, ddfm=c("contain", "residual", "satterthwaite"), tol=1e-12, method.grad="simple", opt=TRUE, items=NULL)
{
	stopifnot(class(obj) == "VCA")
	stopifnot(!is.null(colnames(L)))
	ddfm <- match.arg(ddfm)
	
	if(ddfm == "residual")
	{
		return(obj$Nobs - rankMatrix(getMat(obj, "X")))			# as described in documentation of SAS PROC MIXED
	}
	else if(ddfm == "contain")
	{		
		cn <- colnames(L)										# fe or LS Means names
		cn <- cn[which(L[1,] != 0)]
		
		if(length(cn) == 1 && cn == "int")						# handle intercept fixed effect
		{
			return(DF <- min(obj$aov.org[obj$random, "DF"]))
		}
		fe <- obj$fixed
		tmp <- sapply(fe, function(x) gregexpr(x, cn))			# can names of fixed terms be found in column names of L?
		
		if(class(tmp) == "list")								# there was only a single non-zero column in L --> list returned
			fe <- sapply(tmp, function(x) x != -1)
		else													# mulitple non-zero columns in L
		{
			fe <- apply(tmp, 2, function(y) any(y==1))
		}		
		
		if(any(fe))
		{
			fe <- fe[which(fe)]
			fe <- names(fe)
			rn <- obj$random
			rn <- sapply(rn, function(x) all(sapply(fe, grepl, x)))
			if(any(rn))
			{
				DF <- min(obj$aov.org[names(rn), "DF"])
				return(DF)
			}
			else
			{
				tmpZ <-  getMat(obj, "Z")
				
				if(!is.null(tmpZ))
					return(obj$Nobs - rankMatrix(cbind(as.matrix(getMat(obj, "X")), as.matrix(tmpZ))))
				else
					return(obj$Nobs - rankMatrix(as.matrix(getMat(obj, "X"))))
			}
		}
		else
		{
			return(obj$Nobs - rankMatrix(cbind(as.matrix(getMat(obj, "X")), as.matrix(getMat(obj, "Z")))))	# no random term including the fixed effects term --> residual DFs
		}
	}
	else if(ddfm == "satterthwaite")
	{
		stopifnot(class(obj) == "VCA")		
		if(is.null(dim(L)))
			L <- matrix(L, nrow=1)		
		
		if(obj$balanced == "balanced" && opt && obj$EstMethod != "REML" && FALSE)
		{
			return(obj$aov.tab[2,"DF"])			
		}
		r       <- items$r
		SVD 	<- items$SVD
		nu.m 	<- NULL											# see SAS-help for SAS PROC MIXED -> model /ddfm=sat
		A <- obj$VarCov