R/fun_rc_scmprisk_Cox.R

In CopulaCenR: Copula-Based Regression Models for Multivariate Censored Data

bern_base1_log_lik_rc <- function(p, x1, indata1, b1, b1_d, m1, quantiles = NULL)
{
  
  phi1 <- p[1:(m1+1)]
  ep1 <- cumsum(exp(phi1))
  
  beta1 <- p[(m1+2):length(p)]
  
  t1<-indata1[,"obs_time"]
  
  # baseline
  Lambda1 <- b1%*%ep1
  lambda1 <- b1_d%*%ep1
  
  # survival and density
  u1 <- exp((-1)*Lambda1*exp(x1%*%beta1))
  f1 <- u1 * exp(x1%*%beta1) * lambda1
  
  
  term1 <- ifelse(indata1[,'status']==1, f1, 1)
  term1 <- log(abs(term1))
  
  term2 <- ifelse(indata1[,'status']==0, u1, 1)
  term2 <- log(abs(term2))
  
  logL <- (-1)*sum(term1 + term2)   
  return(logL)
  
}

bern_base2_log_lik_rc <- function(p, x2, indata2, b2, b2_d, m2, quantiles = NULL)
{
  
  phi2 <- p[1:(m2+1)]
  ep2 <- cumsum(exp(phi2))
  
  beta2 <- p[(m2+2):length(p)]
  
  t2<-indata2[,"obs_time"]
  
  # baseline
  Lambda2 <- b2%*%ep2
  lambda2 <- b2_d%*%ep2
  
  # survival and density
  u2 <- exp((-1)*Lambda2*exp(x2%*%beta2))
  f2 <- u2 * exp(x2%*%beta2) * lambda2
  
  
  term1 <- ifelse(indata2[,'status']==1, f2, 1)
  term1 <- log(abs(term1))
  
  term2 <- ifelse(indata2[,'status']==0, u2, 1)
  term2 <- log(abs(term2))
  
  logL <- (-1)*sum(term1 + term2)   
  return(logL)
  
}




bern_copula_pseudo_log_lik_rc <- function(p, fitted, x1, x2, indata1, indata2, b1, b2, b1_d, b2_d, m1, m2, quantiles = NULL, copula)
{
  
  eta <- exp(p[1])
  phi1 <- p[2:(m1+2)]
  ep1 <- cumsum(exp(phi1))
  beta1 <- p[(m1+3):length(p)]
  
  phi2 <- fitted[1:(m2+1)]
  ep2 <- cumsum(exp(phi2))
  beta2 <- fitted[(m2+2):length(fitted)]
  
  t1<-indata1[,"obs_time"]
  t2<-indata2[,"obs_time"]
  
  # survival and density
  u1 <- exp((-1)*(b1%*%ep1)*exp(x1%*%beta1))
  u2 <- exp((-1)*(b2%*%ep2)*exp(x2%*%beta2))
  f1 <- u1 * exp(x1%*%beta1) * (b1_d%*%ep1)
  f2 <- u2 * exp(x2%*%beta2) * (b2_d%*%ep2)
  
  # copula function
  if (copula == "Gumbel") {
    
    c_val<-((-log(u1))^(eta-1)*(-log(u2))^(eta-1)*gh_F(u1,u2,eta)*((-log(u1))^eta+(-log(u2))^eta)^(2/eta-2))/(u1*u2)-gh_F(u1,u2,eta)*((-log(u1))^eta+(-log(u2))^eta)^(1/eta-2)*(1/(u1*u2))*(1/eta-1)*eta*(-log(u1))^(eta-1)*(-log(u2))^(eta-1)
    c_u1_val<-gh_F(u1,u2,eta)*((-log(u1))^eta+(-log(u2))^eta)^(1/eta-1)*(-log(u1))^(eta-1)/u1
    c_u2_val<-gh_F(u1,u2,eta)*((-log(u1))^eta+(-log(u2))^eta)^(1/eta-1)*(-log(u2))^(eta-1)/u2
    C_val<-gh_F(u1,u2,eta)
    
  }
  
  
  if (copula == "Clayton") {
    
    # Clayton Copula function for joint distribution probability
    c_val<-clt_f(u1,u2,eta)
    c_u1_val<-u1^(-eta-1)*(u1^(-eta)+u2^(-eta)-1)^(-1/eta-1) 
    c_u2_val<-u2^(-eta-1)*(u1^(-eta)+u2^(-eta)-1)^(-1/eta-1)
    C_val<-(u1^(-eta)+u2^(-eta)-1)^(-1/eta) # C(u,v)
    
  }
  
  
  C_val <- ifelse((indata1[,"status"] == 0) & (indata2[,"status"] == 0), C_val, 1)
  term1 <- log(abs(C_val))
  
  term2 <- c_u1_val * f1
  term2 <- ifelse((indata1[,"status"] == 1) & (indata2[,"status"] == 0), term2, 1)
  term2 <- log(abs(term2))
  
  term3 <- c_u2_val * f2
  term3 <- ifelse((indata1[,"status"] == 0) & (indata2[,"status"] == 1), term3, 1)
  term3 <- log(abs(term3))
  
  term4 <- c_val * f1 * f2
  term4 <- ifelse((indata1[,"status"] == 1) & (indata2[,"status"] == 1), term4, 1)
  term4[term4 < 0] <- 1
  term4 <- log(abs(term4))
  # cat(paste0('eta: ',eta), '\n')
  logL <- (-1)*sum( term1 + term2 + term3 + term4 )
  
  return(logL)   
}



