R/measurement_error_correction.R

In AFFECT: Accelerated Functional Failure Time Model with Error-Contaminated Survival Times

#' @title Correction of Measurement Error in Survival time and Censoring Status.
#'
#' @description This function aims to correct for measurement error in survival time and
#' misclassification in censoring status. The key strategy in the function \code{ME_correction} includes regression
#' calibration for survival time under additive measurement error models and the unbiased conditional expectation
#' approach for censoring status under misclassification models. With information of parameters in measurement error
#' models implemented, this function will give outputs with corrected survival time and censoring status.
#'
#' @param pi_01 Misclassifcation probability is P(Observed Censoring Status = 0| Actual Censoring Status = 1).
#'
#' @param pi_10 Misclassifcation probability is P(Observed Censoring Status = 1 | Actual Censoring Status = 0).
#'
#' @param gamma0 A scalar that links the observed survival time and true survival time in the classical additive
#'  measurement error model \code{y*=y+gamma0+gamma1*X+v}, where \code{y*}
#'  is observed survival time and \code{y} is true survival time, and \code{x} is covariates and \code{v} is noise
#'  term.
#' @param gamma1 A \code{p}-dimensional vector of parameters in the additive
#'  measurement error model \code{y*=y+gamma0+gamma1*X+v}, where \code{y*} is observed
#'   survival time and \code{y} is true survival time, \code{x} is covariates and \code{v} is
#'  noise term.
#'
#' @param indicator A \code{n}-dimensional vector of misclassified censoring status, such as the second column generated
#' by the function \code{gen_data}.
#'
#' @param yast A \code{n}-dimensional vector of error-prone survival time, such as the first column
#'generated by the function \code{gen_data}.
#'
#' @param cor_covar A \code{c(p,p)} covariance matrix of a \code{p}-dimensional vector of covariates.
#'
#' @param covariate A \code{c(n,p)} matrix of covariates.
#'
#' @return correction_data A \code{c(n,2)} data frame. This first column is the corrected survival time, and the
#'second column is the corrected censoring indicator.
#'
#' @examples
#' ## generate data with misclassification = 0.9 with n = 500,
#' ## p = 50 and variance of noise term is 0.75. The y* is related
#' ## to the first covariate.
#'
#' a <- matrix(0,ncol=50, nrow = 1);a[1,1] <- 1
#' data <- data_gen(n=500, p=50, pi_01 = 0.9, pi_10 = 0.9,
#' gamma0=1, gamma1=a, e_var=0.75)
#'
#' ## Assume that covariates are independent and
#' ## observed survival time is related to first covariate with
#' ## weight equals 1. And the scalar in the classical additive
#' ## measurement error model is 1 and is classifcation probability = 0.9.
#'
#' matrixa <- diag(50)
#' gamma_0 <-  1 ; gamma_1 <- matrix(0,ncol=50, nrow =1); gamma_1[1,1] <- 1
#' corrected_data1 <- ME_correction(pi_10=0.9,pi_01=0.9,gamma0 = gamma_0,
#' gamma1 = gamma_1,
#' cor_covar=matrixa, y=data[,1],
#' indicator=data[,2], covariate = data[,3:52])
#'
#' @export
#'

ME_correction <- function(pi_10,pi_01,gamma0,gamma1,cor_covar,indicator,yast,covariate){

  calculation<- function(x){
    values <- x-mean(x)
    return(values)
  }
  y = yast
  cor_covarw<- matrix(1,ncol=dim(covariate)[2], nrow = 1)
  matrixa <- cor_covarw %*% cor_covar
  co_mean <- as.data.frame(apply(covariate,2, calculation))
  correction_last_part <- c(matrixa%*%t(as.matrix(covariate)))
  mean_of_covariates <- apply(covariate,2,mean)
  estimated_w <- c(gamma0 + gamma1 %*% as.matrix(mean_of_covariates))
  y_hat <- (y -correction_last_part- estimated_w)
  y_hat <-data.frame(y_hat)
  colnames(y_hat) <- c("y_hat")


  indicator_hat_probability <- (indicator-pi_10)/(1-pi_10-pi_01)
  indicator_hat <- NULL
  for (i in c(1:length(indicator_hat_probability))){
    if (indicator_hat_probability[i]<0){
      indicator_hat[i] <- 0
    }else{
      indicator_hat[i] <- 1
    }
  }

  correction_data <- cbind(y_hat,indicator_hat)
  colnames(correction_data) <- c('corrected failure time', 'corrected censoring indicator')
  return(correction_data)
}