R/zeitzeiger_group.R In jakejh/zeitzeiger: Regularized supervised learning for high-dimensional data from an oscillatory system

Documented in zeitzeigerPredictGroup

fxGroup = function(x, timeDiff, time, xFitMean, xFitResid, knots, logArg = FALSE) {
# dim(x): c(m, p)
# length(timeDiff): m
# length(time): n
timeRep = rowSums(expand.grid(time, timeDiff)) %% 1
xRep = x[rep(1:nrow(x), each=length(time)), , drop = FALSE]
logLike = matrix(data = NA, nrow = length(time), ncol = ncol(x))

for (jj in 1:ncol(x)) {
xPredMean = predictIntensity(xFitMean[jj, , drop = FALSE], timeRep,
knots = knots)[, 1]
logLikeTmp = stats::dnorm(xRep[, jj], mean = xPredMean, sd = xFitResid[jj],
log = TRUE)
logLike[, jj] = rowSums(matrix(logLikeTmp, nrow = length(time)))}

if (logArg) {
return(logLike)
} else {
return(exp(logLike))}}

zeitzeigerPredictGivenDensityGroup = function(xTest, groupTest, xFitMean,
xFitResid, knots, timeRange) {
groups = groupTest$group groupsUnique = sort(unique(groups)) negLogLike = matrix(NA, nrow = length(groupsUnique), ncol = length(timeRange)) mleFit = list() for (ii in 1:length(groupsUnique)) { xTestNow = xTest[groups == groupsUnique[ii], , drop=FALSE] timeDiffNow = groupTest[groups == groupsUnique[ii], ]$timeDiff
negLogLike[ii, ] = -fxGroup(xTestNow, timeDiffNow, timeRange, xFitMean,
xFitResid, knots, logArg = TRUE)

timeStart = timeRange[which.min(negLogLike[ii, ])]
names(timeStart) = 'time'
negLogLikeFunc = function(time) -sum(fxGroup(xTestNow, timeDiffNow, time,
xFitMean, xFitResid, knots,
logArg = TRUE))
bbmle::parnames(negLogLikeFunc) = names(timeStart)
warnOrig = getOption('warn')
options(warn = -1)
mleFit[[ii]] = bbmle::mle2(negLogLikeFunc, start = timeStart, vecpar = TRUE,
method = 'L-BFGS-B', lower = 0, upper = 1)
options(warn = warnOrig)}

timePred = sapply(mleFit, function(x) x@coef)
return(list(timeDepLike = exp(-negLogLike), mleFit = mleFit, timePred = timePred))}

#' Predict corresponding time for groups of test observations
#'
#' zeitzeigerPredictGroup predicts the value of the periodic variable
#' for each group of test observations, where the amount of time between each
#' observation in a group is known. This function calls bbmle::mle2.
#'
#' @param xTrain Matrix of measurements for training data, observations in rows
#' and features in columns.
#' @param timeTrain Vector of values of the periodic variable for training observations,
#' where 0 corresponds to the lowest possible value and 1 corresponds to the highest
#' possible value.
#' @param xTest Matrix of measurements for test data, observations in rows
#' and features in columns.
#' @param groupTest data.frame with one row per observation in xTest, and columns
#' for group and timeDiff. Observations in the same group should have the
#' same value of group. Within each group, the value of timeDiff
#' should correspond to the amount of time between that observation and a reference time.
#' Typically, timeDiff will equal zero for one observation per group.
#' @param spcResult Result from zeitzeigerSpc.
#' @param nKnots Number of internal knots to use for the periodic smoothing spline.
#' @param nSpc Vector of the number of SPCs to use for prediction. If NA (default),
#' nSpc will become 1:K, where K is the number of SPCs in spcResult.
#' Each value in nSpc will correspond to one prediction for each test observation.
#' A value of 2 means that the prediction will be based on the first 2 SPCs.
#' @param timeRange Vector of values of the periodic variable at which to calculate likelihood.
#' The time with the highest likelihood is used as the initial value for the
#' MLE optimizer.
#'
#' @return A list with the following elements, where the groups will be sorted by their names.
#' \item{timeDepLike}{3-D array of likelihood, with dimensions for each group of test observations,
#' each element of nSpc, and each element of timeRange.}
#' \item{mleFit}{List (for each element in nSpc) of lists (for each group of test observations)
#' of mle2 objects.}
#' \item{timePred}{Matrix of predicted times for each group of test observations by values of nSpc.}
#'
#' @seealso \link{zeitzeigerPredict}
#'
#' @export
zeitzeigerPredictGroup = function(xTrain, timeTrain, xTest, groupTest, spcResult,
nKnots = 3, nSpc = NA,
timeRange = seq(0, 1 - 0.01, 0.01)) {
zTrain = xTrain %*% spcResult$v zTest = xTest %*% spcResult$v

zFitResult = zeitzeigerFit(zTrain, timeTrain, nKnots)
zFitMean = zFitResult$xFitMean zFitResid = zFitResult$xFitResid

knots = seq(0, 1 - 1 / nKnots, length = nKnots)

if (length(nSpc)==1 && is.na(nSpc)) {
nSpc = 1:length(spcResult$d) } else { if (!all(nSpc %in% 1:length(spcResult$d)) || anyDuplicated(nSpc) > 0) {
stop('nSpc must be unique integers in 1 to the number of singular vectors.')}}

timePred = matrix(NA, nrow = length(unique(groupTest$group)), length(nSpc)) timeDepLike = array(NA, dim = c(length(unique(groupTest$group)), length(nSpc),
length(timeRange)))
mleFit = list() # list (for each nSpc) of lists (for each group)
for (ii in 1:length(nSpc)) {
predResult = zeitzeigerPredictGivenDensityGroup(zTest[, 1:nSpc[ii], drop = FALSE],
groupTest,
zFitMean[1:nSpc[ii], , drop = FALSE],
zFitResid[1:nSpc[ii]],
knots, timeRange)
timeDepLike[, ii, ] = predResult$timeDepLike mleFit[[ii]] = predResult$mleFit
timePred[, ii] = predResult\$timePred}
return(list(timeDepLike = timeDepLike, mleFit = mleFit, timePred = timePred))}

jakejh/zeitzeiger documentation built on Nov. 22, 2018, 6:53 a.m.