R/santaR_CBand.R

In adwolfer/santaR: Short Asynchronous Time-Series Analysis

#' Compute Group Mean Curve Confidence Bands
#'
#' Generate bootstrapped group mean curve Confidence Bands, by resampling of individual curves with replacement. Returns a \emph{SANTAObj} with added Confidence Bands.
#' \itemize{
#'   \item Resampling whole data curves assumes less of the data than resampling of residuals.
#'   \item The resampled distribution is of same size as the original distribution (same number of individuals in each group as in the input data).
#'   \item The degree of freedom for the estimator is identical to the one employed for curve fitting in \code{\link{santaR_fit}}.
#' }
#'
#' @param SANTAObj A fitted \emph{SANTAObj} as generated by \code{\link{santaR_fit}}.
#' @param nBoot (int) Number of bootstrapping rounds. Default 1000.
#' @param alpha (float) Confidence \emph{(0.05 for 95\% Confidence Bands)}. Default 0.05.
#' @param subsampling (int) Number of points to sample in the time range (for the estimator and Confidence Bands). Default is 250.
#'
#' @return A \emph{SANTAObj} with added Confidence Bands for each group.
#'
#' @examples
#' ## 56 measurements, 8 subjects, 7 unique time-points
#' ## Default parameter values decreased to ensure an execution < 2 seconds
#' Yi          <- acuteInflammation$data$var_3
#' ind         <- acuteInflammation$meta$ind
#' time        <- acuteInflammation$meta$time
#' group       <- acuteInflammation$meta$group
#' grouping    <- get_grouping(ind, group)
#' inputMatrix <- get_ind_time_matrix(Yi, ind, time)
#' SANTAObj    <- santaR_fit(inputMatrix, df=5, grouping=grouping, verbose=TRUE)
#' SANTAObj    <- santaR_CBand(SANTAObj, nBoot=100)
#'
#' @family Analysis
#'
#' @export
santaR_CBand              <- function(SANTAObj,nBoot=1000,alpha=0.05,subsampling=250) {  
  ## FUNCTION
  sp.resampler        <- function(sp.frame) {
    n             <- nrow(sp.frame)
    resample.rows <- sample(1:n, size=n, replace=TRUE)  # which rows to choose
    return(sp.frame[resample.rows,])
  }                         # return a resampled sp.frame
  sp.spline.estimator <- function(sp.frame,df,time,grid) {
    
    ## INPUT
    # sp.frame    = the IND x TIME dataset to estimate (data.frame)
    # df          = degree of freedom previously used to fit the dataset (reused for estimation)
    # time        = vector of time used in sp.frame
    # grid        = projection grid
    
    ## Group mean
    mean.val    <- colMeans(sp.frame, na.rm=TRUE)
    notNA       <- !is.na(mean.val)
    
    ## Fit spline on group mean
    if( sum(notNA) >= df ) {
      mean.fit  <- stats::smooth.spline( x=time[notNA], y=mean.val[notNA], df=sum(notNA) )
    } else {
      message("Exit: number of TP < df")
      stop
    }
    ## Return prediction on the grid    
    return( stats::predict( mean.fit, grid )$y )
  }            # return the predicted mean spline over subsampling points on range
  
  ## Initialisation  
  df    <-  SANTAObj$properties$df
  
  ## Bootstrap for each group
  for(nGroup in 1: length(SANTAObj$groups)) {
    if( is.data.frame(SANTAObj$groups[[nGroup]]$groupData.pred) ) {         ## check there is min 1 ind
      # Init
      sp.frame  <- SANTAObj$groups[[nGroup]]$groupData.pred
      time      <- as.numeric(colnames(sp.frame))
      rng       <- range(time)
      grid      <- seq( from=rng[1], to=rng[2], length.out=subsampling )
      
      # Background spline - needed for other version
      #spline.main <- sp.spline.estimator(sp.frame,df,time,grid)
      
      # nBoot bootstrap
      spline.boot <- replicate( n=nBoot, sp.spline.estimator( sp.resampler(sp.frame), df, time, grid)) # nrow=length(grid) ncol=nBoot
      #cband.lower   <- 2*spline.main - apply(spline.boot,1,quantile,probs=1-alpha/2)  # pivotal confidence intervals
      #cband.upper   <- 2*spline.main - apply(spline.boot,1,quantile,probs=alpha/2)    # pivotal confidence intervals
      cband.upper    <- apply(spline.boot,1,stats::quantile,probs=1-alpha/2)  #95% confidence limits -> the 2.5th and 97.5th centiles
      cband.lower    <- apply(spline.boot,1,stats::quantile,probs=alpha/2)
      SANTAObj$groups[[nGroup]]$groupCBand$upperFit   <- stats::smooth.spline(x=grid, y=cband.upper,df=length(grid))
      SANTAObj$groups[[nGroup]]$groupCBand$lowerFit   <- stats::smooth.spline(x=grid, y=cband.lower,df=length(grid))
    } else {
      SANTAObj$groups[[nGroup]]$groupCBand$upperFit   <- NA
      SANTAObj$groups[[nGroup]]$groupCBand$lowerFit   <- NA
    }
  }  
  SANTAObj$properties$CBand$status <- TRUE
  SANTAObj$properties$CBand$alpha  <- alpha
	SANTAObj$properties$CBand$nBoot  <- nBoot
  return(SANTAObj)
}