R/MMFramework.R

Defines functions parserOutputSummary finalModel modelFormula startModel

Documented in finalModel modelFormula parserOutputSummary startModel

## Copyright © 2012-2014 EMBL - European Bioinformatics Institute
##
## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License.
## You may obtain a copy of the License at
##
##     http://www.apache.org/licenses/LICENSE-2.0
##
## Unless required by applicable law or agreed to in writing, software
## distributed under the License is distributed on an "AS IS" BASIS,
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
## See the License for the specific language governing permissions and
## limitations under the License.
##------------------------------------------------------------------------------
## MMFramework.R contains startModel, finalModel, modelFormula
## and parserOutputSummary functions
##------------------------------------------------------------------------------
## startModel creates start model and modifies it after testing of different
## hypothesis (the model effects).
## The model effects are:
## - batch effect (random effect significance),
## - variance effect (TRUE if residual variances for genotype groups are
## homogeneous and FALSE if they are heterogeneous),

## Fixed effects:
## - interaction effect (genotype by sex interaction significance),
## - sex effect (sex significance),
## - weigth effect (weigth significance).

## If user would like to assign other TRUE/FALSE values to the effects of
## the model then he or she has to define keepList argument
## which is vector of TRUE/FALSE values.

## If user has defined fixed effects (keepList argument) then function
## prints out calculated and user defined effects
## (only when outputMessages argument is set to TRUE),
## checks user defined effects for consistency
## (for instance, if there are no "Weight" column in the dataset then weight
## effect can't be assigned to TRUE, etc.)
## and modifies start model according to user defined effects.