bern_copula_log_lik_rc <- function(p, x1, x2, indata1, indata2, b1, b2, b1_d, b2_d, m1, m2, beta1_b_1, quantiles = NULL, copula)
{
  
  eta <- exp(p[1])
  phi1 <- p[2:(m1+2)]
  ep1 <- cumsum(exp(phi1))
  beta1 <- p[(m1+3):(m1+2+length(beta1_b_1))]
  phi2 <- p[(m1+3+length(beta1_b_1)):(m1+3+m2+length(beta1_b_1))]
  ep2 <- cumsum(exp(phi2))
  beta2 <- p[(m1+4+m2+length(beta1_b_1)):length(p)]
  
  t1<-indata1[,"obs_time"]
  t2<-indata2[,"obs_time"]
  
  # survival and density
  u1 <- exp((-1)*(b1%*%ep1)*exp(x1%*%beta1))
  u2 <- exp((-1)*(b2%*%ep2)*exp(x2%*%beta2))
  f1 <- u1 * exp(x1%*%beta1) * (b1_d%*%ep1)
  f2 <- u2 * exp(x2%*%beta2) * (b2_d%*%ep2)
  
  # copula function
  if (copula == "Gumbel") {
    
    c_val<-((-log(u1))^(eta-1)*(-log(u2))^(eta-1)*gh_F(u1,u2,eta)*((-log(u1))^eta+(-log(u2))^eta)^(2/eta-2))/(u1*u2)-gh_F(u1,u2,eta)*((-log(u1))^eta+(-log(u2))^eta)^(1/eta-2)*(1/(u1*u2))*(1/eta-1)*eta*(-log(u1))^(eta-1)*(-log(u2))^(eta-1)
    c_u1_val<-gh_F(u1,u2,eta)*((-log(u1))^eta+(-log(u2))^eta)^(1/eta-1)*(-log(u1))^(eta-1)/u1
    c_u2_val<-gh_F(u1,u2,eta)*((-log(u1))^eta+(-log(u2))^eta)^(1/eta-1)*(-log(u2))^(eta-1)/u2
    C_val<-gh_F(u1,u2,eta)
    
  }
  
  
  if (copula == "Clayton") {
    
    # Clayton Copula function for joint distribution probability
    c_val<-clt_f(u1,u2,eta)
    c_u1_val<-u1^(-eta-1)*(u1^(-eta)+u2^(-eta)-1)^(-1/eta-1) 
    c_u2_val<-u2^(-eta-1)*(u1^(-eta)+u2^(-eta)-1)^(-1/eta-1)
    C_val<-(u1^(-eta)+u2^(-eta)-1)^(-1/eta) # C(u,v)
    
  }
  
  
  C_val <- ifelse((indata1[,"status"] == 0) & (indata2[,"status"] == 0), C_val, 1)
  term1 <- log(abs(C_val))
  
  term2 <- c_u1_val * f1
  term2 <- ifelse((indata1[,"status"] == 1) & (indata2[,"status"] == 0), term2, 1)
  term2 <- log(abs(term2))
  
  term3 <- c_u2_val * f2
  term3 <- ifelse((indata1[,"status"] == 0) & (indata2[,"status"] == 1), term3, 1)
  term3 <- log(abs(term3))
  
  term4 <- c_val * f1 * f2
  term4 <- ifelse((indata1[,"status"] == 1) & (indata2[,"status"] == 1), term4, 1)
  term4[term4 < 0] <- 1
  term4 <- log(abs(term4))
  
  logL<-(-1)*sum( term1 + term2 + term3 + term4 )
  return(logL)   
  
}