R/wbyAgmv.R

In afthd: Accelerated Failure Time for High Dimensional Data with MCMC

#' @title
#'
#' Multivariate estimates of AFT model with weibull distribution using MCMC that supports augmented data.
#'
#' @description Provides estimate of AFT model including Survival time for augmented data with weibull distribution using
#' MCMC for multivariable (maximum 5 covariates of column at a time) in high dimensional gene
#' expression data. It also deals covariates with missing values.
#'
#' @details
#' Here weibull distribution has been used for AFT model with MCMC.
#' This function deals covariates (in data) with missing values. Missing value in any column (covariate) is replaced by mean of that particular covariate.
#'
#' @param m Starting column number of covariates of study in data.
#' @param n Ending column number of covariates of study in data.
#' @param p starting row number for augumented data in entered data.
#' @param q last row number for augumented data in entered data
#' @param t time (same unit as STime in data) after which, estimated STime to be printed (for individuals p to q).
#' @param STime name of survival time in data
#' @param Event name of event in data.  0 is for censored and 1 for occurrence of event.
#' @param nc number of markov chain.
#' @param ni number of iteration for MCMC.
#' @param data High dimensional gene expression data that contains event status, survival time and and set of covariates.
#' @return Posterior estimates of beta's, sigma , tau and deviance are their estimates mean, sd, credible intervals,number of efficient sample (n.eff) and  Rhat. beta's denotes posterior estimates of regression coefficient of the model. beta[1] is for intercept and others are for covariates (which is/are chosen as columns in data).'sigma' is the scale parameter of the distribution.
#' 'STime' in output section provides estimated value of  STime="os" in data only for individual row number p to q.
#' 'Overall_S' in output, provides an overall estimate of STime="os" in data for all individuals nrow(data).
#'  @import R2jags
#'
#' @references Prabhash et al(2016) <doi:10.21307/stattrans-2016-046>
#'
#' @examples
#' ##
#' data(hdata)
#' wbyAgmv(9,13,p=560,q=565,t=200,STime="os",Event="death",2,10,hdata)
#' #
#' ##
#' @export
#' @author Atanu Bhattacharjee, Gajendra Kumar Vishwakarma and Pragya Kumari
#' @seealso wbysmv
#'

