R/ICA.ContCont.MultS.MPC.R

In Surrogate: Evaluation of Surrogate Endpoints in Clinical Trials

ICA.ContCont.MultS.MPC = function (M = 1000, N, Sigma,prob=NULL, Seed = 123, Save.Corr=F,
                                      Show.Progress = FALSE) 
{

  d = nrow(Sigma)[1]
  Vars = diag(Sigma)
  IND = ks::vec(matrix(1:d, ncol = 2, byrow = T), byrow = F)
  Sigma = Sigma[IND, IND]
  R = cov2cor(Sigma)
  if (!.is.PD(R)) {
    alpha = uniroot(function(alpha, R, Rfixed, tol = 1e-04) {
      if (anyNA(R)) {
        R[is.na(R)] = 0
      }
      f = alpha * R + (1 - alpha) * Rfixed
      min(eigen(f)$values) - tol
    }, c(0, 1), R = R, Rfixed = diag(d), tol = 1e-08)$root
    R = R * alpha + (1 - alpha) * diag(d)
    warning(paste("The initial correlation matrix is not PD. TThe matrix was shrunk by a factor alpha=", 
                  alpha, " for correction", sep = ""))
  }
  p = (d-2)/2
  
  if(is.null(prob)){
    prob = choose(p,1:p)/sum(choose(p,1:p))
  }
  
  if(Show.Progress==F){
    opb <- pbapply::pboptions(type="none")
    on.exit(pbapply::pboptions(opb))
  }
  
  Results = pbapply::pblapply(X=1:M,function(X) {
    colnames(R) = rownames(R)=c('T0',paste('S',1:p,0,sep=''),'T1',paste('S',1:p,1,sep=''))
    sum.r = 0
    while(sum.r==0){
      r.num = sample(1:p,1,prob=prob)
      r = sample(c(rep(1,r.num),rep(0,p-r.num))) 
      sum.r=sum(r)
    }
    
    if(sum(r) < p){
      var.p = (1:p)[r == 1] 
      fixZero.p = (1:p)[r == 0]
      colr = 2+p + (1:p)[r==0]
      rowr = (1:p)[r==0] + 1
      R[1:(p+1),colr] = R[colr,1:(p+1)] = 0
      R[rowr,(2+p):d] = R[(2+p):d,rowr] = 0
      
      if(length(var.p)==1){
        if(length(fixZero.p)==1){
          p.index = c(var.p,fixZero.p)
        }else{
          p.index = c(var.p,sample(fixZero.p))        
        }
      }else{
        if(length(fixZero.p)==1){
          p.index = c(sample(var.p),fixZero.p)
        }else{
          p.index = c(sample(var.p),sample(fixZero.p))    
        }
      }
    }else{
      p.index=1:p
    }
    
    IND.2 =  c(1,1+p.index,p+2+p.index[p:1],p+2)
    R.test = R[IND.2,IND.2]
    parm.a = sum(r)+1
    .r.random = tryCatch(.Correlation.matrix.PC(R.test,Range = c(-1, 1),parm.a),error=function(e){NULL})
    if(is.null(.r.random)){
      .r.random = matrix(NA,d,d)
      return(c(NA,NA, r,.r.random[lower.tri(.r.random)]))
    }
    .r.random = .r.random[order(IND.2),order(IND.2)]
    IND = ks::vec(matrix(1:d, ncol = 2), byrow = TRUE)
    .r.random = .r.random[IND,IND]
    ICA = .MultivarICA.fun(.r.random, Vars, N)
    if(Save.Corr){
      return(c(ICA,r, .r.random[lower.tri(.r.random)]))      
    }else{
      return(c(ICA,r))      
    }
    
  }, cl=NULL)
  Results = do.call('rbind',Results)
  R2_H = Results[,1]
  Corr.R2_H = Results[,2 ]
  r = Results[,3:(p+2)]
  if(Save.Corr){
    Lower.Dig.Corrs.All = Results[,-c(1:(p+2))]
    Outcome = list(R2_H = R2_H, Corr.R2_H = Corr.R2_H, Lower.Dig.Corrs.All = Lower.Dig.Corrs.All,surr.eval.r=r)    
  }else{
    Outcome = list(R2_H = R2_H, Corr.R2_H = Corr.R2_H,surr.eval.r=r)    
  }
  class(Outcome) <- "ICA.ContCont.MultS"
  return(Outcome)
}

