Nothing
## 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.
##------------------------------------------------------------------------------
## Diagnostictest.R contains testFinalModel function
##------------------------------------------------------------------------------
## Diagnostic test output for MM quality of fit. There are no arguments checks
## assuming that function is called internally from the finalModel function.
## Otherwise should be used with precaution.
testFinalModel <- function(phenTestResult)
{
result <- phenTestResult
x <- analysedDataset(result)
depVariable <- getVariable(result)
analysisResults <- analysisResults(phenTestResult)
equation <- analysisResults$equation
keep_weight <- analysisResults$model.effect.weight
keep_sex <- analysisResults$model.effect.sex
keep_interaction <- analysisResults$model.effect.interaction
keep_batch <- analysisResults$model.effect.batch
keep_equalvar <- analysisResults$model.effect.variance
a <- levels(x$Genotype)
numberofsexes <- analysisResults$numberSexes
if (!keep_weight && equation == "withWeight") {
testresults <- c(a[1], NA, a[2], NA, NA, NA)
} else if ('Batch' %in% colnames(x)) {
res <- resid(analysisResults$model.output)
data_all <- data.frame(x, res)
genotype_no <- length(a)
data_all[, c("Sex", "Batch")] <- lapply(data_all[, c("Sex",
"Batch")], factor)
No_batches <- nlevels(data_all$Batch)
outputnumeric <-
is.numeric(analysisResults$model.output$apVar)
Gp1 <- subset(data_all, data_all$Genotype == a[1])
Gp2 <- subset(data_all, data_all$Genotype == a[2])
No_Gp1 <- sum(is.finite(Gp1[, c("res")]))
No_Gp2 <- sum(is.finite(Gp2[, c("res")]))
if (keep_batch &&
No_batches > 7 &&
outputnumeric &&
!is.null(analysisResults$model.output) &&
is(analysisResults$model.output,'lme')) {
blups <- ranef(analysisResults$model.output)
blups_test <-
suppressWarnings(cvm.test(blups [, 1])$p.value)
sdests <- exp(attr(analysisResults$model.output$apVar,
"Pars")) #extract variance estimates
Zbat <- model.matrix( ~ Batch,
model.frame(~ Batch,
analysisResults$model.output$groups)) #create random effects design matrix
ycov <-
(Zbat %*% t(Zbat)) * sdests["reStruct.Batch"] ^ 2 + diag(rep
(1, nrow(
analysisResults$model.output$groups
))) * sdests["lSigma"] ^ 2 #create estimated cov(y)
Lt <-
chol(solve(ycov)) #Cholesky decomposition of inverse of cov
# (y) (see Houseman '04 eq. (2))
rotres <-
Lt %*% analysisResults$model.output$residuals[, "fixed"] #rotated residuals
rotated_residual_test <-
suppressWarnings(cvm.test(rotres)$p.value)
} else{
blups_test <- NA
rotated_residual_test <- NA
}
if (No_Gp1 > 7) {
gp1_norm_res <- suppressWarnings(cvm.test(Gp1$res)$p.value)
} else{
gp1_norm_res <- NA
}
if (No_Gp2 > 7) {
gp2_norm_res <- suppressWarnings(cvm.test(Gp2$res)$p.value)
} else{
gp2_norm_res <- NA
}
testresults <-
c(a[1],
gp1_norm_res,
a[2],
gp2_norm_res,
blups_test,
rotated_residual_test)
}
else{
testresults <- c(a[1], NA, a[2], NA, NA, NA)
}
return(testresults)
}
##------------------------------------------------------------------------------
testFinalLRModel <- function(phenTestResult)
{
result <- phenTestResult
x <- analysedDataset(result)
depVariable <- getVariable(result)
analysisResults <- analysisResults(phenTestResult)
keep_sex <- analysisResults$model.effect.sex
keep_interaction <- analysisResults$model.effect.interaction
keep_batch <- analysisResults$model.effect.batch
a <- levels(x$Genotype)
numberofsexes <- analysisResults$numberSexes
blups_test <- NA
rotated_residual_test <- NA
gp1_norm_res <- NA
gp2_norm_res <- NA
testResults <-
c(a[1],
gp1_norm_res,
a[2],
gp2_norm_res,
blups_test,
rotated_residual_test)
return(testResults)
}
##------------------------------------------------------------------------------
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.