# Creates two lists of lists. First has training samples, second has test samples for a range
# of different cross-validation schemes.
#' @import utils
.samplesSplits <- function(crossValParams, outcome)
{
if(crossValParams@samplesSplits %in% c("k-Fold", "Permute k-Fold"))
{
nPermutations <- ifelse(crossValParams@samplesSplits == "k-Fold", 1, crossValParams@permutations)
nFolds <- crossValParams@folds
samplesFolds <- lapply(1:nPermutations, function(permutation)
{
# Create maximally-balanced folds, so class balance is about the same in all.
allFolds <- vector(mode = "list", length = nFolds)
foldsIndexes <- rep(1:nFolds, length.out = length(outcome))
foldsIndex = 1
# Dummy encoding for when outcome is not a class.
if(is(outcome, "Surv")) outcome <- factor(rep("a", length(outcome)))
for(outcomeName in levels(outcome))
{
# Permute the indexes of samples in the class.
whichSamples <- sample(which(outcome == outcomeName))
whichFolds <- foldsIndexes[foldsIndex:(foldsIndex + length(whichSamples) - 1)]
# Put each sample into its fold.
for(sampleIndex in 1:length(whichSamples))
{
allFolds[[whichFolds[sampleIndex]]] <- c(allFolds[[whichFolds[sampleIndex]]], whichSamples[sampleIndex])
}
# Move the beginning index to the first new index.
foldsIndex <- foldsIndex + length(whichSamples)
}
list(train = lapply(1:nFolds, function(index) unlist(allFolds[setdiff(1:nFolds, index)])),
test = allFolds
)
})
# Reorganise into two separate lists, no more nesting.
list(train = unlist(lapply(samplesFolds, '[[', 1), recursive = FALSE),
test = unlist(lapply(samplesFolds, '[[', 2), recursive = FALSE))
} else if(crossValParams@samplesSplits == "Permute Percentage Split") {
# Take the same percentage of samples from each class to be in training set.
percent <- crossValParams@percentTest
samplesTrain <- round((100 - percent) / 100 * table(outcome))
samplesTest <- round(percent / 100 * table(outcome))
samplesLists <- lapply(1:crossValParams@permutations, function(permutation)
{
trainSet <- unlist(mapply(function(outcomeName, number)
{
sample(which(outcome == outcomeName), number)
}, levels(outcome), samplesTrain))
testSet <- setdiff(1:length(classes), trainSet)
list(trainSet, testSet)
})
# Reorganise into two lists: training, testing.
list(train = lapply(samplesLists, "[[", 1), test = lapply(samplesLists, "[[", 2))
} else if(crossValParams@samplesSplits == "Leave-k-Out") { # leave k out.
testSamples <- as.data.frame(utils::combn(length(outcome), crossValParams@leave))
trainingSamples <- lapply(testSamples, function(sample) setdiff(1:length(outcome), sample))
list(train = as.list(trainingSamples), test = as.list(testSamples))
}
}
# Creates a two-column table for tracking the permutation, fold number, or subset of each set
# of test samples (i.e. for leave-k-out scheme).
.splitsTestInfo <- function(crossValParams, samplesSplits)
{
permutationIDs <- NULL
foldIDs <- NULL
subsetIDs <- NULL
if(crossValParams@samplesSplits %in% c("k-Fold", "Permute k-Fold"))
{
foldsSamples <- lengths(samplesSplits[[2]][1:crossValParams@folds])
totalSamples <- sum(foldsSamples)
if(crossValParams@samplesSplits == "Permute k-Fold")
permutationIDs <- rep(1:crossValParams@permutations, each = totalSamples)
times <- ifelse(is.null(crossValParams@permutations), 1, crossValParams@permutations)
foldIDs <- rep(rep(1:crossValParams@folds, foldsSamples), times = times)
} else if(crossValParams@samplesSplits == "Permute Percentage Split") {
permutationIDs <- rep(1:crossValParams@permutations, each = length(samplesSplits[[2]][[1]]))
} else { # Leave-k-out
totalSamples <- length(unique(unlist(samplesSplits[[2]])))
subsetIDs <- rep(1:choose(totalSamples, crossValParams@leave), each = crossValParams@leave)
}
summaryTable <- cbind(permutation = permutationIDs, fold = foldIDs, subset = subsetIDs)
}
# Add extra variables from within runTest functions to function specified by a params object.
.addIntermediates <- function(params)
{
intermediateName <- params@intermediate
intermediates <- list(dynGet(intermediateName, inherits = TRUE))
if(is.null(names(params@intermediate))) names(intermediates) <- intermediateName else names(intermediates) <- names(params@intermediate)
params@otherParams <- c(params@otherParams, intermediates)
params
}
# Carries out one iteration of feature selection. Basically, a ranking function is used to rank
# the features in the training set from best to worst and different top sets are used either for
# predicting on the training set (resubstitution) or nested cross-validation of the training set,
# to find the set of top features which give the best (user-specified) performance measure.
.doSelection <- function(measurementsTrain, outcomeTrain, crossValParams, modellingParams, verbose)
{
tuneParams <- modellingParams@selectParams@tuneParams
performanceType <- tuneParams[["performanceType"]]
topNfeatures <- tuneParams[["nFeatures"]]
tuneParams <- tuneParams[-match(c("performanceType", "nFeatures"), names(tuneParams))] # Only used as evaluation metric.
# Make selectParams NULL, since we are currently doing selection and it shouldn't call
# itself infinitely, but save the parameters of the ranking function for calling the ranking
# function directly using do.call below.
featureRanking <- modellingParams@selectParams@featureRanking
otherParams <- modellingParams@selectParams@otherParams
doSubset <- modellingParams@selectParams@subsetToSelections
modellingParams@selectParams <- NULL
betterValues <- .ClassifyRenvir[["performanceInfoTable"]][.ClassifyRenvir[["performanceInfoTable"]][, "type"] == performanceType, "better"]
if(is.function(featureRanking)) # Not a list for ensemble selection.
{
paramList <- list(measurementsTrain, outcomeTrain, verbose = verbose)
paramList <- append(paramList, otherParams) # Used directly by a feature ranking function for rankings of features.
if(length(tuneParams) == 0) tuneParams <- list(None = "none")
tuneCombosSelect <- expand.grid(tuneParams, stringsAsFactors = FALSE)
# Generate feature rankings for each one of the tuning parameter combinations.
rankings <- lapply(1:nrow(tuneCombosSelect), function(rowIndex)
{
tuneCombo <- tuneCombosSelect[rowIndex, , drop = FALSE]
if(tuneCombo != "none") # Add real parameters before function call.
paramList <- append(paramList, tuneCombo)
do.call(featureRanking, paramList)
})
if(attr(featureRanking, "name") == "previousSelection") # Actually selection not ranking.
return(list(NULL, rankings[[1]], NULL))
if(crossValParams@tuneMode == "none") # No parameters to choose between.
return(list(NULL, rankings[[1]], NULL))
tuneParamsTrain <- list(topN = topNfeatures)
tuneParamsTrain <- append(tuneParamsTrain, modellingParams@trainParams@tuneParams)
tuneCombosTrain <- expand.grid(tuneParamsTrain, stringsAsFactors = FALSE)
modellingParams@trainParams@tuneParams <- NULL
allPerformanceTables <- lapply(rankings, function(rankingsVariety)
{
# Creates a matrix. Columns are top n features, rows are varieties (one row if None).
performances <- sapply(1:nrow(tuneCombosTrain), function(rowIndex)
{
whichTry <- 1:tuneCombosTrain[rowIndex, "topN"]
if(doSubset)
{
topFeatures <- rankingsVariety[whichTry]
measurementsTrain <- measurementsTrain[, topFeatures, drop = FALSE] # Features in columns
} else { # Pass along features to use.
modellingParams@trainParams@otherParams <- c(modellingParams@trainParams@otherParams, setNames(list(rankingsVariety[whichTry]), names(modellingParams@trainParams@intermediate)))
}
if(ncol(tuneCombosTrain) > 1) # There are some parameters for training.
modellingParams@trainParams@otherParams <- c(modellingParams@trainParams@otherParams, tuneCombosTrain[rowIndex, 2:ncol(tuneCombosTrain), drop = FALSE])
modellingParams@trainParams@intermediate <- character(0)
# Do either resubstitution classification or nested-CV classification and calculate the resulting performance metric.
if(crossValParams@tuneMode == "Resubstitution")
{
# Specify measurementsTrain and outcomeTrain for testing, too.
result <- runTest(measurementsTrain, outcomeTrain, measurementsTrain, outcomeTrain,
crossValParams = NULL, modellingParams = modellingParams,
verbose = verbose, .iteration = "internal")
predictions <- result[["predictions"]]
# Classifiers will use a column "class" and survival models will use a column "risk".
if(class(predictions) == "data.frame")
predictedOutcome <- predictions[, na.omit(match(c("class", "risk"), colnames(predictions)))]
else
predictedOutcome <- predictions
calcExternalPerformance(outcomeTrain, predictedOutcome, performanceType)
} else {
result <- runTests(measurementsTrain, outcomeTrain, crossValParams, modellingParams, verbose = verbose)
result <- calcCVperformance(result, performanceType)
median(performance(result)[[performanceType]])
}
})
bestOne <- ifelse(betterValues == "lower", which.min(performances)[1], which.max(performances)[1])
list(data.frame(tuneCombosTrain, performance = performances), bestOne)
})
tablesBestMetrics <- sapply(allPerformanceTables, function(tableIndexPair) tableIndexPair[[1]][tableIndexPair[[2]], "performance"])
tunePick <- ifelse(betterValues == "lower", which.min(tablesBestMetrics)[1], which.max(tablesBestMetrics)[1])
if(verbose == 3)
message("Features selected.")
tuneDetails <- allPerformanceTables[[tunePick]] # List of length 2.
rankingUse <- rankings[[tunePick]]
selectionIndices <- rankingUse[1:(tuneDetails[[1]][tuneDetails[[2]], "topN"])]
names(tuneDetails) <- c("tuneCombinations", "bestIndex")
colnames(tuneDetails[[1]])[ncol(tuneDetails[[1]])] <- performanceType
list(ranked = rankingUse, selected = selectionIndices, tune = tuneDetails)
} else if(is.list(featureRanking)) { # It is a list of functions for ensemble selection.
featuresIndiciesLists <- mapply(function(selector, selParams)
{
paramList <- list(measurementsTrain, outcomeTrain, trainParams = trainParams,
predictParams = predictParams, verbose = verbose)
paramList <- append(paramList, selParams)
do.call(selector, paramList)
}, modellingParams@selectParams@featureRanking, modellingParams@selectParams@otherParams, SIMPLIFY = FALSE)
performances <- sapply(topNfeatures, function(topN)
{
topIndices <- unlist(lapply(featuresIndiciesLists, function(featuresIndicies) featuresIndicies[1:topN]))
topIndicesCounts <- table(topIndices)
keep <- names(topIndicesCounts)[topIndicesCounts >= modellingParams@selectParams@minPresence]
measurementsTrain <- measurementsTrain[, as.numeric(keep), drop = FALSE] # Features in columns
if(crossValParams@tuneMode == "Resubstitution")
{
result <- runTest(measurementsTrain, outcomeTrain,
measurementsTrain, outcomeTrain,
crossValParams = NULL, modellingParams,
verbose = verbose, .iteration = "internal")
predictions <- result[["predictions"]]
if(class(predictions) == "data.frame")
predictedOutcome <- predictions[, "class"]
else
predictedOutcome <- predictions
calcExternalPerformance(outcomeTrain, predictedOutcome, performanceType)
} else {
result <- runTests(measurementsSubset, outcomeTrain, crossValParams, modellingParams, verbose = verbose)
result <- calcCVperformance(result, performanceType)
median(performance(aResult)[[performanceType]])
}
})
bestOne <- ifelse(betterValues == "lower", which.min(performances)[1], which.max(performances)[1])
selectionIndices <- unlist(lapply(featuresLists, function(featuresList) featuresList[1:topNfeatures[bestOne]]))
names(table(selectionIndices))[table(selectionIndices) >= modellingParams@selectParams@minPresence]
list(NULL, selectionIndices, NULL)
} else { # Previous selection
selectedFeatures <- list(NULL, selectionIndices, NULL)
}
}
# Only for transformations that need to be done within cross-validation.
.doTransform <- function(measurementsTrain, measurementsTest, transformParams, verbose)
{
paramList <- list(measurementsTrain, measurementsTest)
if(length(transformParams@otherParams) > 0)
paramList <- c(paramList, transformParams@otherParams)
paramList <- c(paramList, verbose = verbose)
do.call(transformParams@transform, paramList)
}
# Code to create a function call to a training function. Might also do training and testing
# within the same function, so test samples are also passed in case they are needed.
.doTrain <- function(measurementsTrain, outcomeTrain, measurementsTest, outcomeTest, crossValParams, modellingParams, verbose)
{
tuneDetails <- NULL
if(!is.null(modellingParams@trainParams@tuneParams) && is.null(modellingParams@selectParams))
{
performanceType <- modellingParams@trainParams@tuneParams[["performanceType"]]
modellingParams@trainParams@tuneParams <- modellingParams@trainParams@tuneParams[-match("performanceType", names(modellingParams@trainParams@tuneParams))]
tuneCombos <- expand.grid(modellingParams@trainParams@tuneParams, stringsAsFactors = FALSE)
modellingParams@trainParams@tuneParams <- NULL
performances <- sapply(1:nrow(tuneCombos), function(rowIndex)
{
modellingParams@trainParams@otherParams <- c(modellingParams@trainParams@otherParams, as.list(tuneCombos[rowIndex, ]))
if(crossValParams@tuneMode == "Resubstitution")
{
result <- runTest(measurementsTrain, outcomeTrain, measurementsTrain, outcomeTrain,
crossValParams = NULL, modellingParams,
verbose = verbose, .iteration = "internal")
predictions <- result[["predictions"]]
if(class(predictions) == "data.frame")
predictedOutcome <- predictions[, colnames(predictions) %in% c("class", "risk")]
else
predictedOutcome <- predictions
calcExternalPerformance(outcomeTrain, predictedOutcome, performanceType)
} else {
result <- runTests(measurementsTrain, outcomeTrain,
crossValParams, modellingParams,
verbose = verbose, .iteration = "internal")
result <- calcCVperformance(result, performanceType)
median(performances(result)[[performanceType]])
}
})
allPerformanceTable <- data.frame(tuneCombos, performances)
colnames(allPerformanceTable)[ncol(allPerformanceTable)] <- performanceType
betterValues <- .ClassifyRenvir[["performanceInfoTable"]][.ClassifyRenvir[["performanceInfoTable"]][, "type"] == performanceType, "better"]
bestOne <- ifelse(betterValues == "lower", which.min(performances)[1], which.max(performances)[1])
tuneChosen <- tuneCombos[bestOne, , drop = FALSE]
tuneDetails <- list(tuneCombos, bestOne)
names(tuneDetails) <- c("tuneCombinations", "bestIndex")
modellingParams@trainParams@otherParams <- tuneChosen
}
if (!"previousTrained" %in% attr(modellingParams@trainParams@classifier, "name"))
# Don't name these first two variables. Some classifier functions might use classesTrain and others use outcomeTrain.
paramList <- list(measurementsTrain, outcomeTrain)
else # Don't pass the measurements and classes, because a pre-existing classifier is used.
paramList <- list()
if(is.null(modellingParams@predictParams)) # One function does both training and testing.
paramList <- c(paramList, measurementsTest)
if(length(modellingParams@trainParams@otherParams) > 0)
paramList <- c(paramList, modellingParams@trainParams@otherParams)
paramList <- c(paramList, verbose = verbose)
trained <- do.call(modellingParams@trainParams@classifier, paramList)
if(verbose >= 2)
message("Training completed.")
list(model = trained, tune = tuneDetails)
}
# Creates a function call to a prediction function.
.doTest <- function(trained, measurementsTest, predictParams, verbose)
{
if(!is.null(predictParams@predictor))
{
paramList <- list(trained, measurementsTest)
if(length(predictParams@otherParams) > 0) paramList <- c(paramList, predictParams@otherParams)
paramList <- c(paramList, verbose = verbose)
prediction <- do.call(predictParams@predictor, paramList)
} else { prediction <- trained } # Trained is actually the predictions because only one function, not two.
if(verbose >= 2)
message("Prediction completed.")
prediction
}
# Converts the characteristics of cross-validation into a pretty character string.
.validationText <- function(crossValParams)
{
switch(crossValParams@samplesSplits,
`Permute k-Fold` = paste(crossValParams@permutations, "Permutations,", crossValParams@folds, "Folds"),
`k-Fold` = paste(crossValParams@folds, "-fold cross-validation", sep = ''),
`Leave-k-Out` = paste("Leave", crossValParams@leave, "Out"),
`Permute Percentage Split` = paste(crossValParams@permutations, "Permutations,", crossValParams@percentTest, "% Test"),
independent = "Independent Set")
}
# Used for ROC area under curve calculation.
# PRtable is a data frame with columns FPR, TPR and class.
# distinctClasses is a vector of all of the class names.
.calcArea <- function(PRtable, distinctClasses)
{
do.call(rbind, lapply(distinctClasses, function(aClass)
{
classTable <- subset(PRtable, class == aClass)
areaSum <- 0
for(index in 2:nrow(classTable))
{
# Some samples had identical predictions but belong to different classes.
if(classTable[index, "FPR"] != classTable[index - 1, "FPR"] && classTable[index, "TPR"] != classTable[index - 1, "TPR"])
{
newArea <- (classTable[index, "FPR"] - classTable[index - 1, "FPR"]) * classTable[index - 1, "TPR"] + # Rectangle part
0.5 * (classTable[index, "FPR"] - classTable[index - 1, "FPR"]) * (classTable[index, "TPR"] - classTable[index - 1, "TPR"]) # Triangle part on top.
} else { # Only one sample with predicted score. Line went either up or right, but not both.
newArea <- (classTable[index, "FPR"] - classTable[index - 1, "FPR"]) * classTable[index, "TPR"]
}
areaSum <- areaSum + newArea
}
data.frame(classTable, AUC = round(areaSum, 2), check.names = FALSE)
}))
}
# Converts features into strings to be displayed in plots.
.getFeaturesStrings <- function(importantFeatures)
{
# Do nothing if only a simple vector of feature IDs.
if(!is.null(ncol(importantFeatures[[1]]))) # Data set and feature ID columns.
importantFeatures <- lapply(importantFeatures, function(features) paste(features[, 1], features[, 2]))
else if("Pairs" %in% class(importantFeatures[[1]]))
importantFeatures <- lapply(importantFeatures, function(features) paste(first(features), second(features), sep = '-'))
importantFeatures
}
# Function to overwrite characteristics which are automatically derived from function names
# by user-specified values.
.filterCharacteristics <- function(characteristics, autoCharacteristics)
{
# Overwrite automatically-chosen names with user's names.
if(nrow(autoCharacteristics) > 0 && nrow(characteristics) > 0)
{
overwrite <- na.omit(match(characteristics[, "characteristic"], autoCharacteristics[, "characteristic"]))
if(length(overwrite) > 0)
autoCharacteristics <- autoCharacteristics[-overwrite, ]
}
# Merge characteristics tables and return.
rbind(characteristics, autoCharacteristics)
}
# Don't just plot groups alphabetically, but do so in a meaningful order.
.addUserLevels <- function(plotData, orderingList)
{
for(orderingIndex in 1:length(orderingList))
{
plotData[, names(orderingList)[orderingIndex]] <- factor(plotData[, names(orderingList)[orderingIndex]],
levels = orderingList[[orderingIndex]])
}
plotData
}
# Function to identify the parameters of an S4 method.
.methodFormals <- function(f, signature) {
tryCatch({
fdef <- getGeneric(f)
method <- selectMethod(fdef, signature)
genFormals <- base::formals(fdef)
b <- body(method)
if(is(b, "{") && is(b[[2]], "<-") && identical(b[[2]][[2]], as.name(".local"))) {
local <- eval(b[[2]][[3]])
if(is.function(local))
return(formals(local))
warning("Expected a .local assignment to be a function. Corrupted method?")
}
genFormals
},
error = function(error) {
formals(f)
})
}
# Find the x-axis positions where a set of density functions cross-over.
# The trivial cross-overs at the beginning and end of the data range are removed.
# Used by the mixtures of normals and naive Bayes classifiers.
.densitiesCrossover <- function(densities) # A list of densities created by splinefun.
{
if(!all(table(unlist(lapply(densities, function(density) density[['x']]))) == length(densities)))
stop("x positions are not the same for all of the densities.")
lapply(1:length(densities), function(densityIndex) # All crossing points with other class densities.
{
unlist(lapply(setdiff(1:length(densities), densityIndex), function(otherIndex)
{
allDifferences <- densities[[densityIndex]][['y']] - densities[[otherIndex]][['y']]
crosses <- which(diff(sign(allDifferences)) != 0)
crosses <- sapply(crosses, function(cross) # Refine location for plateaus.
{
isSmall <- rle(allDifferences[(cross+1):length(allDifferences)] < 0.000001)
if(isSmall[["values"]][1] == "TRUE")
cross <- cross + isSmall[["lengths"]][1] / 2
cross
})
if(length(crosses) > 1 && densities[[densityIndex]][['y']][crosses[1]] < 0.000001 && densities[[densityIndex]][['y']][crosses[length(crosses)]] < 0.000001)
crosses <- crosses[-c(1, length(crosses))] # Remove crossings at ends of densities.
densities[[densityIndex]][['x']][crosses]
}))
})
}
# Samples in the training set are upsampled or downsampled so that the class imbalance is
# removed.
.rebalanceTrainingClasses <- function(measurementsTrain, classesTrain, balancing)
{
samplesPerClassTrain <- table(classesTrain)
downsampleTo <- min(samplesPerClassTrain)
upsampleTo <- max(samplesPerClassTrain)
trainBalanced <- unlist(mapply(function(classSize, className)
{
if(balancing == "downsample" && classSize > downsampleTo)
sample(which(classesTrain == className), downsampleTo)
else if(balancing == "upsample" && classSize < upsampleTo)
sample(which(classesTrain == className), upsampleTo, replace = TRUE)
else
which(classesTrain == className)
}, samplesPerClassTrain, names(samplesPerClassTrain), SIMPLIFY = FALSE))
measurementsTrain <- measurementsTrain[trainBalanced, ]
classesTrain <- classesTrain[trainBalanced]
list(measurementsTrain = measurementsTrain, classesTrain = classesTrain)
}
.transformKeywordToFunction <- function(keyword)
{
switch(
keyword,
"none" = NULL,
"diffLoc" = subtractFromLocation
)
}
.selectionKeywordToFunction <- function(keyword)
{
switch(
keyword,
"none" = NULL,
"t-test" = differentMeansRanking,
"limma" = limmaRanking,
"edgeR" = edgeRranking,
"Bartlett" = bartlettRanking,
"Levene" = leveneRanking,
"DMD" = DMDranking,
"likelihoodRatio" = likelihoodRatioRanking,
"KS" = KolmogorovSmirnovRanking,
"KL" = KullbackLeiblerRanking,
"CoxPH" = coxphRanking,
"selectMulti" = selectMulti
)
}
.classifierKeywordToParams <- function(keyword)
{
switch(
keyword,
"randomForest" = RFparams(),
"randomSurvivalForest" = RSFparams(),
"GLM" = GLMparams(),
"elasticNetGLM" = elasticNetGLMparams(),
"SVM" = SVMparams(),
"NSC" = NSCparams(),
"DLDA" = DLDAparams(),
"naiveBayes" = naiveBayesParams(),
"mixturesNormals" = mixModelsParams(),
"kNN" = kNNparams(),
"CoxPH" = coxphParams(),
"CoxNet" = coxnetParams()
)
}
.dlda <- function(x, y, prior = NULL){ # Remove this once sparsediscrim is reinstated to CRAN.
obj <- list()
obj$labels <- y
obj$N <- length(y)
obj$p <- ncol(x)
obj$groups <- levels(y)
obj$num_groups <- nlevels(y)
est_mean <- "mle"
# Error Checking
if (!is.null(prior)) {
if (length(prior) != obj$num_groups) {
stop("The number of 'prior' probabilities must match the number of classes in 'y'.")
}
if (any(prior <= 0)) {
stop("The 'prior' probabilities must be nonnegative.")
}
if (sum(prior) != 1) {
stop("The 'prior' probabilities must sum to one.")
}
}
if (any(table(y) < 2)) {
stop("There must be at least 2 observations in each class.")
}
# By default, we estimate the 'a priori' probabilities of class membership with
# the MLEs (the sample proportions).
if (is.null(prior)) {
prior <- as.vector(table(y) / length(y))
}
# For each class, we calculate the MLEs (or specified alternative estimators)
# for each parameter used in the DLDA classifier. The 'est' list contains the
# estimators for each class.
obj$est <- tapply(seq_along(y), y, function(i) {
stats <- list()
stats$n <- length(i)
stats$xbar <- colMeans(x[i, , drop = FALSE])
stats$var <- with(stats, (n - 1) / n * apply(x[i, , drop = FALSE], 2, var))
stats
})
# Calculates the pooled variance across all classes.
obj$var_pool <- Reduce('+', lapply(obj$est, function(x) x$n * x$var)) / obj$N
# Add each element in 'prior' to the corresponding obj$est$prior
for(k in seq_len(obj$num_groups)) {
obj$est[[k]]$prior <- prior[k]
}
class(obj) <- "dlda"
obj
}
.predict <- function(object, newdata, ...) { # Remove once sparsediscrim is reinstated to CRAN.
if (!inherits(object, "dlda")) {
stop("object not of class 'dlda'")
}
if (is.vector(newdata)) {
newdata <- as.matrix(newdata)
}
scores <- apply(newdata, 1, function(obs) {
sapply(object$est, function(class_est) {
with(class_est, sum((obs - xbar)^2 / object$var_pool) + log(prior))
})
})
if (is.vector(scores)) {
min_scores <- which.min(scores)
} else {
min_scores <- apply(scores, 2, which.min)
}
# Posterior probabilities via Bayes Theorem
means <- lapply(object$est, "[[", "xbar")
covs <- replicate(n=object$num_groups, object$var_pool, simplify=FALSE)
priors <- lapply(object$est, "[[", "prior")
posterior <- .posterior_probs(x=newdata,
means=means,
covs=covs,
priors=priors)
class <- factor(object$groups[min_scores], levels = object$groups)
list(class = class, scores = scores, posterior = posterior)
}
.posterior_probs <- function(x, means, covs, priors) { # Remove once sparsediscrim is reinstated to CRAN.
if (is.vector(x)) {
x <- matrix(x, nrow = 1)
}
x <- as.matrix(x)
posterior <- mapply(function(xbar_k, cov_k, prior_k) {
if (is.vector(cov_k)) {
post_k <- apply(x, 1, function(obs) {
.dmvnorm_diag(x=obs, mean=xbar_k, sigma=cov_k)
})
} else {
post_k <- dmvnorm(x=x, mean=xbar_k, sigma=cov_k)
}
prior_k * post_k
}, means, covs, priors)
if (is.vector(posterior)) {
posterior <- posterior / sum(posterior)
posterior <- matrix(posterior, nrow = 1) # Ensure it's always matrix, like just below.
colnames(posterior) <- names(priors)
} else {
posterior <- posterior / rowSums(posterior)
}
posterior
}
.dmvnorm_diag <- function(x, mean, sigma) { # Remove once sparsediscrim is reinstated to CRAN.
exp(sum(dnorm(x, mean=mean, sd=sqrt(sigma), log=TRUE)))
}