# function to evaluate whether a correlation matrix is PD based on the eigenvalues
.is.PD = function(X, tol = 1e-08) {
  X[is.na(X)] = 0
  min.lambda = min(eigen(X, only.values = T, symmetric = T)$values)
  return(min.lambda > tol)
}
# function to estimate a single correlation r_[j,j+k] based on partial correlations
.r.random = function(j, k, R) {
  d = ncol(R)
  r1 = R[j, (j + 1):(j + k - 1)]
  r3 = R[(j + 1):(j + k - 1), j + k]
  R2.inv = R[(j + 1):(j + k - 1), (j + 1):(j + k - 1)]
  R2 = solve(R2.inv)
  #r.c = extraDistr::rnsbeta(1, parm.b, parm.b, -1, 1)
  r.c = extraDistr::rnsbeta(1, 1 + 0.5 * (d - 1 - k), 1 + 
                              0.5 * (d - 1 - k), -1, 1)
  D2 = (1 - tcrossprod(crossprod(r1, R2), r1)) * (1 - tcrossprod(crossprod(r3, 
                                                                           R2), r3))
  if (D2 < 0 & D2 > -1e-08) {
    D2 = 0
  }
  r = tcrossprod(crossprod(r1, R2), r3) + r.c * sqrt(D2)
  return(r)
}
# function to generate a random correlation matrix based on partial correlations
.Correlation.matrix.PC = function(R, Range = c(-1, 1),parm.a) {
  Rf = R
  Rf[is.na(Rf)] = 0
  d = ncol(R)
  j.ind = do.call("c", lapply(1:(d - 1), function(x) {
    x:1
  }))
  k.ind = do.call("c", lapply(1:(d - 1), function(x) {
    1:x
  }))
  for (i in 1:(length(j.ind))) {
    j = j.ind[i]
    k = k.ind[i]
    if (is.na(R[j, j + k])) {
      if (k == 1) {
        R[j, j + k] = R[j + k, j] = extraDistr::rnsbeta(1, parm.a, parm.a, -1, 1)
      }
      else {
        R[j, j + k] = R[j + k, j] = .r.random(j, k, R)
      }
      if (is.nan(R[j, j + k])) {
        stop("error")
      }
    }
  }
  return(R)
}
# function to compute the ICA
.MultivarICA.fun = function(R, Sigma, N) {
  d = nrow(R)
  sdMat = diag(sqrt(Sigma))
  rtn = sdMat %*% R %*% t(sdMat)
  var_diff <- function(cov_mat) {
    cov_val <- cov_mat[1, 1] + cov_mat[2, 2] - (2 * cov_mat[1, 
                                                            2])
    fit <- c(cov_val)
    fit
  }
  cov_2_diffs <- function(cov_mat) {
    cov_val <- (cov_mat[2, 2] - cov_mat[1, 2]) - (cov_mat[2, 
                                                          1] - cov_mat[1, 1])
    fit <- c(cov_val)
    fit
  }
  A <- matrix(var_diff(cov_mat = rtn[1:2, 1:2]), nrow = 1)
  B <- NULL
  Aantal <- (dim(R)[1] - 2)/2
  rtn_part <- rtn[c(3:dim(rtn)[1]), c(1, 2)]
  for (z in 1:Aantal) {
    cov_mat_here <- rtn_part[c((z * 2) - 1, z * 2), c(1:2)]
    B <- rbind(B, cov_2_diffs(cov_mat_here))
  }
  Dim <- dim(R)[1]
  Sub_mat_var <- rtn[c(3:Dim), c(3:Dim)]
  C <- matrix(NA, nrow = Aantal, ncol = Aantal)
  for (l in 1:Aantal) {
    for (k in 1:Aantal) {
      Sub_mat_var_hier <- Sub_mat_var[c((k * 2) - 1, 
                                        k * 2), c((l * 2) - 1, l * 2)]
      C[k, l] <- cov_2_diffs(cov_mat = Sub_mat_var_hier)
    }
  }
  Delta <- cbind(rbind(A, B), rbind(t(B), C))
  ICA <- (t(B) %*% solve(C) %*% B)/A
  Adj.ICA <- 1 - (1 - ICA) * ((N - 1)/(N - Aantal - 1))
  return(c(ICA, Adj.ICA))
}