## As the result PhenTestResult object that contains calculated or user
## defined model effects and MM start model is created.
startModel <- function(phenList,
											 depVariable,
											 equation = "withWeight",
											 outputMessages = TRUE,
											 pThreshold = 0.05,
											 keepList = NULL,
											 modelWeight = NULL,
											 threshold = 10 ^ -18,
											 check = 1)
{
	x <- getDataset(phenList)
	mw = checkWeights(
		w = modelWeight,
		n = nrow(x),
		threshold = threshold,
		date      =  x$Batch,
		check = check
	)
	WindowingActive = !(is.null(modelWeight) || is.null(mw$w))
	if (WindowingActive) {
		message('heterogeneous residual variances for genotype groups is replaced by model weights')
		x              = x[mw$wInd,]
		x$ModelWeight  = mw$w
		FixV  = CombV = nlme::varFixed( ~ 1 / ModelWeight)
		# CombV = varComb(varFixed( ~ 1 / ModelWeight), varIdent(form =  ~ 1 |
		# 																											 	Genotype))
	} else{
		FixV  = NULL
		CombV = nlme::varIdent(form =  ~ 1 | Genotype)
	}

	numberofsexes <- length(levels(x$Sex))
	mean_list <- NULL
	if (!is.null(keepList)) {
		## User's values for effects
		user_keep_weight <- keepList[3]
		user_keep_sex <- keepList[4]
		user_keep_interaction <- keepList[5]
		user_keep_batch <- keepList[1]
		user_keep_equalvar <- keepList[2]
		if (!('Weight' %in% colnames(x)) && user_keep_weight) {
			if (outputMessages)
				message("Warning:\nWeight column is missed in the dataset. 'keepWeight' is set to FALSE.")
			user_keep_weight <- FALSE
		}
		if (!('Batch' %in% colnames(x)) && user_keep_batch) {
			if (outputMessages)
				message("Warning:\nBatch column is missed in the dataset. 'keepBatch' is set to FALSE.")
			user_keep_batch <- FALSE
		}
	}
	# Averages for percentage changes - is the ratio of the genotype effect for a sex relative to
	# the wildtype signal for that variable for that sex - calculation
	# WT <- subset(x, x$Genotype == refGenotype(phenList))
	#mean_all <- mean(WT[,c(depVariable)],na.rm=TRUE)
	mean_all <- mean(x[, c(depVariable)], na.rm = TRUE)
	mean_list <- c(mean_all)
	if (numberofsexes == 2) {
		#WT_f <- subset(WT,WT$Sex=="Female")
		#WT_m <- subset(WT,WT$Sex=="Male")
		#mean_f <- mean(WT_f[,c(depVariable)],na.rm=TRUE)
		#mean_m <- mean(WT_m[,c(depVariable)],na.rm=TRUE)
		all_f <- subset(x, x$Sex == "Female")
		all_m <- subset(x, x$Sex == "Male")
		mean_f <- mean(all_f[, c(depVariable)], na.rm = TRUE)
		mean_m <- mean(all_m[, c(depVariable)], na.rm = TRUE)
		mean_list <- c(mean_all, mean_f, mean_m)
	}

	# end of percentage change calculations

	## Start model formula: homogenous residual variance,
	## genotype and sex interaction included
	model.formula  <-
		modelFormula(equation, numberofsexes, depVariable)


	#START OF tryCatch
	finalResult <- tryCatch({
		if (batchIn(phenList)) {
			## GLS fit of model formula (no random effects)
			## Model 1A (model_withoutbatch)
			model_withoutbatch <- do.call("gls",
																		args = list(
																			model = model.formula,
																			data  = x,
																			weights = FixV,
																			na.action = "na.omit"
																		))
			## MM fit of model formula (with random effects)
			## Model 1 (model_MM)
			model_MM <-
				do.call(
					"lme",
					args = list(
						fixed  = model.formula,
						random =  ~ 1 |
							Batch,
						data = x,
						weights =  FixV,
						na.action = "na.omit",
						method = "REML"
					)
				)
			## Test: the random effects associated with batch intercepts can be
			## ommited from model
			## Hypothesis 1
			## Null Hypothesis: variance of batch = 0
			## Alternative Hypothesis: variance of batch > 0
			## For the division by 2 explanations see p.80 of "Linear Mixed Models"

			p.value.batch <-
				(anova(model_MM, model_withoutbatch)$"p-value"[2]) / 2
			## The result of the test for Hypothesis 1 will help to select
			## the structure for random effects
			keep_batch <- p.value.batch < pThreshold

			## MM fit of model formula with heterogeneous residual variances for
			## genotype groups
			## Model 1 assumes homogeneous residual variances
			## Model 2 with heterogeneous residual variances
			model_hetvariance <-
				do.call(
					"lme",
					args = list(
						fixed = model.formula,
						random =  ~ 1 |
							Batch,
						data    = x,
						weights = CombV,
						na.action = "na.omit",
						method = "REML"
					)
				)


			## Test: the variance of the residuals is the same (homogeneous)
			## for all genotype groups
			## Hypothesis 2
			## Null Hypothesis: all residual variances are equal
			## Alternative Hypothesis: the residue variance is not equal
			if (WindowingActive) {
				p.value.variance = NA
				keep_equalvar = FALSE
			} else{
				p.value.variance <-
					(anova(model_MM, model_hetvariance)$"p-value"[2])
				## The result of the test for Hypothesis 2 will help to select a
				## covariance structure for the residuals
				keep_equalvar <- p.value.variance > pThreshold
			}
		} else {
			## No Batch effects
			keep_batch <- FALSE

			## Model 1A (model_withoutbatch)
			model_MM <- do.call("gls",
													args = list(
														model = model.formula,
														data = x,
														weights = FixV,
														na.action = "na.omit"
													))

			model_withoutbatch <- model_MM

			## MM fit of model formula with heterogeneous residual variances for
			## genotype groups
			## Model 1 assumes homogeneous residual variances
			## Model 2 with heterogeneous residual variances
			model_hetvariance <-
				do.call(
					"gls",
					args = list(
						model = model.formula,
						data  = x,
						weights = CombV,
						na.action = "na.omit"
					)
				)

			## Test: the variance of the residuals is the same (homogeneous)
			## for all genotype groups
			## Hypothesis 2
			## Null Hypothesis: all residual variances are equal
			## Alternative Hypothesis: the residue variance is not equal
			if (WindowingActive) {
				p.value.variance = NA
				keep_equalvar = FALSE
			} else{
				p.value.variance <-
					(anova(model_MM, model_hetvariance)$"p-value"[2])
				## The result of the test for Hypothesis 2 will help to select a
				## covariance structure for the residuals
				keep_equalvar <- p.value.variance > pThreshold
			}
		}


		## Model fit is selected according to test results
		if (keep_batch && keep_equalvar) {
			## Model 1
			model <- model_MM
		} else if (keep_batch && !keep_equalvar) {
			## Model 2
			model <- model_hetvariance
		} else if (!keep_batch && keep_equalvar) {
			## Model 1A
			model = model_withoutbatch
		} else if (!keep_batch && !keep_equalvar) {
			## Modify model 1A to heterogeneous residual variances
			model <- do.call(
				"gls",
				args = list(
					model = model.formula,
					data = x,
					weights = CombV,
					na.action = "na.omit"
				)
			)
		}


		## Tests for significance of fixed effects using TypeI F-test from anova
		## functionality by using selected model
		anova_results <-
			anova(model, type = "marginal")$"p-value" < pThreshold

		if (numberofsexes == 2) {
			## Result of the test for sex significance (fixed effect 1.)
			keep_sex <- anova_results[3]
			## Eq.2
			if (equation == "withWeight") {
				## Result of the test for weight significance  (fixed effect 3.)

				keep_weight <- anova_results[4]
				## Result of the test for genotype by sex interaction
				## significance (fixed effect 2.)
				keep_interaction <- anova_results[5]

				## Technical results needed for the output
				## Interaction test results are kept for the output
				interactionTest <-
					anova(model, type = "marginal")$"p-value"[5]

			} else{
				## Result of the test for weight significance  (fixed effect 3.)
				## It's FALSE since here the equation 1 is used - without weight
				## effect
				keep_weight <- FALSE
				## Result of the test for genotype by sex interaction
				## significance (fixed effect 2.)
				keep_interaction <- anova_results[4]
				## Interaction test results are kept for the output
				interactionTest <-
					anova(model, type = "marginal")$"p-value"[4]
			}
		} else {
			keep_sex <- FALSE
			keep_interaction <- FALSE
			interactionTest <- NA
			if (equation == "withWeight")
				keep_weight <- anova_results[3]
			else
				keep_weight <- FALSE
		}

		if (!keep_weight && equation == "withWeight") {
			equation = "withoutWeight"
			if (outputMessages)
				message(
					paste(
						"Since weight effect is not significant the equation ",
						"'withoutWeight' should be used instead.",
						sep = ""
					)
				)
		}

		if (outputMessages)
			message(paste("Information:\nEquation: '", equation, "'.\n", sep =
											""))

		if (outputMessages)
			message(
				paste(
					"Information:\nCalculated values for model effects are: ",
					"keepBatch=",
					keep_batch,
					", keepEqualVariance=",
					keep_equalvar,
					", keepWeight=",
					keep_weight,
					", keepSex=",
					keep_sex,
					", keepInteraction=",
					keep_interaction,
					".\n",
					sep = ""
				)
			)

		## Results for user defined model effects values
		if (!is.null(keepList)) {
			if (outputMessages)
				message(
					paste(
						"Information:\nUser's values for model effects are: ",
						"keepBatch=",
						user_keep_batch,
						", keepEqualVariance=",
						user_keep_equalvar,
						", keepWeight=",
						user_keep_weight,
						", keepSex=",
						user_keep_sex,
						", keepInteraction=",
						user_keep_interaction,
						".\n",
						sep = ""
					)
				)
			## Model fit is selected according to user defined model effects
			if (user_keep_batch && user_keep_equalvar) {
				## Model 1
				model <- model_MM
			} else if (user_keep_batch &&
								 !user_keep_equalvar) {
				## Model 2
				model <- model_hetvariance
			} else if (!user_keep_batch &&
								 user_keep_equalvar) {
				## Model 1A
				model <- model_withoutbatch
			} else if (!user_keep_batch &&
								 !user_keep_equalvar) {
				## Modify model 1A to heterogeneous residual variances
				model <-
					do.call(
						"gls",
						args = list(
							model = model.formula,
							data = x,
							weights = CombV,
							na.action = "na.omit"
						)
					)
			}

			if (numberofsexes == 2) {
				if (equation == "withWeight") {
					interactionTest <- anova(model, type = "marginal")$"p-value"[5]
				}        else{
					interactionTest <- anova(model, type = "marginal")$"p-value"[4]

				}
			}  else {
				interactionTest <- NA
			}

			compList <-
				(
					keepList == c(
						keep_batch,
						keep_equalvar,
						keep_weight,
						keep_sex,
						keep_interaction
					)
				)

			if (length(compList[compList == FALSE]) > 0 &&
					outputMessages)
				message(
					"Warning:\nCalculated values differ from user defined values for model effects.\n"
				)

			keep_weight <- user_keep_weight
			keep_sex <- user_keep_sex
			keep_interaction <- user_keep_interaction
			keep_batch <- user_keep_batch
			keep_equalvar <- user_keep_equalvar

		}

		#MM modelling output
		linearRegressionOutput <- list(
			model.output = model,
			equation = equation,
			model.effect.batch = keep_batch,
			model.effect.variance = keep_equalvar,
			model.effect.interaction = keep_interaction,
			model.output.interaction = interactionTest,
			model.effect.sex = keep_sex,
			model.effect.weight = keep_weight,
			numberSexes = numberofsexes,
			pThreshold = pThreshold,
			model.formula.genotype = model.formula,
			model.output.averageRefGenotype = mean_list
		)

		finalResult <- new(
			"PhenTestResult",
			analysedDataset = x,
			depVariable = depVariable,
			refGenotype = refGenotype(phenList),
			testGenotype = testGenotype(phenList),
			method = "MM",
			transformationRequired = FALSE,
			lambdaValue = integer(0),
			scaleShift = integer(0),
			analysisResults = linearRegressionOutput,
			modelWeight = mw,
			phenList = phenList,
			equation = equation,
			outputMessages = outputMessages,
			pThreshold = pThreshold#,
			#keepList = keepList
		)

	},
	#END OF tryCatch
	#error = function(e) {}
	error = function(error_mes) {
		message(error_mes)
		finalResult <- NULL
		if (equation == "withWeight")
			stop_message <-
			paste(
				"Error:\nCan't fit the model ",
				format(model.formula),
				". Try MM with equation 'withoutWeight'. ",
				"Another option is jitter.\n",
				sep = ""
			)
		else
			stop_message <-
			paste(
				"Error:\nCan't fit the model ",
				format(model.formula),
				". Try to add jitter or RR plus method.\n",
				sep = ""
			)

		if (outputMessages) {
			message(stop_message)
			opt <- options(show.error.messages = FALSE)
			on.exit(options(opt))
			stop()
		}  else {
			stop(stop_message)
		}
	},
	warning = function(w) {
		message(w)
	})
	# agg = c(as.list(environment()), list())
	# save(agg, file = file.path(getwd(), 'superDebugAnalysis.Rdata'))
	return(finalResult)

}

