R/bayessurvreg3Para.R

################################################
#### AUTHOR:     Arnost Komarek             ####
####             (2005)                     ####
####                                        ####
#### FILE:       bayessurvreg3Para.R        ####
####                                        ####
#### FUNCTIONS:  bayessurvreg3Para          ####
################################################

###
###  This is a copy of bayessurvreg3 function from 20150326
###  where call to bayessurvreg.design function
###  is replaced by its copy to make this work in parallel
###  when running simulations.
###

### ======================================
### bayessurvreg3Para
### ======================================
bayessurvreg3Para <- function
(  formula,
   random,
   formula2,
   random2,
   data = parent.frame(),
   classification,
   classParam = list(Model = c("Examiner", "Factor:Examiner"),
                     a.sens = 1, b.sens = 1, a.spec = 1, b.spec = 1,
                     init.sens = NULL, init.spec = NULL),
   na.action = na.fail,
   onlyX = FALSE,
   nsimul = list(niter = 10, nthin = 1, nburn = 0, nwrite = 10),   
   prior,
   prior.beta,
   prior.b,
   init = list(iter = 0),
   mcmc.par = list(type.update.a = "slice",   k.overrelax.a = 1,   k.overrelax.sigma = 1,   k.overrelax.scale = 1,
                   type.update.a.b = "slice", k.overrelax.a.b = 1, k.overrelax.sigma.b = 1, k.overrelax.scale.b = 1),
   prior2,
   prior.beta2,
   prior.b2,
   init2,
   mcmc.par2 = list(type.update.a = "slice",   k.overrelax.a = 1,   k.overrelax.sigma = 1,   k.overrelax.scale = 1,
                    type.update.a.b = "slice", k.overrelax.a.b = 1, k.overrelax.sigma.b = 1, k.overrelax.scale.b = 1),
   priorinit.Nb,
   rho = list(type.update = "fixed.zero", init = 0, sigmaL = 0.1),
   store = list(a = FALSE, a2 = FALSE, y = FALSE, y2 = FALSE, r = FALSE, r2 = FALSE, b = FALSE, b2 = FALSE,
                a.b = FALSE, a.b2 = FALSE, r.b = FALSE, r.b2 = FALSE), 
   dir = getwd())
{  
  thispackage = "bayesSurv"
  #thispackage = NULL
  store  <- bayessurvreg3.checkStore(store)
  nsimul <- bayessurvreg.checknsimul(nsimul)
  
  transform = function(t){log(t)}
  dtransform = function(t){1/t}
  transform2 = function(t){t}
  dtransform2 = function(t){1}
  
  ## Give a function call to be recorded in a resulting object.
  call <- match.call(expand.dots = TRUE)
  
  ## Extract all the design information from the function call
  ## Do not transform the response at this moment
  doubly <- ifelse(missing(formula2), FALSE, TRUE)
  m <- match.call(expand.dots = FALSE)

  ##### ---------------------------------------------------------------------
  ##### bayessurvreg.design function
  #####  --> currently only for not doubly-censored data 
  ##### ---------------------------------------------------------------------
  ## des <- bayessurvreg.design(m = m, formula = formula, random = random, data = data, transform = transform2, dtransform = transform2)   ## This call replaced by the code below.

   tempF <- c("", "formula", "data", "subset", "na.action")
   mF <- m[match(tempF, names(m), nomatch=0)]
   mF[[1]] <- as.name("model.frame")
   special <- c("cluster")
   TermsF <- if(missing(data)) terms(formula, special)
             else              terms(formula, special, data = data)
   mF$formula <- TermsF
   mF <- eval(mF, parent.frame())
 
   ## Response matrix
   Y <- model.extract(mF, "response")
   Yinit <- Y
   if (!inherits(Y, "Surv"))
      stop("Response must be a survival object. ")
   type <- attr(Y, "type")
   if (type == 'counting') stop ("Invalid survival type ('counting' is not implemented). ")
   n <- nrow(Y)
   nY <- ncol(Y)
   
   ## Cluster indicators   
   cluster <- attr(TermsF, "specials")$cluster
   dropx <- NULL
   if (length(cluster)) {
     tempc <- untangle.specials(TermsF, "cluster", 1:10)
     ord <- attr(TermsF, "order")[tempc$terms]
     if (any(ord > 1))
         stop("Cluster can not be used in an interaction")
     cluster <- strata(mF[, tempc$vars], shortlabel = TRUE)
     dropx <- tempc$terms
   }
   if (length(dropx)){
     newTermsF <- TermsF[-dropx]
     attr(newTermsF, "intercept") <- attr(TermsF, "intercept")  ## I do not why but the command on the previous row
                                                                 ## sets attr(newTermsF, "intercept") always to 1,
                                                                 ## irrespective of what attr(TermsF, "intercept") was...
   }
   else{
     newTermsF <- TermsF
   }      

   ## Design matrix for both fixed and random effects X
   ## (finally, always without the intercept)
     ## Temporarily, include always at least the intercept
     ##  (to get nice rownames and initial estimates)  
   attr(newTermsF, "intercept") <- 1
   Xinit <- model.matrix(newTermsF, mF)
   rnamesX <- row.names(Xinit)
   cnamesX <- colnames(Xinit)
   
     ## Finally, intercept will be always removed
   nX <- ncol(Xinit) - 1
   cnamesX <- cnamesX[-1]
   if (nX){
     X <- Xinit
     X <- X[,-1]                                      ## removal of the intercept
     attr(X, "assign") <- attr(Xinit, "assign")[-1]   ## removal of the intercept
     attr(X, "contrasts") <- attr(Xinit, "contrasts")
   }
   else{
     X <- NULL
   }     
   indb <- if (nX) rep(-1, nX)                          ## initially, all effects are only fixed
           else    0
      
   ## Design matrix for random effects 
   randomInt <- FALSE
   nrandom <- 0
   if (!missing(random)){
     if (!length(cluster)) stop ("You have to indicate clusters when you wnat to include some random effects. ")
     tempR <- c("", "random", "data", "subset", "na.action")
     mR <- m[match(tempR, names(m), nomatch=0)]
     mR[[1]] <- as.name("model.frame")
     names(mR)[2] <- "formula"
     TermsR <- if(missing(data)) terms(random)
               else              terms(random, data = data)
     lTR <- length(attr(TermsR, "variables"))
     if (lTR == 1 & !attr(TermsR, "intercept")){        ## do nothing, in reality no random terms in the model
       names.random <- character(0)
     }
     else{
       if (lTR == 1 & attr(TermsR, "intercept")){        ## the only random term is the intercept
         randomInt <- TRUE
         names.random <- character(0)
       }
       else{
         mR$formula <- TermsR
         mR <- eval(mR, parent.frame())
         if (attr(TermsR, "intercept")){
           randomInt <- TRUE
           names.random <- colnames(model.matrix(TermsR, mR))[-1]
           ## attr(TermsR, "intercept") <- 0     ## remove it from the design matrix of random effects
         }
         else{
           names.random <- colnames(model.matrix(TermsR, mR))
         }           
       }
       nrandom <- 1*randomInt + length(names.random)
       if (sum(names.random %in% cnamesX) != nrandom - 1*randomInt) stop("Each random effect has to have also its fixed counterpart.")
       find.indeces <- function(all.eff){     ## Find indeces of random effects in a design matrix to be passed to C++
         where <- names.random %in% all.eff
         if (!sum(where)) return (-1)
         if (sum(where) > 1) stop("Error, contact the author.")
         index <- (1:length(names.random))[where]
         if (!randomInt) index <- index - 1
         return(index)
       }
       if (nX) indb <- as.numeric(apply(matrix(cnamesX, ncol = 1), 1, find.indeces))
     }       
   }
   else{
     names.random <- character(0)
   }     
   if (randomInt) names.random <- c("(Intercept)", names.random)
   nfixed <- nX - (nrandom - 1*randomInt)
   
   ## Give indeces of factors in the design matrix, it was used to define MH blocks in the earlier version of this program
   n.factors <- 0
   n.in.factors <- NULL
   factors <- NULL
   if (nX){
     temp <- attr(X, "assign")
     if (length(temp) == 1) factors <- 0
     else{
       factors <- numeric(length(temp))
       n.in.factors <- numeric(0)
       temp  <- temp - c(0, temp[1:(length(temp)-1)])
       i <- length(temp)
       while (i >= 1){
         if (temp[i] == 0){
           n.factors <- n.factors + 1
           factors[i] <- n.factors
           n.in.factor <- 1
           while (temp[i-1] == 0){
             i <- i - 1
             factors[i] <- n.factors
             n.in.factor <- n.in.factor + 1
           }
           i <- i - 1
           factors[i] <- n.factors
           n.in.factors <- c(n.in.factor + 1, n.in.factors)
         }
         else
           n.in.factors <- c(1, n.in.factors)         
         i <- i - 1
       }         
     }
     if (length(temp) != nX) stop("Something is wrong, contact the author.")
   }

   ## Sort everything according to the cluster indicator
   ##   and find the numbers of observations per cluster
   if (length(cluster)){
     ordering <- order(cluster)
     Y <- Y[ordering, ]
     cluster <- cluster[ordering]
     rnamesX <- rnamesX[ordering]
     if (nX){
       namesX <- cnamesX
       if (nX == 1) X <- matrix(X[ordering], ncol = 1)
       else         X <- as.matrix(X[ordering, ])
       colnames(X) <- namesX
     }
     ncluster <- length(attr(cluster, "levels"))
     helpf <- function(cl){return(sum(cluster %in% attr(cluster, "levels")[cl]))}
     nwithin <- apply(matrix(1:ncluster, ncol = 1), 1, "helpf")
   }
   else{
     if (nX){
       namesX <- cnamesX
       if (nX == 1) X <- matrix(X, ncol=1)
       else         X <-  as.matrix(X[,])
       colnames(X) <- namesX
     }
     cluster <- 1:n
     ncluster <- n
     nwithin <- rep(1, n)
   }     
   
   ## Transform the response
   if (type == 'interval') {
      if (any(Y[,3]==3)) Y <- cbind(eval(call("transform2", Y[,1:2])), Y[,3])
      else               Y <- cbind(eval(call("transform2", Y[,1])), Y[,3])
   }
   else if (type == 'left'){
          Y <- cbind(eval(call("transform2", Y[,1])), 2-Y[,2])   ## change 0 indicator into 2 indicating left censoring
        }
        else  ## type = 'right' or 'interval2'
           Y <- cbind(eval(call("transform2", Y[,1])), Y[,2])

   if (!all(is.finite(Y))) stop("Invalid survival times for this distribution (infinity on log-scale not allowed). ")

   des <- list(n = n, ncluster = ncluster, nwithin = nwithin, nY = nY, nX = nX,
               nfixed = nfixed, nrandom = nrandom, randomInt = randomInt,
               Y = Y, X = X, Yinit = Yinit, Xinit = Xinit, cluster = cluster, indb = indb,
               rnamesX = rnamesX, names.random = names.random, factors = factors, n.factors = n.factors, n.in.factors = n.in.factors)
  
  ##### ---------------------------------------------------------------------
  ##### End of bayessurvreg.design function
  ##### ---------------------------------------------------------------------
  
  if (doubly) des2 <- bayessurvreg.design(m = m, formula = formula2, random = random2, data = data, transform = transform2, dtransform = transform2)
  else        des2 <- list(nX = 0, n = des$n, nrandom = 0, randomInt = FALSE)
  
  if (onlyX){
    if (doubly) return(list(X = des$X, X2 = des2$X))
    else        return(des$X)
  } 

  if (!des$nrandom){
    store$b   <- FALSE
    store$a.b <- FALSE
    store$r.b <- FALSE
  }
  if (!des2$nrandom){
    store$b2   <- FALSE
    store$a.b2 <- FALSE
    store$r.b2 <- FALSE
  }  

  ## Perform some checks
  if (des$nrandom > 1) stop("Only intercept may appear in 'random'")
  if (des$nrandom & !des$randomInt) stop("Only intercept may appear in 'random'")
  
  if (des2$nrandom > 1) stop("Only intercept may appear in 'random2'")
  if (des2$nrandom & !des2$randomInt) stop("Only intercept may appear in 'random2'")

  nobs <- des$n
  if (doubly){
    nobs2 <- des2$n
    if (nobs != nobs2) stop("Inconsistent formula and formula2 (different number of observations indicated)")
  }

  
  ### ------------------------------------------------------------------------------------
  ### 201305: Code related to misclassification models
  ### ------------------------------------------------------------------------------------
  if (missing(classification)){
    mclass <- list(Model = "NONE", nModel = 0, nvisit = rep(1, nobs),
                   nExaminer = 0, labelExaminer = "0", Examiner = rep(0, nobs),
                   nFactor   = 0, labelFactor   = "0", Factor   = rep(0, nobs),
                   Prior     = rep(1, 4),
                   logTime   = 0, Status = 0,
                   init.sens = 0, init.spec = 0)

    ## mclass$nModel == 0 ==> no misclassification model 
    
  }else{
    if (doubly) stop("misclassification model not (yet) implemented for doubly-interval-censored data.")    
    mclass <- list()        ## To keep the information being then passed to C++
    
    if (!is.data.frame(classification)) stop("classification must be a data.frame")
    if (ncol(classification) < 4) stop("too low number of columns in the classification data.frame")    
    if (ncol(classification) == 4){
      colnames(classification) <- c("idUnit", "Time", "Examiner", "Status")
      classification[, "Factor"] <- 0
    }else{
      classification <- classification[, 1:5]
      colnames(classification) <- c("idUnit", "Time", "Examiner", "Status", "Factor")
    }  

    ### Some checks
    if (any(is.na(classification[, "idUnit"])))          stop("NA's not allowed in the idUnit column of classification.")    
    if (any(is.na(classification[, "Time"])))            stop("NA's not allowed in the Time column of classification.")
    if (any(classification[, "Time"] <= 0))              stop("all values in the Time column of classification must be positive.")
    if (any(is.na(classification[, "Examiner"])))        stop("NA's not allowed in the Examiner column of classification.")
    if (any(is.na(classification[, "Status"])))          stop("NA's not allowed in the Status column of classification.")
    if (any(!(classification[, "Status"] %in% c(0, 1)))) stop("all values in the Status column of classification must be 0/1.")
    
    ### idUnit values that appear in the data
    IDUnit <- unique(classification[, "idUnit"])                                             ### This maintains original ordering

    ### Consistency check    
    if (length(IDUnit) != nobs) stop("classification data.frame not consistent with the survival data (different number of experimental units)")

    ### Misclassification model
    classParam$Model <- classParam$Model[1]
    mclass$Model <- classParam$Model
    mclass$nModel <- pmatch(mclass$Model, table = c("Examiner", "Factor:Examiner"))
    if (is.na(mclass$nModel)) stop("unknown classParam$Model specified.")
        
    ### Number of visits per unit (tooth)
    mclass$nvisit <- as.numeric(table(factor(classification[, "idUnit"], levels = IDUnit)))  ### Also here, the original ordering is maintained

    ### Number of examinators
    mclass$nExaminer <- length(unique(classification[, "Examiner"]))

    ### Examiners labeled 0, ..., nExaminer - 1 and their original labels
    mclass$labelExaminer <- levels(factor(classification[, "Examiner"]))
    mclass$Examiner <- as.numeric(factor(classification[, "Examiner"])) - 1

    ### Factor related variables (if needed)
    switch (mclass$Model,
      "Examiner" = {
        mclass$nFactor <- 1

        mclass$labelFactor <- "0"
        mclass$Factor <- rep(0, sum(mclass$nvisit))
      },
      "Factor:Examiner" = {
        ### Number of unique values of Factor
        mclass$nFactor <- length(unique(classification[, "Factor"]))

        ### Factor values labeled 0, ..., nFactor - 1 and their original labels
        mclass$labelFactor <- levels(factor(classification[, "Factor"]))
        mclass$Factor <- as.numeric(factor(classification[, "Factor"])) - 1            
      },
      stop("some part of the code not yet implemented for this option.")             
    )
    
    ### Prior parameters
    if (is.null(classParam$a.sens)) classParam$a.sens <- 1
    if (is.null(classParam$b.sens)) classParam$b.sens <- 1    
    if (is.null(classParam$a.spec)) classParam$a.spec <- 1
    if (is.null(classParam$b.spec)) classParam$b.spec <- 1
    if (classParam$a.sens <= 0 | classParam$b.sens <= 0 | classParam$a.spec <= 0 | classParam$b.spec <= 0) stop("incorrect prior parameter for misclassification sensitivities/specificities")
    mclass$Prior <- c(classParam$a.sens, classParam$b.sens, classParam$a.spec, classParam$b.spec)
    names(mclass$Prior) <- c("a.sens", "b.sens", "a.spec", "b.spec")    

    ### Transformation of the visit times
    mclass$logTime <- log(classification[, "Time"])

    ### Status
    mclass$Status <- classification[, "Status"]
    
    ### Initial sensitivities
    if (is.null(classParam$init.sens)){
      switch (mclass$Model,
        "Examiner" = {
          classParam$init.sens <- runif(mclass$nExaminer, min = 0.80, max = 0.90)
          names(classParam$init.sens) <- mclass$labelExaminer
        },
        "Factor:Examiner" = {
          classParam$init.sens <- matrix(runif(mclass$nFactor * mclass$nExaminer, min = 0.80, max = 0.90), nrow = mclass$nFactor, ncol = mclass$nExaminer)
          rownames(classParam$init.sens) <- mclass$labelFactor
          colnames(classParam$init.sens) <- mclass$labelExaminer
        },
        stop("some part of the code not yet implemented for this option.") 
      )
    }
    if (length(classParam$init.sens) != mclass$nFactor * mclass$nExaminer) stop("incorrect classParam$init.sens specified")
    if (any(is.na(classParam$init.sens))) stop("NA's not allowed in classParam$init.sens")
    if (any(classParam$init.sens <= 0) | any(classParam$init.sens >= 1)) stop("all classParam$init.sens must lie in (0, 1)")

    mclass$init.sens <- as.numeric(classParam$init.sens)
    
    ### Initial specificities
    if (is.null(classParam$init.spec)){
      switch (mclass$Model,
        "Examiner" = {
          classParam$init.spec <- runif(mclass$nExaminer, min = 0.80, max = 0.90)
          names(classParam$init.spec) <- mclass$labelExaminer
        },
        "Factor:Examiner" = {
          classParam$init.spec <- matrix(runif(mclass$nFactor * mclass$nExaminer, min = 0.80, max = 0.90), nrow = mclass$nFactor, ncol = mclass$nExaminer)
          rownames(classParam$init.spec) <- mclass$labelFactor
          colnames(classParam$init.spec) <- mclass$labelExaminer
        },
        stop("some part of the code not yet implemented for this option.") 
      )
    }        
    if (length(classParam$init.spec) != mclass$nFactor * mclass$nExaminer) stop("incorrect classParam$init.spec specified")
    if (any(is.na(classParam$init.spec))) stop("NA's not allowed in classParam$init.spec")
    if (any(classParam$init.spec <= 0) | any(classParam$init.spec >= 1)) stop("all classParam$init.spec must lie in (0, 1)")

    mclass$init.spec <- as.numeric(classParam$init.spec)    
  }  
  ### ------------------------------------------------------------------------------------
    
  ## Priors and inits for beta parameters
  if (missing(init))        init <- list()
  if (missing(init2))       init2 <- list()
  if (missing(mcmc.par))    mcmc.par <- list()
  if (missing(mcmc.par2))   mcmc.par2 <- list()
  if (!doubly)              mcmc.par2 <- list()

  if (missing(prior.beta)) prior.beta <- list()
  betadi <- bayessurvreg3.priorBeta(prior.beta, init, des)
  init$beta <- attr(betadi, "init")
  prior.beta <- attr(betadi, "prior.beta")
  
  if (missing(prior.beta2)) prior.beta2 <- list()
  betadi2 <- bayessurvreg3.priorBeta(prior.beta2, init2, des2)
  init2$beta <- attr(betadi2, "init")
  prior.beta2 <- attr(betadi2, "prior.beta")
  
  ## Priors and inits for rho (correlation coefficient between two random intercepts)
  if (missing(priorinit.Nb)){
    rho <- bayessurvreg3.checkrho(rho=rho, doubly=doubly)
    if (rho$type.update == "fixed.zero") version <- 3
    else                                 version <- 31
  }
  else{
    rho <- list(type.update = "fixed.zero", init=0, sigmaL=0.1)
    rho <- bayessurvreg3.checkrho(rho=rho, doubly=doubly)    
    version <- 32
  }  
  
  ## Priors and inits for random intercept
  if (version %in% c(3, 31)){
    if (missing(prior.b)) prior.b <- list()  
    reffdi <- bayessurvreg3.priorb(prior.b=prior.b, init=init, design=des, mcmc.par=mcmc.par)
    prior.b  <- attr(reffdi, "prior.b")
    init     <- attr(reffdi, "init")
    mcmc.par <- attr(reffdi, "mcmc.par")
  
    ## Priors and inits for random intercept
    if (missing(prior.b2)) prior.b2 <- list()  
    reffdi2 <- bayessurvreg3.priorb(prior.b=prior.b2, init=init2, design=des2, mcmc.par=mcmc.par2)
    prior.b2  <- attr(reffdi2, "prior.b")
    init2     <- attr(reffdi2, "init")
    mcmc.par2 <- attr(reffdi2, "mcmc.par")
  }
  else{
    if (version == 32){
      Reff <- bayessurvreg3.priorinitNb(priorinit.Nb=priorinit.Nb, init=init, init2=init2, design=des, design2=des2, doubly=doubly)
      reffdi  <- Reff$reffdi
      reffdi2 <- Reff$reffdi2      
      priorinit.Nb <- attr(Reff, "priorinit.Nb")
      init         <- attr(Reff, "init")
      init2        <- attr(Reff, "init2")
      rm(list="Reff")
    }
    else{
      stop("It is strange but I cannot determine the version to be used")
    }  
  }  
  
  ## Priors and inits for error G-spline, (censored) observations and observational allocations
  if (!doubly) prior2 <- list()
  prinit <- bayessurvreg3.priorInit(prior, init, des, mcmc.par, prior2, init2, des2, mcmc.par2, doubly)
  init      <- attr(prinit, "init")
  prior     <- attr(prinit, "prior")  
  mcmc.par  <- attr(prinit, "mcmc.par")
  init2     <- attr(prinit, "init2")
  prior2    <- attr(prinit, "prior2")  
  mcmc.par2 <- attr(prinit, "mcmc.par2")


  ## Compute quantities to determine the space needed to be allocated
  ##   and numbers of iterations in different phases
  if (nsimul$nburn >= nsimul$niter) nsimul$nburn <- nsimul$niter - 1
  if (nsimul$nburn < 0) nsimul$nburn <- 0
 
  if (nsimul$nburn == 0) nruns <- 1
  else                   nruns <- 2

  nrun <- numeric(2)
  nrun[2] <- nsimul$niter - nsimul$nburn
  nrun[1] <- nsimul$nburn

  nwrite.run <- nrun
  nwrite.run[nsimul$nwrite <= nrun] <- nsimul$nwrite   
  max.nwrite <- max(nwrite.run)

  ## Write headers to files with stored values  
  bayessurvreg3.writeHeaders(dir = dir, doubly = doubly, prior.init = prinit,
                             priorb.di = reffdi, priorb2.di = reffdi2, store = store, design = des, design2 = des2, version = version,
                             mclass = mclass)  

  ## Combine similar parameters into one vector  
  dims <- c(nobs, as.numeric(doubly))
  storeV <- c(store$a, store$y, store$r, store$b, store$a2, store$y2, store$r2, store$b2, store$a.b, store$r.b, store$a.b2, store$r.b2)
  nsimul.run1 <- c(nrun[1], nsimul$nthin, nwrite.run[1])
  nsimul.run2 <- c(nrun[2], nsimul$nthin, nwrite.run[2])
  names(nsimul.run1) <- names(nsimul.run2) <- c("niter", "nthin", "nwrite")
  nsample1 <- nsimul.run1["niter"] %/% nsimul.run1["nthin"]
  nsample2 <- nsimul.run2["niter"] %/% nsimul.run2["nthin"]  
  
  cat("Simulation started on                       ", date(), "\n", sep = "")
  fit <- .C(C_bayessurvreg2, as.character(dir),
                             dims = as.integer(dims),
                             X1 = as.double(if(des$nX) t(des$X) else 0),
                             X2 = as.double(if(des2$nX) t(des2$X) else 0),
                             y1.left = as.double(prinit$y.left),
                             y1.right = as.double(prinit$y.right),
                             status1 = as.integer(prinit$status),
                             t2.left = as.double(prinit$t2.left),
                             t2.right = as.double(prinit$t2.right),
                             status2 = as.integer(prinit$status2),
                             iPML = double(dims[1]),
                             Ys1 = as.double(prinit$y),
                             Ys2 = as.double(prinit$y2),
                             r1 = as.integer(prinit$r),
                             r2 = as.integer(prinit$r2),
                             specif = as.integer(prinit$specification),
                             r1.b = as.integer(reffdi$r),
                             r2.b = as.integer(reffdi2$r),
                             specif.b = as.integer(c(reffdi$specification, reffdi2$specification)),
                             GsplI1 = as.integer(prinit$Gparmi),
                             GsplD1 = as.double(prinit$Gparmd),
                             GsplI2 = as.integer(prinit$Gparmi2),
                             GsplD2 = as.double(prinit$Gparmd2),
                             priorBetaI1 = as.integer(betadi$parmI),
                             priorBetaD1 = as.double(betadi$parmD),
                             priorBetaI2 = as.integer(betadi2$parmI),
                             priorBetaD2 = as.double(betadi2$parmD),
                             priorbI1 = as.integer(reffdi$bparmI),
                             priorbD1 = as.double(reffdi$bparmD),
                             priorbI2 = as.integer(reffdi2$bparmI),
                             priorbD2 = as.double(reffdi2$bparmD),
                             priorCovMatI1 = as.integer(reffdi$GsplI),
                             priorCovMatD1 = as.double(reffdi$GsplD),
                             priorCovMatI2 = as.integer(reffdi2$GsplI),
                             priorCovMatD2 = as.double(reffdi2$GsplD),
                             rho             = as.double(rho$init),
                             rho.accept      = integer(nsample1),
                             rho.type.update = as.integer(rho$typeI),
                             rho.sigmaL      = as.double(rho$sigmaL),
                             mclass.sens.spec = as.double(c(mclass$init.sens, mclass$init.spec)),            
                             mclass.logtime   = as.double(mclass$logTime),
                             mclass.status    = as.integer(mclass$Status),
                             mclass.paramI    = as.integer(c(mclass$nModel, mclass$nExaminer, mclass$nFactor, mclass$nvisit, mclass$Examiner, mclass$Factor)),
                             mclass.paramD    = as.double(mclass$Prior),
                             iter = as.integer(prinit$iter),
                             nsimul = as.integer(nsimul.run1),
                             store = as.integer(storeV),
                             version = as.integer(version),
                             mainSimul = as.integer(0),            
                             err = integer(1),
           PACKAGE = thispackage)
  
  if (fit$err != 0) stop ("Something went wrong during the simulation.")
  cat("Burn-up finished on                         ", date(), "   (iteration ", fit$iter, ")", "\n", sep = "")
  
  ## Rewrite sampled values by new files
  bayessurvreg3.writeHeaders(dir = dir, doubly = doubly, prior.init = prinit,
                             priorb.di = reffdi, priorb2.di = reffdi2, store = store, design = des, design2 = des2, version = version,
                             mclass = mclass)  
  
  ## Main simulation
  fit <- .C(C_bayessurvreg2, as.character(dir),
                             dims = as.integer(dims),
                             X1 = as.double(fit$X1),
                             X2 = as.double(fit$X2),
                             y1.left = as.double(fit$y1.left),
                             y1.right = as.double(fit$y1.right),
                             status1 = as.integer(fit$status1),
                             t2.left = as.double(fit$t2.left),
                             t2.right = as.double(fit$t2.right),
                             status2 = as.integer(fit$status2),
                             iPML = double(dims[1]),            
                             Ys1 = as.double(fit$Ys1),
                             Ys2 = as.double(fit$Ys2),
                             r1 = as.integer(fit$r1),
                             r2 = as.integer(fit$r2),
                             specif = as.integer(fit$specif),
                             r1.b = as.integer(fit$r1.b),
                             r2.b = as.integer(fit$r2.b),
                             specif.b = as.integer(fit$specif.b),            
                             GsplI1 = as.integer(fit$GsplI1),
                             GsplD1 = as.double(fit$GsplD1),
                             GsplI2 = as.integer(fit$GsplI2),
                             GsplD2 = as.double(fit$GsplD2),
                             priorBetaI1 = as.integer(fit$priorBetaI1),
                             priorBetaD1 = as.double(fit$priorBetaD1),
                             priorBetaI2 = as.integer(fit$priorBetaI2),
                             priorBetaD2 = as.double(fit$priorBetaD2),
                             priorbI1 = as.integer(fit$priorbI1),
                             priorbD1 = as.double(fit$priorbD1),
                             priorbI2 = as.integer(fit$priorbI2),
                             priorbD2 = as.double(fit$priorbD2),
                             priorCovMatI1 = as.integer(fit$priorCovMatI1),
                             priorCovMatD1 = as.double(fit$priorCovMatD1),
                             priorCovMatI2 = as.integer(fit$priorCovMatI2),
                             priorCovMatD2 = as.double(fit$priorCovMatD2),
                             rho             = as.double(fit$rho),
                             rho.accept      = integer(nsample2),
                             rho.type.update = as.integer(fit$rho.type.update),
                             rho.sigmaL      = as.double(fit$rho.sigmaL),
                             mclass.sens.spec = as.double(fit$mclass.sens.spec),            
                             mclass.logtime   = as.double(mclass$logTime),
                             mclass.status    = as.integer(mclass$Status),
                             mclass.paramI    = as.integer(c(mclass$nModel, mclass$nExaminer, mclass$nFactor, mclass$nvisit, mclass$Examiner, mclass$Factor)),
                             mclass.paramD    = as.double(mclass$Prior),            
                             iter = as.integer(fit$iter),
                             nsimul = as.integer(nsimul.run2),
                             store = as.integer(fit$store),
                             version = as.integer(fit$version),
                             mainSimul = as.integer(1),
                             err = integer(1),
           PACKAGE = thispackage)
  
  if (fit$err != 0) stop ("Something went wrong during the simulation.")  
  cat("Simulation finished on                      ", date(), "   (iteration ", fit$iter, ")", "\n", sep = "")     
  
  toreturn <- fit$iter
  attr(toreturn, "call") <- call

  attr(toreturn, "init") <- init
  attr(toreturn, "prior") <- prior
  attr(toreturn, "prior.beta") <- prior.beta
  attr(toreturn, "mcmc.par") <- mcmc.par

  if (doubly){
    attr(toreturn, "init2") <- init2
    attr(toreturn, "prior2") <- prior2
    attr(toreturn, "prior.beta2") <- prior.beta2
    attr(toreturn, "mcmc.par2") <- mcmc.par2
  }

  if (version %in% c(3, 31)){
    attr(toreturn, "prior.b") <- prior.b
    if (doubly) attr(toreturn, "prior.b2") <- prior.b2
  }  
  if (version == 31){
    attr(toreturn, "rho.accept") <- fit$rho.accept
  }  

  if (version == 32){
    attr(toreturn, "priorinit.Nb") <- fit$priorinit.Nb
  }

  if (mclass$nModel){
    attr(toreturn, "prior.classification") <- mclass$Prior
    attr(toreturn, "init.sens") <- classParam$init.sens
    attr(toreturn, "init.spec") <- classParam$init.spec

    ##### Calculate model fit statistics
    ##### -----------------------------------
    
    ### Deviance sample
    devGJK <- scan(paste(dir, "/devianceGJK.sim", sep = ""), skip = 1, quiet = TRUE)
    D.bar <- mean(devGJK, na.rm = TRUE)

    ### Posterior means of sensitivities/specificities
    sensspec <- matrix(scan(paste(dir, "/sens_spec.sim", sep = ""), skip = 1, quiet = TRUE), ncol = 2 * mclass$nExaminer * mclass$nFactor, byrow = TRUE)
    colnames(sensspec) <- scan(paste(dir, "/sens_spec.sim", sep = ""), what = character(), nlines = 1, quiet = TRUE)
    sensspec.bar <- apply(sensspec, 2, mean, na.rm = TRUE)
    rm(list = "sensspec")

    ### Init for posterior means of linear predictors
    eta.bar <- rep(0, des$n)
    
    ### Posterior means of random effects
    if (des$nrandom > 0){
      if (store$b){        
        bb <- matrix(scan(paste(dir, "/b.sim", sep = ""), skip = 1, quiet = TRUE), ncol = des$ncluster, byrow = TRUE)
        bb.bar <- apply(bb, 2, mean, na.rm = TRUE)
        rm(list = "bb")
      }else{
        bb.bar <- rep(NA, des$ncluster)
        warning("to have DIC calculated, store$b must be TRUE.")
      }
      eta.bar <- eta.bar + rep(bb.bar, des$nwithin)
    }

    ### Posterior means of fixed effects
    if (des$nX){
      beta <- matrix(scan(paste(dir, "/beta.sim", sep = ""), skip = 1, quiet = TRUE), ncol = des$nX, byrow = TRUE)
      colnames(beta) <- scan(paste(dir, "/beta.sim", sep = ""), what = character(), nlines = 1, quiet = TRUE)
      beta.bar <- apply(beta, 2, mean, na.rm = TRUE)
      rm(list = "beta")
      eta.bar <- eta.bar + as.numeric(des$X %*% beta.bar)
    }  
    
    ### Posterior means of G-spline (error term) parameters
    gspline <- matrix(scan(paste(dir, "/gspline.sim", sep = ""), skip = 1, quiet = TRUE), ncol = 5, byrow = TRUE)
    colnames(gspline) <- scan(paste(dir, "/gspline.sim", sep = ""), what = character(), nlines = 1, quiet = TRUE)
    gspline.bar <- apply(gspline, 2, mean, na.rm = TRUE)
    gspline.sigmaSq.bar <- mean(gspline[, "sigma1"]^2, na.rm = TRUE)
    gspline.scaleSq.bar <- mean(gspline[, "scale1"]^2, na.rm = TRUE)    
    rm(list = "gspline")

    mweight <- matrix(scan(paste(dir, "/mweight.sim", sep = ""), skip = 1, quiet = TRUE), ncol = 2 * prinit$Gparmi["K1"] + 1, byrow = TRUE)
    mweight.bar <- apply(mweight, 2, mean, na.rm = TRUE)
    rm(list = "mweight")
    
    ### Input for D.in.bar to be calculated in C++
    if (any(is.na(eta.bar))){
      D.in.bar <- NA      
    }else{
      maxnvisit <- max(mclass$nvisit)
      nss <- mclass$nExaminer * mclass$nFactor
      inD <- .C(C_iPML_misclass_GJK,     
                iPML  = double(des$n),
                dwork = double(3 * (1 + maxnvisit)),
                min_etaM  = as.double((-1) * eta.bar),
                sens      = as.double(sensspec.bar[1:nss]),
                spec      = as.double(sensspec.bar[(nss + 1):(2 * nss)]),
                logvtime  = as.double(mclass$logTime),
                status    = as.integer(mclass$Status),
                nExaminer = as.integer(mclass$nExaminer),
                nFactor   = as.integer(mclass$nFactor),
                nvisit    = as.integer(mclass$nvisit),
                maxnvisit = as.integer(maxnvisit),
                Examiner  = as.integer(mclass$Examiner),
                Factor    = as.integer(mclass$Factor),
                gg_K      = as.integer(prinit$Gparmi["K1"]),
                gg_gamma  = as.double(gspline.bar["gamma1"]),
                gg_delta  = as.double(gspline.bar["delta1"]),
                gg_sigma  = as.double(sqrt(gspline.sigmaSq.bar)),
                #gg_sigma  = as.double(gspline.bar["sigma1"]),
                gg_intcpt = as.double(gspline.bar["intercept1"]),
                gg_scale  = as.double(sqrt(gspline.scaleSq.bar)),
                #gg_scale  = as.double(gspline.bar["scale1"]),
                gg_w      = as.double(mweight.bar),
                nP        = as.integer(des$n),
                PACKAGE = thispackage)
      D.in.bar <- -2 * sum(log(inD$iPML))      
    }
           
    ### DIC and log(PML)
    pD <- D.bar - D.in.bar
    DIC <- D.bar + pD
    logPML <- sum(log(fit$iPML))
    
    attr(toreturn, "fitStat") <- c(logPML, DIC, D.bar, D.in.bar, pD)
    names(attr(toreturn, "fitStat")) <- c("logPML", "DIC", "D.bar", "D.in.bar", "pD")
    attr(toreturn, "iPML")   <- fit$iPML
    
  }
  
  class(toreturn) <- "bayessurvreg3"
  return(toreturn)    
}

Try the bayesSurv package in your browser

Any scripts or data that you put into this service are public.

bayesSurv documentation built on May 2, 2019, 3:26 a.m.