wbyAgmv=function(m,n,p,q,t,STime,Event,nc,ni,data){

  nr<-nrow(data)

  if(STime!="os"){
    names(data)[names(data) == STime] <- "os"
  }
  if(Event!="death"){
    names(data)[names(data) == Event] <- "death"
  }

  d11 <- subset(data, select = c(m:n))
  le<-length(d11)
  for(i in 1:nr) {
    for(j in 1:le) {
      d11[i,j] = ifelse(is.na(d11[i,j])=="TRUE", mean(d11[,j], na.rm=TRUE), d11[i,j])
    }
  }

  pnt<-NULL
  for(j in 1:le)
  {
    if(sum(d11[,j])==0) {
      pnt<-c(pnt,j)
    }
  }

  if(is.null(pnt)==F){
    d11 <- d11[,-pnt]
  }

  len<-length(d11)

  d12<-data.frame(data[,c('death','os')],d11)

  mx<-max(d12$os) + 10
  surt<-ifelse(d12$death == 1, d12$os, NA)
  stcen<-ifelse(d12$death == 0, d12$os, mx)
  d12$os<-surt
  cen<-as.numeric(is.na(surt))
  d12<-data.frame(d12,stcen,cen)

  if(len>5){
    message("Outcome for first 5 covariates :")
    vv<-subset(d11,select = c(1:5))
  } else {
    vv<-d11
  }
  vname<-colnames(vv)

  if(len==1){
    data1<-list(os=d12$os, stcen=d12$stcen, cen=d12$cen, v1=vv[,1], N = nr, p=p, q=q, t=t)
    modelj1<-function(){
      for (i in 1:N) {
        sV1[i] <- (v1[i]-mean(v1[]))/sd(v1[])
        os[i] ~ dweib(alpha,lambda[i])
        cen[i] ~ dinterval(os[i],stcen[i])
        lambda[i] <-  log(2)*exp(-mu[i]*sqrt(tau))
        mu[i] <-  beta[1] + beta[2]*sV1[i]
        S[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))
      }
      alpha <- sqrt(tau)

      for(i in p:q) {
        #estimated survival S for augmented data or desired row number from p to q, 1<=p<=q<=N, where N = nrow(data)
        STime[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))
      }

      for(i in 1:2){
        beta[i] ~ dnorm(0,0.000001)
        rm[i] <- exp(beta[i])
        prob[i] <- step(beta[i])
      }
      tau ~ dgamma(0.001,0.001)
      sigma <- sqrt(1/tau)
    }

    inits1 <- function() {
      list(beta=c(0,0), tau=1)
    }
    jagsft <- jags(model.file=modelj1, data=data1, inits = inits1,
                   parameters.to.save = c('beta','tau','sigma','S','STime'), n.chains=nc, n.iter = ni)
    cat("Estimates for covariates : ", vname,"\n")
    f<-data.frame(jagsft$BUGSoutput$summary)
    ost<-list()
    for(i in 1:9){
      ost[i] <- mean(f[1:nr,i])
    }
    ost<-data.frame(ost)
    mtx<-matrix(nrow=7+q-p, ncol = 9)
    colnames(mtx)<-colnames(f)
    mtx<-data.frame(mtx)
    rownames(mtx)[1:(6+q-p)]<-rownames(f)[(1+nr):(nr+6+q-p)]
    rownames(mtx)[7+q-p]<-c("Overall_S")
    for(i in 1:(6+q-p)){
      mtx[i,]<-f[i+nr,]
    }
    mtx[7+q-p,]<-ost[1,]
    return(mtx)
  } else if(len==2){
    data2<-list(os=d12$os, stcen=d12$stcen, cen=d12$cen, v1=vv[,1], v2=vv[,2], N = nr , p=p, q=q, t=t)
    modelj2<-function(){
      for (i in 1:N) {
        sV1[i] <- (v1[i]-mean(v1[]))/sd(v1[])
        sV2[i] <- (v2[i]-mean(v2[]))/sd(v2[])
        os[i] ~ dweib(alpha,lambda[i])
        cen[i] ~ dinterval(os[i],stcen[i])
        lambda[i] <-  log(2)*exp(-mu[i]*sqrt(tau))
        mu[i] <-  beta[1] + beta[2]*sV1[i] + beta[3]*sV2[i]
        S[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))

      }
      alpha <- sqrt(tau)
      for(i in p:q) {
        #estimated survival S for augmented data or desired row number from p to q, 1<=p<=q<=N, where N = nrow(data)
        STime[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))
      }
      for(i in 1:3){
        beta[i] ~ dnorm(0,0.000001)
        rm[i] <- exp(beta[i])
        prob[i] <- step(beta[i])
      }
      tau ~ dgamma(0.001,0.001)
      sigma <- sqrt(1/tau)
    }

    inits2 <- function() {
      list(beta=c(0,0,0), tau=1)
    }
    jagsft <- jags(model.file=modelj2, data=data2, inits = inits2,
                   parameters.to.save = c('beta','tau','sigma','S','STime'), n.chains=nc, n.iter = ni)
    cat("Estimates for covariates : ", vname,"\n")
    f<-data.frame(jagsft$BUGSoutput$summary)
    ost<-list()
    for(i in 1:9){
      ost[i] <- mean(f[1:nr,i])
    }
    ost<-data.frame(ost)
    mtx<-matrix(nrow=8+q-p, ncol = 9)
    colnames(mtx)<-colnames(f)
    mtx<-data.frame(mtx)
    rownames(mtx)[1:(7+q-p)]<-rownames(f)[(1+nr):(nr+7+q-p)]
    rownames(mtx)[8+q-p]<-c("Overall_S")
    for(i in 1:(7+q-p)){
      mtx[i,]<-f[i+nr,]
    }
    mtx[8+q-p,]<-ost[1,]
    return(mtx)
  } else if(len==3){
    data3<-list(os=d12$os, stcen=d12$stcen, cen=d12$cen, v1=vv[,1], v2=vv[,2], v3=vv[,3], N = nr, p=p, q=q, t=t )
    modelj3<-function(){
      for (i in 1:N) {
        sV1[i] <- (v1[i]-mean(v1[]))/sd(v1[])
        sV2[i] <- (v2[i]-mean(v2[]))/sd(v2[])
        sV3[i] <- (v3[i]-mean(v3[]))/sd(v3[])
        os[i] ~ dweib(alpha,lambda[i])
        cen[i] ~ dinterval(os[i],stcen[i])
        lambda[i] <-  log(2)*exp(-mu[i]*sqrt(tau))
        mu[i] <-  beta[1] + beta[2]*sV1[i] + beta[3]*sV2[i] + beta[4]*sV3[i]
        S[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))
      }
      alpha <- sqrt(tau)
      for(i in p:q) {
        #estimated survival S for augmented data or desired row number from p to q, 1<=p<=q<=N, where N = nrow(data)
        STime[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))
      }
      for(i in 1:4){
        beta[i] ~ dnorm(0,0.000001)
        rm[i] <- exp(beta[i])
        prob[i] <- step(beta[i])
      }
      tau ~ dgamma(0.001,0.001)
      sigma <- sqrt(1/tau)
    }
    inits3 <- function() {
      list(beta=c(0,0,0,0), tau=1)
    }
    jagsft <- jags(model.file=modelj3, data=data3, inits = inits3,
                   parameters.to.save = c('beta','tau','sigma','S','STime'), n.chains=nc, n.iter = ni)
    cat("Estimates for covariates : ", vname,"\n")
    f<-data.frame(jagsft$BUGSoutput$summary)
    ost<-list()
    for(i in 1:9){
      ost[i] <- mean(f[1:nr,i])
    }
    ost<-data.frame(ost)
    mtx<-matrix(nrow=9+q-p, ncol = 9)
    colnames(mtx)<-colnames(f)
    mtx<-data.frame(mtx)
    rownames(mtx)[1:(8+q-p)]<-rownames(f)[(1+nr):(nr+8+q-p)]
    rownames(mtx)[9+q-p]<-c("Overall_S")
    for(i in 1:(8+q-p)){
      mtx[i,]<-f[i+nr,]
    }
    mtx[9+q-p,]<-ost[1,]
    return(mtx)
  } else if(len==4){
    data4<-list(os=d12$os, stcen=d12$stcen, cen=d12$cen, v1=vv[,1], v2=vv[,2], v3=vv[,3], v4=vv[,4], N = nr , p=p, q=q, t=t)
    modelj4<-function(){
      for (i in 1:N) {
        sV1[i] <- (v1[i]-mean(v1[]))/sd(v1[])
        sV2[i] <- (v2[i]-mean(v2[]))/sd(v2[])
        sV3[i] <- (v3[i]-mean(v3[]))/sd(v3[])
        sV4[i] <- (v4[i]-mean(v4[]))/sd(v4[])
        os[i] ~ dweib(alpha,lambda[i])
        cen[i] ~ dinterval(os[i],stcen[i])
        lambda[i] <-  log(2)*exp(-mu[i]*sqrt(tau))
        mu[i] <-  beta[1] + beta[2]*sV1[i] + beta[3]*sV2[i] + beta[4]*sV3[i]
        + beta[5]*sV4[i]
        S[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))

      }
      alpha <- sqrt(tau)
      for(i in p:q) {
        #estimated survival S for augmented data or desired row number from p to q, 1<=p<=q<=N, where N = nrow(data)
        STime[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))
      }
      for(i in 1:5){
        beta[i] ~ dnorm(0,0.000001)
        rm[i] <- exp(beta[i])
        prob[i] <- step(beta[i])
      }
      tau ~ dgamma(0.001,0.001)
      sigma <- sqrt(1/tau)
    }

    inits4 <- function() {
      list(beta=c(0,0,0,0,0), tau=1)
    }
    jagsft <- jags(model.file=modelj4, data=data4, inits = inits4,
                   parameters.to.save = c('beta','tau','sigma','S','STime'), n.chains=nc, n.iter = ni)
    message("Estimates using covariates : ", vname,"\n")
    f<-data.frame(jagsft$BUGSoutput$summary)
    ost<-list()
    for(i in 1:9){
      ost[i] <- mean(f[1:nr,i])
    }
    ost<-data.frame(ost)
    mtx<-matrix(nrow=10+q-p, ncol = 9)
    colnames(mtx)<-colnames(f)
    mtx<-data.frame(mtx)
    rownames(mtx)[1:(9+q-p)]<-rownames(f)[(1+nr):(nr+9+q-p)]
    rownames(mtx)[10+q-p]<-c("Overall_S")
    for(i in 1:(9+q-p)){
      mtx[i,]<-f[i+nr,]
    }
    mtx[10+q-p,]<-ost[1,]
    return(mtx)

  } else {
    data5<-list(os=d12$os, stcen=d12$stcen, cen=d12$cen, v1=vv[,1], v2=vv[,2], v3=vv[,3], v4=vv[,4], v5=vv[,5], N = nr, p=p, q=q,t=t)
    modelj5<-function(){
      for (i in 1:N) {
        sV1[i] <- (v1[i]-mean(v1[]))/sd(v1[])
        sV2[i] <- (v2[i]-mean(v2[]))/sd(v2[])
        sV3[i] <- (v3[i]-mean(v3[]))/sd(v3[])
        sV4[i] <- (v4[i]-mean(v4[]))/sd(v4[])
        sV5[i] <- (v5[i]-mean(v5[]))/sd(v5[])
        os[i] ~ dweib(alpha,lambda[i])
        cen[i] ~ dinterval(os[i],stcen[i])
        lambda[i] <- log(2)*exp(-mu[i]*sqrt(tau))
        mu[i] <-  beta[1] + beta[2]*sV1[i] + beta[3]*sV2[i] + beta[4]*sV3[i] + beta[5]*sV4[i] + beta[6]*sV5[i]
        S[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))
      }
      alpha <- sqrt(tau)
      for(i in p:q) {
        #estimated survival S for augmented data or desired row number from p to q, 1<=p<=q<=N, where N = nrow(data)
        STime[i] <- exp(-log(2)*exp(( log(t) - mu[i])*sqrt(tau)))
      }

      for(i in 1:6){
        beta[i] ~ dnorm(0,0.000001)
        rm[i] <- exp(beta[i])
        prob[i] <- step(beta[i])
      }
      tau ~ dgamma(0.001,0.001)
      sigma <- sqrt(1/tau)
    }

    inits5 <- function() {
      list(beta=c(0,0,0,0,0,0), tau=1)
    }
    jagsft <- jags(model.file=modelj5, data=data5, inits = inits5,
                   parameters.to.save = c('beta','tau','sigma','S','STime'), n.chains=2, n.iter = 10)
    cat("Estimates for covariates : ", vname,"\n")
    f<-data.frame(jagsft$BUGSoutput$summary)
    ost<-list()
    for(i in 1:9){
      ost[i] <- mean(f[1:nr,i])
    }
    ost<-data.frame(ost)
    mtx<-matrix(nrow=11+q-p, ncol = 9)
    colnames(mtx)<-colnames(f)
    mtx<-data.frame(mtx)
    rownames(mtx)[1:(10+q-p)]<-rownames(f)[(1+nr):(nr+10+q-p)]
    rownames(mtx)[11+q-p]<-c("OverallS")
    for(i in 1:(10+q-p)){
      mtx[i,]<-f[i+nr,]
    }
    mtx[11+q-p,]<-ost[1,]
    return(mtx)
  }

}


utils::globalVariables(c("N","v1","sd","v2","v3","v4","v5","tau","step"))