##------------------------------------------------------------------------------
## Creates formula for the start model based on equation and number of sexes
## in the data
modelFormula <- function(equation, numberofsexes, depVariable)
{
	model.formula <- switch(equation,
													## Eq.2
													withWeight = {
														## Fixed effects: 1) Genotype 2) Sex 3) Genotype by Sex
														## interaction 4) Weight
														if (numberofsexes == 2) {
															as.formula(paste(
																depVariable,
																"~",
																paste("Genotype",
																			"Sex", "Genotype*Sex", "Weight", sep = "+")
															))
														} else{
															as.formula(paste(depVariable, "~", paste("Genotype",
																																			 "Weight",  sep = "+")))
														}
													},
													## Eq.1
													withoutWeight = {
														## Fixed effects: 1) Genotype 2) Sex 3) Genotype by Sex
														## interaction
														if (numberofsexes == 2) {
															as.formula(paste(
																depVariable,
																"~",
																paste("Genotype", "Sex", "Genotype*Sex", sep = "+")
															))
														} else{
															as.formula(paste(depVariable, "~",
																							 paste("Genotype",  sep = "+")))
														}
													})
	return(model.formula)
}

##------------------------------------------------------------------------------
## Works with PhenTestResult object created by testDataset function.
## Builds final model based on the significance of different model effects,
## depVariable and equation stored in phenTestResult object (see testDataset.R).
finalModel <-
	function(phenTestResult,
					 outputMessages = TRUE,
					 modelWeight = NULL)
	{
		## Checks and stop messages
		stop_message <- ""

		## Check PhenTestResult object
		if (is(phenTestResult, "PhenTestResult")) {
			finalResult <- phenTestResult
			linearRegressionOutput <- analysisResults(phenTestResult)
			x <- analysedDataset(finalResult)
			#####
			mw = x$ModelWeight
			if (!is.null(modelWeight) && !is.null(mw)) {
				FixV  = CombV = nlme::varFixed( ~ 1 / ModelWeight)
				# CombV = varComb(varFixed(~ 1 / ModelWeight) ,varIdent(form =  ~ 1 |
				# 												 	Genotype)
				#								)
			} else{
				FixV  = NULL
				CombV = nlme::varIdent(form =  ~ 1 | Genotype)
			}
			#####
			depVariable <- getVariable(finalResult)
			keep_weight <- linearRegressionOutput$model.effect.weight
			keep_sex <- linearRegressionOutput$model.effect.sex
			keep_interaction <-
				linearRegressionOutput$model.effect.interaction
			keep_batch <- linearRegressionOutput$model.effect.batch
			keep_equalvar <-
				linearRegressionOutput$model.effect.variance

			## Stop function if there are no datasets to work with
			if (is.null(x))
				stop_message <-
				"Error:\nPlease create a PhenList object first and run function 'testDataset'.\n"

			## Stop function if there are no enough input parameters
			if (is.null(linearRegressionOutput$equation) ||
					is.null(depVariable) || is.null(keep_batch)
					|| is.null(keep_equalvar)
					|| is.null(keep_sex) || is.null(keep_interaction))
				stop_message <-
				"Error:\nPlease run function 'testDataset' first.\n"
		}  else{
			stop_message <-
				"Error:\nPlease create a PhenTestResult object first.\n"
		}



		if (nchar(stop_message) > 0) {
			if (outputMessages) {
				message(stop_message)
				opt <- options(show.error.messages = FALSE)
				on.exit(options(opt))
				stop()
			}    else {
				stop(stop_message)
			}
		}

		## END Checks and stop messages


		## Build final null model
		## Goal:  to test fixed effects of the model and based on the output
		## build the final null model formula for later
		## testing - as a null model it automatically excludes genotype and
		## interaction term.
		## The anova function tests the fixed effects associated by treatment
		## with a null hypothesis that the regression
		## coefficients are equal to zero  and an alternative hypothesis that the
		## regression coefficient are not equal to zero.
		## If the p-values of these tests are less than 0.05 we reject the null
		## and accept the alternative that the are
		## significant components of the model and should be included.
		## If no terms are significant a model can be build with just an
		## intercept element this is specified as
		## "model.formula <- as.formula(paste(depVariable, "~", "1"))"

		## Null model: genotype is not significant
		model_null.formula <- switch(
			linearRegressionOutput$equation,
			withWeight = {
				## Eq.2
				if (linearRegressionOutput$numberSexes == 2) {
					if (!keep_sex) {
						as.formula(paste(depVariable, "~", "Weight"))

					} else{
						as.formula(paste(depVariable, "~",
														 paste("Sex", "Weight", sep = "+")))
					}
				} else{
					as.formula(paste(depVariable, "~", "Weight"))
				}
			},
			withoutWeight = {
				## Eq.1
				if (linearRegressionOutput$numberSexes == 2) {
					if (!keep_sex && !keep_interaction) {
						as.formula(paste(depVariable, "~", "1"))
					} else{
						as.formula(paste(depVariable, "~", "Sex"))
					}
				} else{
					as.formula(paste(depVariable, "~", "1"))
				}
			}
		)

		## Alternative model: genotype is significant
		model_genotype.formula <-
			switch(
				linearRegressionOutput$equation,
				withWeight = {
					## Eq.2
					if (linearRegressionOutput$numberSexes == 2) {
						if ((keep_sex && keep_weight && keep_interaction) |
								(!keep_sex &&
								 keep_weight && keep_interaction)) {
							as.formula(paste(
								depVariable,
								"~",
								paste("Sex", "Genotype:Sex", "Weight", sep = "+")
							))

						} else if (keep_sex &&
											 keep_weight && !keep_interaction) {
							as.formula(paste(
								depVariable,
								"~",
								paste("Genotype", "Sex", "Weight", sep = "+")
							))

						} else if (!keep_sex &&
											 keep_weight && !keep_interaction) {
							as.formula(paste(depVariable, "~",
															 paste("Genotype", "Weight", sep = "+")))
						}

					} else{
						as.formula(paste(depVariable, "~",
														 paste("Genotype", "Weight", sep = "+")))
					}
				},
				withoutWeight = {
					## Eq.1
					if (linearRegressionOutput$numberSexes == 2) {
						if (!keep_sex  && !keep_interaction) {
							as.formula(paste(depVariable, "~", "Genotype"))
						} else if ((keep_sex && keep_interaction) |
											 (!keep_sex && keep_interaction)) {
							as.formula(paste(
								depVariable,
								"~",
								paste("Sex",  "Genotype:Sex", sep = "+")
							))
						} else if (keep_sex && !keep_interaction) {
							as.formula(paste(depVariable, "~",
															 paste("Genotype", "Sex", sep = "+")))
						}
					} else{
						as.formula(paste(depVariable, "~", paste("Genotype")))
					}
				}
			)

		finalResult <- tryCatch({
			## Test: genotype groups association with dependent variable
			## Null Hypothesis: genotypes are not associated with dependent variable
			## Alternative Hypothesis: genotypes are associated with dependent
			## variable
			if (keep_batch && keep_equalvar) {
				model_genotype <-
					do.call(
						"lme",
						args = list(
							fixed = model_genotype.formula,
							random =  ~ 1 |
								Batch,
							data = x,
							weights = FixV,
							na.action = "na.omit",
							method = "ML"
						)
					)

				model_null <-
					do.call(
						"lme",
						args = list(
							fixed = model_null.formula,
							data  = x,
							weights =  FixV,
							random =  ~ 1 |
								Batch,
							na.action = "na.omit",
							method = "ML"
						)
					)
				p.value <-
					(anova(model_genotype, model_null)$"p-value"[2])

			}
			if (keep_batch && !keep_equalvar) {
				model_genotype <- do.call(
					"lme",
					args = list(
						fixed = model_genotype.formula,
						random =  ~ 1 |
							Batch,
						data = x,
						weights = CombV,
						na.action = "na.omit",
						method = "ML"
					)
				)

				model_null <-
					do.call(
						"lme",
						args = list(
							fixed  = model_null.formula,
							random =  ~ 1 |
								Batch,
							data = x,
							weights = CombV,
							na.action = "na.omit",
							method = "ML"
						)
					)

				p.value <-
					(anova(model_genotype, model_null)$"p-value"[2])
			} else if (!keep_batch && !keep_equalvar) {
				model_genotype <- do.call(
					"gls",
					args = list(
						model = model_genotype.formula,
						data  = x,
						weights = CombV,
						method = "ML",
						na.action = "na.omit"
					)
				)

				model_null <-
					do.call(
						"gls",
						args = list(
							model = model_null.formula,
							data  = x,
							weights = CombV,
							method = "ML",
							na.action = "na.omit"
						)
					)

				p.value <-
					(anova(model_genotype, model_null)$"p-value"[2])
			} else if (!keep_batch && keep_equalvar) {
				model_genotype <- do.call(
					"gls",
					args = list(
						model = model_genotype.formula,
						data  = x,
						weights = FixV,
						method = "ML",
						na.action = "na.omit"
					)
				)

				model_null <-
					do.call(
						"gls",
						args = list(
							model = model_null.formula,
							data = x,
							weights = FixV,
							method = "ML",
							na.action = "na.omit"
						)
					)

				p.value <-
					(anova(model_genotype, model_null)$"p-value"[2])
			}



			## Final model version with na.exclude and REML method
			if (keep_batch && keep_equalvar) {
				## Model 1
				model_genotype <-
					do.call(
						"lme",
						args = list(
							fixed = model_genotype.formula,
							random =  ~ 1 |
								Batch,
							data = x,
							weights =  FixV,
							na.action = "na.exclude",
							method = "REML"
						)
					)
			} else if (keep_batch && !keep_equalvar) {
				## Model 2
				model_genotype <-
					do.call(
						"lme",
						args = list(
							fixed  = model_genotype.formula,
							random =  ~ 1 |
								Batch,
							data = x,
							weights = CombV,
							na.action = "na.exclude",
							method = "REML"
						)
					)
			} else if (!keep_batch && !keep_equalvar) {
				## Model 2A
				model_genotype <-
					do.call(
						"gls",
						args = list(
							model = model_genotype.formula,
							data  = x,
							weights = CombV,
							na.action = "na.exclude"
						)
					)
			} else if (!keep_batch && keep_equalvar) {
				## Model 1A
				model_genotype <-
					do.call(
						"gls",
						args = list(
							model = model_genotype.formula,
							data  = x,
							weights   = FixV,
							na.action = "na.exclude"
						)
					)
			}


			## Store the results after the final MM modelling step
			linearRegressionOutput$model.output <-
				model_genotype
			linearRegressionOutput$model.null <- model_null
			linearRegressionOutput$model.output.genotype.nulltest.pVal <-
				p.value
			linearRegressionOutput$model.formula.null <-
				model_null.formula
			linearRegressionOutput$model.formula.genotype <-
				model_genotype.formula
			#linearRegressionOutput$model.effect.variance <- keep_equalvar



			## Parse modeloutput and choose output depending on model
			linearRegressionOutput$model.output.summary <-
				parserOutputSummary(linearRegressionOutput)

			# Percentage changes - is the ratio of the genotype effect for a sex relative to
			# the wildtype signal for that variable for that sex - calculation
			mean_list <-
				linearRegressionOutput$model.output.averageRefGenotype
			denominator <- mean_list[1]

			if (linearRegressionOutput$numberSexes == 2) {
				# without weight
				if (is.na(linearRegressionOutput$model.output.summary['weight_estimate'])) {
					if (!is.na(linearRegressionOutput$model.output.summary['sex_estimate']) &&
							!is.na(linearRegressionOutput$model.output.summary['sex_FvKO_estimate']))
					{
						denominator_f <-
							linearRegressionOutput$model.output.summary['intercept_estimate']
						denominator_m <-
							linearRegressionOutput$model.output.summary['intercept_estimate'] +
							linearRegressionOutput$model.output.summary['sex_estimate']
						denominator_f <- denominator
						denominator_m <- denominator
						ratio_f <-
							linearRegressionOutput$model.output.summary['sex_FvKO_estimate'] / denominator_f
						ratio_m <-
							linearRegressionOutput$model.output.summary['sex_MvKO_estimate'] / denominator_m
					}
					else if (!is.na(linearRegressionOutput$model.output.summary['sex_FvKO_estimate']))
					{
						#denominator <- linearRegressionOutput$model.output.summary['intercept_estimate']
						ratio_f <-
							linearRegressionOutput$model.output.summary['sex_FvKO_estimate'] / denominator
						ratio_m <-
							linearRegressionOutput$model.output.summary['sex_MvKO_estimate'] / denominator
					}
					else if (!is.na(linearRegressionOutput$model.output.summary['sex_estimate']))
					{
						denominator_f <-
							linearRegressionOutput$model.output.summary['intercept_estimate']
						denominator_m <-
							linearRegressionOutput$model.output.summary['intercept_estimate'] +
							linearRegressionOutput$model.output.summary['sex_estimate']
						denominator_f <- denominator
						denominator_m <- denominator
						ratio_f <-
							linearRegressionOutput$model.output.summary['genotype_estimate'] / denominator_f
						ratio_m <-
							linearRegressionOutput$model.output.summary['genotype_estimate'] / denominator_m
					}
					else
					{
						#denominator <- linearRegressionOutput$model.output.summary['intercept_estimate']
						ratio_f <-
							linearRegressionOutput$model.output.summary['genotype_estimate'] / denominator
						ratio_m <- ratio_f
					}
				}
				# with weight
				else{
					mean_list <- linearRegressionOutput$model.output.averageRefGenotype
					denominator_f <- mean_list[2]
					denominator_m <- mean_list[3]
					denominator_f <- denominator
					denominator_m <- denominator
					if (!is.na(linearRegressionOutput$model.output.summary['sex_estimate']) &&
							!is.na(linearRegressionOutput$model.output.summary['sex_FvKO_estimate']))
					{
						ratio_f <-
							linearRegressionOutput$model.output.summary['sex_FvKO_estimate'] / denominator_f
						ratio_m <-
							linearRegressionOutput$model.output.summary['sex_MvKO_estimate'] / denominator_m
					}
					else if (!is.na(linearRegressionOutput$model.output.summary['sex_FvKO_estimate']))
					{
						ratio_f <-
							linearRegressionOutput$model.output.summary['sex_FvKO_estimate'] / denominator_f
						ratio_m <-
							linearRegressionOutput$model.output.summary['sex_MvKO_estimate'] / denominator_m
					}
					else
					{
						ratio_f <-
							linearRegressionOutput$model.output.summary['genotype_estimate'] / denominator_f
						ratio_m <-
							linearRegressionOutput$model.output.summary['genotype_estimate'] / denominator_m
					}
				}

				linearRegressionOutput$model.output.percentageChanges <-
					c(ratio_f * 100, ratio_m * 100)
				names(linearRegressionOutput$model.output.percentageChanges) <-
					c('female*genotype ratio', 'male*genotype ratio')
				#[email protected] <- linearRegressionOutput
				#finalResult <- finalOutput
			}
			else{
				# without weight
				if (is.na(linearRegressionOutput$model.output.summary['weight_estimate'])) {
					#denominator <- linearRegressionOutput$model.output.summary['intercept_estimate']
					ratio_f <-
						linearRegressionOutput$model.output.summary['genotype_estimate'] / denominator
				}
				# with weight
				else{
					#mean_list <- linearRegressionOutput$model.output.averageRefGenotype
					#denominator <- mean_list[1]
					ratio_f <-
						linearRegressionOutput$model.output.summary['genotype_estimate'] / denominator
				}

				linearRegressionOutput$model.output.percentageChanges <-
					c(ratio_f * 100)
				names(linearRegressionOutput$model.output.percentageChanges) <-
					c('all*genotype ratio')

			}
			# end of percentage changes calculation


			finalResult@analysisResults <-
				linearRegressionOutput
			## Assign MM quality of fit
			finalResult@analysisResults$model.output.quality <-
				testFinalModel(finalResult)
			finalResult
		},
		# End of tryCatch statement - if fails try to suggest smth useful for the user
		#error = function(e){}
		error = function(error_mes) {
			message(error_mes)
			finalResult <- NULL
			if (linearRegressionOutput$equation == "withWeight")
				stop_message <-
				paste(
					"Error:\nCan't fit the model ",
					format(model_genotype.formula),
					". Try MM with equation 'withoutWeight'. ",
					"Another option is jitter\n",
					sep = ""
				)
			else
				stop_message <-
				paste(
					"Error:\nCan't fit the model ",
					format(model_genotype.formula),
					". Try to add jitter or RR plus method.\n",
					sep = ""
				)

			if (outputMessages) {
				message(stop_message)
				opt <- options(show.error.messages = FALSE)
				on.exit(options(opt))
				stop()
			}
			else {
				stop(stop_message)
			}

		},
		warning = function(w) {
			message(w)
		})
		return(finalResult)
	}
