Nothing
R/MMFramework.R
In PhenStat: Statistical analysis of phenotypic data
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)
}
equalvar_pvalue = batch_pvalue = NA #HAMED 2/1/2019
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
batch_pvalue <- p.value.batch #HAMED 2/1/2019
## 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
equalvar_pvalue <- p.value.variance #HAMED 2/1/2019
}
} 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: ",
"keep_Batch=",
keep_batch,
", keep_EqualVariance=",
keep_equalvar,
", keep_Weight=",
keep_weight,
", keep_Sex=",
keep_sex,
", keep_Interaction=",
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: ",
"keep_Batch=",
user_keep_batch,
", keep_EqualVariance=",
user_keep_equalvar,
", keep_Weight=",
user_keep_weight,
", keep_Sex=",
user_keep_sex,
", keep_Interaction=",
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.batch.pvalue = batch_pvalue, #HAMED 2/1/2019
model.effect.variance = keep_equalvar,
model.variance.pvalue = equalvar_pvalue, #HAMED 2/1/2019
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')
#finalOutput@result <- 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 Nov. 8, 2020, 8:13 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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)
}
equalvar_pvalue = batch_pvalue = NA #HAMED 2/1/2019
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
batch_pvalue <- p.value.batch #HAMED 2/1/2019
## 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
equalvar_pvalue <- p.value.variance #HAMED 2/1/2019
}
} 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: ",
"keep_Batch=",
keep_batch,
", keep_EqualVariance=",
keep_equalvar,
", keep_Weight=",
keep_weight,
", keep_Sex=",
keep_sex,
", keep_Interaction=",
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: ",
"keep_Batch=",
user_keep_batch,
", keep_EqualVariance=",
user_keep_equalvar,
", keep_Weight=",
user_keep_weight,
", keep_Sex=",
user_keep_sex,
", keep_Interaction=",
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.batch.pvalue = batch_pvalue, #HAMED 2/1/2019
model.effect.variance = keep_equalvar,
model.variance.pvalue = equalvar_pvalue, #HAMED 2/1/2019
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')
#finalOutput@result <- 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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.