##------------------------------------------------------------------------------
## Parses model output summary and returns in readable vector format
parserOutputSummary <- function(linearRegressionOutput)
{
	result <- linearRegressionOutput
	modeloutput_summary <- summary(result$model.output)

	# Set all values to NA initially prior to selecting those relevant to the model
	genotype_estimate <- NA
	genotype_estimate_SE <- NA
	genotype_p_value <- NA
	sex_estimate <- NA
	sex_estimate_SE <- NA
	sex_p_value <- NA
	intercept_estimate <- NA
	intercept_estimate_SE <- NA
	weight_estimate <- NA
	weight_estimate_SE <- NA
	weight_p_value <- NA
	sex_FvKO_estimate <- NA
	sex_FvKO_SE <- NA
	sex_FvKO_p_value <- NA
	sex_MvKO_estimate <- NA
	sex_MvKO_SE <- NA
	sex_MvKO_p_value <- NA

	# Calculate the length of the table of the fitted model summary
	lengthoftable <- {
		table_length <- NA
		if (result$equation == "withWeight") {
			if (result$numberSexes == 2) {
				if ((result$model.effect.sex
						 && result$model.effect.interaction) |
						(!result$model.effect.sex &&
						 result$model.effect.interaction)) {
					table_length <- 5
				} else if (result$model.effect.sex &&
									 !result$model.effect.interaction) {
					table_length <- 4
				} else {
					table_length <- 3
				}
			} else{
				table_length <- 3
			}
		}
		## Eq.1
		else{
			if (result$numberSexes == 2) {
				if ((result$model.effect.sex
						 && result$model.effect.interaction) |
						(!result$model.effect.sex &&
						 result$model.effect.interaction)) {
					table_length <- 4
				} else if (!result$model.effect.sex &&
									 !result$model.effect.interaction) {
					table_length <- 2
				} else{
					table_length <- 3
				}
			} else{
				table_length <- 2
			}
		}

		table_length
	}

	# Decision tree to pull the information depending on the final model fitted
	# note position is dependent on equation and presence of effects
	switch(result$equation,
				 withoutWeight = {
				 	sex_index <-
				 		match(c("SexMale"), row.names(modeloutput_summary[["tTable"]]))
				 	sex_FvKO_index <- 3
				 	sex_MvKO_index <- 4
				 	if (is.na(sex_index)) {
				 		sex_index <-
				 			match(c("SexFemale"), row.names(modeloutput_summary[["tTable"]]))
				 		sex_FvKO_index <- 4
				 		sex_MvKO_index <- 3
				 	}
				 	if (result$model.effect.batch) {
				 		## for mixed model
				 		intercept_estimate = modeloutput_summary[["tTable"]][[1]]
				 		intercept_estimate_SE = modeloutput_summary[["tTable"]][[(1 +
				 																																lengthoftable)]]
				 		if ((result$model.effect.sex &&
				 				 result$model.effect.interaction)
				 				|
				 				(!result$model.effect.sex &&
				 				 result$model.effect.interaction)) {
				 			sex_estimate = modeloutput_summary[["tTable"]][[sex_index]]
				 			sex_estimate_SE = modeloutput_summary[["tTable"]][[(sex_index +
				 																														lengthoftable)]]
				 			sex_p_value = modeloutput_summary[["tTable"]][[(sex_index +
				 																												4 * lengthoftable)]]
				 			sex_FvKO_estimate = modeloutput_summary[["tTable"]][[sex_FvKO_index]]
				 			sex_FvKO_SE = modeloutput_summary[["tTable"]][[(sex_FvKO_index +
				 																												lengthoftable)]]
				 			sex_FvKO_p_value = modeloutput_summary[["tTable"]][[(sex_FvKO_index +
				 																													 	4 * lengthoftable)]]
				 			sex_MvKO_estimate = modeloutput_summary[["tTable"]][[sex_MvKO_index]]
				 			sex_MvKO_SE = modeloutput_summary[["tTable"]][[(sex_MvKO_index +
				 																												lengthoftable)]]
				 			sex_MvKO_p_value = modeloutput_summary[["tTable"]][[(sex_MvKO_index +
				 																													 	4 * lengthoftable)]]
				 		} else if (!result$model.effect.sex &&
				 							 !result$model.effect.interaction) {
				 			genotype_estimate = modeloutput_summary[["tTable"]][[2]]
				 			genotype_estimate_SE = modeloutput_summary[["tTable"]][[(2 +
				 																															 	lengthoftable)]]
				 			genotype_p_value =  modeloutput_summary[["tTable"]][[(2 +
				 																															4 * lengthoftable)]]
				 		} else{
				 			genotype_estimate = modeloutput_summary[["tTable"]][[2]]
				 			genotype_estimate_SE = modeloutput_summary[["tTable"]][[(2 +
				 																															 	lengthoftable)]]
				 			genotype_p_value =  modeloutput_summary[["tTable"]][[(2 +
				 																															4 * lengthoftable)]]
				 			sex_estimate = modeloutput_summary[["tTable"]][[sex_index]]
				 			sex_estimate_SE = modeloutput_summary[["tTable"]][[(sex_index +
				 																														lengthoftable)]]
				 			sex_p_value = modeloutput_summary[["tTable"]][[(sex_index +
				 																												4 * lengthoftable)]]
				 		}
				 	} else{
				 		## Adaption for being a linear model rather than a mixed model
				 		intercept_estimate = modeloutput_summary[["tTable"]][[1]]
				 		intercept_estimate_SE = modeloutput_summary[["tTable"]][[(1 +
				 																																lengthoftable)]]
				 		if ((result$model.effect.sex &&
				 				 result$model.effect.interaction)
				 				|
				 				(!result$model.effect.sex &&
				 				 result$model.effect.interaction)) {
				 			sex_estimate = modeloutput_summary[["tTable"]][[sex_index]]
				 			sex_estimate_SE = modeloutput_summary[["tTable"]][[(sex_index +
				 																														lengthoftable)]]
				 			sex_p_value = modeloutput_summary[["tTable"]][[(sex_index +
				 																												3 * lengthoftable)]]
				 			sex_FvKO_estimate = modeloutput_summary[["tTable"]][[sex_FvKO_index]]
				 			sex_FvKO_SE = modeloutput_summary[["tTable"]][[(sex_FvKO_index +
				 																												lengthoftable)]]
				 			sex_FvKO_p_value = modeloutput_summary[["tTable"]][[(sex_FvKO_index +
				 																													 	3 * lengthoftable)]]
				 			sex_MvKO_estimate = modeloutput_summary[["tTable"]][[sex_MvKO_index]]
				 			sex_MvKO_SE = modeloutput_summary[["tTable"]][[(sex_MvKO_index +
				 																												lengthoftable)]]
				 			sex_MvKO_p_value = modeloutput_summary[["tTable"]][[(sex_MvKO_index +
				 																													 	3 * lengthoftable)]]
				 		} else if (!result$model.effect.sex &&
				 							 !result$model.effect.interaction) {
				 			genotype_estimate = modeloutput_summary[["tTable"]][[2]]
				 			genotype_estimate_SE = modeloutput_summary[["tTable"]][[(2 +
				 																															 	lengthoftable)]]
				 			genotype_p_value =  modeloutput_summary[["tTable"]][[(2 +
				 																															3 * lengthoftable)]]
				 		} else{
				 			genotype_estimate = modeloutput_summary[["tTable"]][[2]]
				 			genotype_estimate_SE = modeloutput_summary[["tTable"]][[(2 +
				 																															 	lengthoftable)]]
				 			genotype_p_value =  modeloutput_summary[["tTable"]][[(2 +
				 																															3 * lengthoftable)]]
				 			sex_estimate = modeloutput_summary[["tTable"]][[sex_index]]
				 			sex_estimate_SE = modeloutput_summary[["tTable"]][[(sex_index +
				 																														lengthoftable)]]
				 			sex_p_value = modeloutput_summary[["tTable"]][[(3 +
				 																												3 * lengthoftable)]]
				 		}
				 	}
				 },
				 withWeight = {
				 	sex_index <-
				 		match(c("SexMale"), row.names(modeloutput_summary[["tTable"]]))
				 	sex_FvKO_index <- 4
				 	sex_MvKO_index <- 5
				 	if (is.na(sex_index)) {
				 		sex_index <-
				 			match(c("SexFemale"), row.names(modeloutput_summary[["tTable"]]))
				 		sex_FvKO_index <- 5
				 		sex_MvKO_index <- 4
				 	}
				 	if (!result$model.effect.weight) {
				 		## If weight is not significant then the output is the
				 		## same as fitting model Eq1 and so no output is needed.
				 		result$model.effect.batch = NA
				 	} else{
				 		if (result$model.effect.batch) {
				 			## For mixed model
				 			intercept_estimate = modeloutput_summary[["tTable"]][[1]]
				 			intercept_estimate_SE = modeloutput_summary[["tTable"]][[(1 +
				 																																	lengthoftable)]]
				 			if ((
				 				result$model.effect.weight && result$model.effect.sex &&
				 				result$model.effect.interaction
				 			) |
				 			(
				 				result$model.effect.weight &&
				 				!result$model.effect.sex &&
				 				result$model.effect.interaction
				 			)) {
				 				sex_estimate = modeloutput_summary[["tTable"]][[sex_index]]
				 				sex_estimate_SE = modeloutput_summary[["tTable"]][[(sex_index +
				 																															lengthoftable)]]
				 				sex_p_value = modeloutput_summary[["tTable"]][[(sex_index +
				 																													4 * lengthoftable)]]
				 				sex_FvKO_estimate = modeloutput_summary[["tTable"]][[sex_FvKO_index]]
				 				sex_FvKO_SE = modeloutput_summary[["tTable"]][[(sex_FvKO_index +
				 																													lengthoftable)]]
				 				sex_FvKO_p_value = modeloutput_summary[["tTable"]][[(sex_FvKO_index +
				 																														 	4 * lengthoftable)]]
				 				sex_MvKO_estimate = modeloutput_summary[["tTable"]][[sex_MvKO_index]]
				 				sex_MvKO_SE = modeloutput_summary[["tTable"]][[(sex_MvKO_index +
				 																													lengthoftable)]]
				 				sex_MvKO_p_value = modeloutput_summary[["tTable"]][[(sex_MvKO_index +
				 																														 	4 * lengthoftable)]]
				 				weight_estimate = modeloutput_summary[["tTable"]][[3]]
				 				weight_estimate_SE = modeloutput_summary[["tTable"]][[(3 +
				 																															 	lengthoftable)]]
				 				weight_p_value = modeloutput_summary[["tTable"]][[(3 +
				 																													 	4 * lengthoftable)]]
				 			} else if (result$model.effect.weight &&
				 								 !result$model.effect.sex &&
				 								 !result$model.effect.interaction) {
				 				genotype_estimate = modeloutput_summary[["tTable"]][[2]]
				 				genotype_estimate_SE = modeloutput_summary[["tTable"]][[(2 +
				 																																 	lengthoftable)]]
				 				genotype_p_value =  modeloutput_summary[["tTable"]][[(2 +
				 																																4 * lengthoftable)]]
				 				weight_estimate = modeloutput_summary[["tTable"]][[3]]
				 				weight_estimate_SE = modeloutput_summary[["tTable"]][[(3 +
				 																															 	lengthoftable)]]
				 				weight_p_value = modeloutput_summary[["tTable"]][[(3 +
				 																													 	4 * lengthoftable)]]
				 			} else if (result$model.effect.weight &&
				 								 result$model.effect.sex &&
				 								 !result$model.effect.interaction) {
				 				genotype_estimate = modeloutput_summary[["tTable"]][[2]]
				 				genotype_estimate_SE = modeloutput_summary[["tTable"]][[(2 +
				 																																 	lengthoftable)]]
				 				genotype_p_value =  modeloutput_summary[["tTable"]][[(2 +
				 																																4 * lengthoftable)]]
				 				sex_estimate = modeloutput_summary[["tTable"]][[3]]
				 				sex_estimate_SE = modeloutput_summary[["tTable"]][[(3 +
				 																															lengthoftable)]]
				 				sex_p_value = modeloutput_summary[["tTable"]][[(3 +
				 																													4 * lengthoftable)]]
				 				weight_estimate = modeloutput_summary[["tTable"]][[4]]
				 				weight_estimate_SE = modeloutput_summary[["tTable"]][[(4 +
				 																															 	lengthoftable)]]
				 				weight_p_value = modeloutput_summary[["tTable"]][[(4 +
				 																													 	4 * lengthoftable)]]
				 			}
				 		} else{
				 			## Adaption for being a linear model rather than a mixed model
				 			intercept_estimate = modeloutput_summary[["tTable"]][[1]]
				 			intercept_estimate_SE = modeloutput_summary[["tTable"]][[(1 +
				 																																	lengthoftable)]]
				 			if ((
				 				result$model.effect.weight &&
				 				result$model.effect.sex &&
				 				result$model.effect.interaction
				 			) |
				 			(
				 				result$model.effect.weight &&
				 				!result$model.effect.sex &&
				 				result$model.effect.interaction
				 			)) {
				 				sex_estimate = modeloutput_summary[["tTable"]][[sex_index]]
				 				sex_estimate_SE = modeloutput_summary[["tTable"]][[(sex_index +
				 																															lengthoftable)]]
				 				sex_p_value = modeloutput_summary[["tTable"]][[(sex_index +
				 																													3 * lengthoftable)]]
				 				weight_estimate = modeloutput_summary[["tTable"]][[3]]
				 				weight_estimate_SE = modeloutput_summary[["tTable"]][[(3 +
				 																															 	lengthoftable)]]
				 				weight_p_value = modeloutput_summary[["tTable"]][[(3 +
				 																													 	3 * lengthoftable)]]
				 				sex_FvKO_estimate = modeloutput_summary[["tTable"]][[sex_FvKO_index]]
				 				sex_FvKO_SE = modeloutput_summary[["tTable"]][[(sex_FvKO_index +
				 																													lengthoftable)]]
				 				sex_FvKO_p_value = modeloutput_summary[["tTable"]][[(sex_FvKO_index +
				 																														 	3 * lengthoftable)]]
				 				sex_MvKO_estimate = modeloutput_summary[["tTable"]][[sex_MvKO_index]]
				 				sex_MvKO_SE = modeloutput_summary[["tTable"]][[(sex_MvKO_index +
				 																													lengthoftable)]]
				 				sex_MvKO_p_value = modeloutput_summary[["tTable"]][[(sex_MvKO_index +
				 																														 	3 * lengthoftable)]]

				 			} else if (result$model.effect.weight &&
				 								 result$model.effect.sex &&
				 								 !result$model.effect.interaction) {
				 				genotype_estimate = modeloutput_summary[["tTable"]][[2]]
				 				genotype_estimate_SE = modeloutput_summary[["tTable"]][[(2 +
				 																																 	lengthoftable)]]
				 				genotype_p_value =  modeloutput_summary[["tTable"]][[(2 +
				 																																3 * lengthoftable)]]
				 				sex_estimate = modeloutput_summary[["tTable"]][[sex_index]]
				 				sex_estimate_SE = modeloutput_summary[["tTable"]][[(sex_index +
				 																															lengthoftable)]]
				 				sex_p_value = modeloutput_summary[["tTable"]][[(sex_index +
				 																													3 * lengthoftable)]]
				 				weight_estimate = modeloutput_summary[["tTable"]][[4]]
				 				weight_estimate_SE = modeloutput_summary[["tTable"]][[(4 +
				 																															 	lengthoftable)]]
				 				weight_p_value = modeloutput_summary[["tTable"]][[(4 +
				 																													 	3 * lengthoftable)]]
				 			} else if (result$model.effect.weight &&
				 								 !result$model.effect.sex &&
				 								 !result$model.effect.interaction) {
				 				genotype_estimate = modeloutput_summary[["tTable"]][[2]]
				 				genotype_estimate_SE = modeloutput_summary[["tTable"]][[(2 +
				 																																 	lengthoftable)]]
				 				genotype_p_value =  modeloutput_summary[["tTable"]][[(2 +
				 																																3 * lengthoftable)]]
				 				weight_estimate = modeloutput_summary[["tTable"]][[3]]
				 				weight_estimate_SE = modeloutput_summary[["tTable"]][[(3 +
				 																															 	lengthoftable)]]
				 				weight_p_value = modeloutput_summary[["tTable"]][[(3 +
				 																													 	3 * lengthoftable)]]
				 			}
				 		}
				 	}
				 })
	output <-
		c(
			genotype_estimate,
			genotype_estimate_SE,
			genotype_p_value,
			sex_estimate,
			sex_estimate_SE,
			sex_p_value,
			weight_estimate,
			weight_estimate_SE,
			weight_p_value,
			intercept_estimate,
			intercept_estimate_SE,
			sex_FvKO_estimate,
			sex_FvKO_SE,
			sex_FvKO_p_value,
			sex_MvKO_estimate,
			sex_MvKO_SE,
			sex_MvKO_p_value
		)
	names(output) <- c(
		"genotype_estimate",
		"genotype_estimate_SE",
		"genotype_p_value",
		"sex_estimate",
		"sex_estimate_SE",
		"sex_p_value",
		"weight_estimate",
		"weight_estimate_SE",
		"weight_p_value",
		"intercept_estimate",
		"intercept_estimate_SE",
		"sex_FvKO_estimate",
		"sex_FvKO_SE",
		"sex_FvKO_p_value",
		"sex_MvKO_estimate",
		"sex_MvKO_SE",
		"sex_MvKO_p_value"
	)

	return(output)
}
##------------------------------------------------------------------------------

Try the PhenStat package in your browser

Any scripts or data that you put into this service are public.

PhenStat documentation built on Sept. 8, 2018, 6:01 p.m.