Nothing
#' @title Fit Cox Proportional Hazard Regression Model under dependent right censoring via MPL and Archimedean Copulas
#'
#' @description Simultaneously estimate the regression coefficients and the baseline hazard function of proportional
#' hazard Cox models under dependent right censoring using maximum penalised likelihood (MPL) and Archimedean Copulas
#' @usage coxph_mpl_dc(surv, cova, control,...)
#'
#' @param surv the outcome of survival data, with the first column X (observed time), second column
#' del (event indicator) and third column eta (dependent right censoring indicator).
#' @param cova the covariate matrix, with dimension of n rows and p columns, where 'n' is the sample size and 'p' is the number of covariates
#' Default is \code{cova=matrix(0,n,1)}, which the covariates are not considered.
#' @param control object of class \code{coxph_mpl_dc.control} specifying control options like basis choice,
#' refer to \code{coxph_mpl_dc.control} to see the defaults.
#' @param ... other arguments. In \code{coxph_mpl_dc}, these elements, will be passed to \code{\link{coxph_mpl_dc.control}}.
#'
#'
#' @details coxph_mpl_dc allows to simultaneously estimate the regression coefficients and baseline hazard function of Cox proportional hazard models,
#' with dependent and independent right censored data, by maximizing a penalized likelihood,
#' in which a penalty function is used to smooth the baseline hazard estimates.
#' Note that the dependence between event and censoring times is modelled by an Archimedean copula.
#'
#' Optimization is achieved using an iterative algorithm, which combines Newton’s method and the multiplicative iterative algorithm proposed by Ma (2010),
#' and respects the non-negativity constraints on the baseline hazard estimate (refer to Ma et al (2014) and Xu et al (2018)).
#'
#' The centered covariate matrix \eqn{Z} is used in the optimization process to get a better shaped (penalized) log-likelihood.
#' Baseline hazard parameter estimates and covariance matrix are then respectively corrected using a correction factor and the delta method.
#'
#' The estimates of zero are possible for baseline hazard parameters (e.g., when number of knots is relatively large to sample size) and
#' will correspond to active constraints as defined by Moore and Sadler (2008). Inference, as described by Ma et al (2014) or Xu et al (2018),
#' is then corrected accordingly (refer to Moore and Sadler (2008)) by adequately cutting the corresponding covariance matrix.
#'
#' There are currently three ways to perform inference on model parameters: Let \eqn{H} denote the negative of Hessian matrix of the non-penalized likelihood,
#' \eqn{Q} denote the product of the first order derivative of the penalized likelihood by its transpose,
#' and \eqn{M_2} denote the negative of the second order derivative of the penalized likelihood.
#' Then the three estimated covariance matrices for the MPL estimates are \eqn{M_{2}^{-1}}, \eqn{M_{2}^{-1}HM_{2}^{-1}} and \eqn{M_{2}^{-1}QM_{2}^{-1}}.
#'
#' Simulation studies the coverage levels of confidence intervals for the regression parameters seem to indicate
#' \eqn{M_{2}^{-1}HM_{2}^{-1}} performs better when using the piecewise constant basis,
#' and that \eqn{M_{2}^{-1}QM_{2}^{-1}} performs better when using other bases.
#'
#'
#' @return \item{mpl_theta}{MPL estimates of the regression coefficient for the basis functions of the baseline hazard of T, i.e. theta}
#' \item{mpl_gamma}{MPL estimates of the regression coefficient for the basis functions of the baseline hazard of C, i.e. gamma}
#' \item{mpl_h0t}{MPL estimates of the baseline hazard for T, i.e. h_{0T}(x_i)}
#' \item{mpl_h0c}{MPL estimates of the baseline hazard for C, i.e. h_{0C}(x_i)}
#' \item{mpl_H0t}{MPL estimates of the baseline cumulative hazard for T, i.e. H_{0T}(x_i)}
#' \item{mpl_H0c}{MPL estimates of the baseline cumulative hazard for C, i.e. H_{0c}(x_i)}
#' \item{mpl_S0t}{MPL estimates of the baseline survival for T, i.e. S_{0T}(x_i)}
#' \item{mpl_S0c}{MPL estimates of the baseline survival for C, i.e. S_{0C}(x_i)}
#' \item{mpl_beta}{MPL estimates of beta}
#' \item{mpl_phi}{MPL estimates of phi}
#' \item{penloglik}{the penalized log-likelihood function given the MPL estimates}
#' \item{mpl_Ubeta}{the first derivative of penalized log-likelihood function with respect to beta given the MPL estimates}
#' \item{mpl_Uphi}{the first derivative of penalized log-likelihood unction with respect to phi given the MPL estimates}
#' \item{mpl_Utheta}{the first derivative of penalized log-likelihood function with respect to theta given the MPL estimates}
#' \item{mpl_Ugamma}{the first derivative of penalized log-likelihood function with respect to gamma given the MPL estimates}
#' \item{iteration}{a vector of length 3 indicating the number of iterations used to estimate the smoothing parameter
#' (first value, equal to 1 when the user specified a chosen value), the beta, phi, theta and gamma parameters during the entire process (second value),
#' and beta, phi, theta and gamma parameters during the last smoothing parameter iteration (third value)}
#' \item{mpl_cvl}{the cross validation value given the MPL estimates}
#' \item{mpl_aic}{the AIC value given the MPL estimates}
#' \item{mpl_beta_sd}{the asymptotic standard deviation of the MPL estimated beta}
#' \item{mpl_phi_sd}{the asymptotic standard deviation of the MPL estimated phi}
#' \item{mpl_h0t_sd}{the asymptotic standard deviation of the MPL estimates for the baseline hazard coefficient of T, i.e. theta}
#' \item{mpl_h0c_sd}{the asymptotic standard deviation of the MPL estimates for the baseline hazard coefficient of C, i.e. gamma}
#'
#' \item{mpl_ht0_sd}{the asymptotic standard deviation of the MPL estimates for the baseline hazard of T}
#' \item{mpl_hc0_sd}{the asymptotic standard deviation of the MPL estimates for the baseline hazard of C}
#' \item{mpl_Ht0_sd}{the asymptotic standard deviation of the MPL estimates for the cumulative baseline hazard of T}
#' \item{mpl_Hc0_sd}{the asymptotic standard deviation of the MPL estimates for the cumulative baseline hazard of C}
#' \item{mpl_St0_sd}{the asymptotic standard deviation of the MPL estimates for the baseline survival of T}
#' \item{mpl_Sc0_sd}{the asymptotic standard deviation of the MPL estimates for the baseline survival of C}
#'
#' \item{mpl_est_cov}{the asymptotic covariance matrix of the MPL estimates}
#' \item{mpl_beta_phi_zp}{the MPL estimates for regression coefficient with their corresponding standard deviations, z scores and p-values}
#' \item{binwv}{the width of each discretized bin of the observed times when piecewise constant approximation applied to the baseline hazards}
#' \item{ID}{the bin ID for each individual of the sample when piecewise constant approximation applied to the baseline hazards}
#' \item{binedg}{the edge for each discretized bin among the observed time vector X, which are the internal knots and boundaries}
#' \item{psix}{basis function matrix psi(x_i) with dimension of n by m for baseline hazard, where m=number of internal knots+ordSp}
#' \item{Psix}{basis function matrix Psi(x_i) with dimension of n by m for baseline cumulative hazard}
#' Inputs defined in coxph_mpl_dc.control
#'
#' @author Jing Xu, Jun Ma, Thomas Fung
#'
#' @references Ma, J. (2010). \emph{"Positively constrained multiplicative iterative algorithm for maximum penalised likelihood tomographic reconstruction"}.
#' IEEE Transactions On Signal Processing 57, 181-192.
#'
#' Ma, J. and Heritier, S. and Lo, S. (2014). \emph{"On the Maximum Penalised Likelihood Approach forProportional Hazard Models with Right Censored Survival Data"}.
#' Computational Statistics and Data Analysis 74, 142-156.
#'
#' Xu J, Ma J, Connors MH, Brodaty H. (2018). \emph{"Proportional hazard model estimation under dependent censoring using copulas and penalized likelihood"}.
#' Statistics in Medicine 37, 2238–2251.
#'
#'
#' @seealso \code{\link{plot.coxph_mpl_dc}}, \code{\link{coxph_mpl_dc.control}}, \code{\link{coef.coxph_mpl_dc}}
#'
#' @importFrom copula indepCopula claytonCopula gumbelCopula frankCopula iTau rCopula
#' @importFrom survival coxph Surv
#' @importFrom splines2 iSpline mSpline
#' @importFrom matrixcalc is.positive.definite
#'
#' @examples
#' \donttest{
#' ##-- Copula types
#' copula3 <- 'frank'
#'
#' ##-- Marginal distribution for T, C, and A
#' a <- 2
#' lambda <- 2
#' cons7 <- 0.2
#' cons9 <- 10
#' tau <- 0.8
#' betas <- c(-0.5, 0.1)
#' phis <- c(0.3, 0.2)
#' distr.ev <- 'weibull'
#' distr.ce <- 'exponential'
#'
#' ##-- Sample size
#' n <- 200
#'
#' ##-- One sample Monte Carlo dataset
#' cova <- cbind(rbinom(n, 1, 0.5), runif(n, min=-10, max=10))
#' surv <- surv_data_dc(n, a, cova, lambda, betas, phis, cons7, cons9, tau, copula3,
#' distr.ev, distr.ce)
#' n <- nrow(cova)
#' p <- ncol(cova)
#' ##-- event and dependent censoring proportions
#' colSums(surv)[c(2,3)]/n
#' X <- surv[,1] # Observed time
#' del<-surv[,2] # failure status
#' eta<-surv[,3] # dependent censoring status
#'
#' ##-- control inputs for the coxph_mpl_dc function
#' control <- coxph_mpl_dc.control(ordSp = 4,
#' binCount = 100,
#' tau = 0.8, copula = copula3,
#' pent = 'penalty_mspl', smpart = 'REML',
#' penc = 'penalty_mspl', smparc = 'REML',
#' cat.smpar = 'No' )
#'
#' ##-- Fitting cox ph hazard model for T using MPL and an correct copula
#' #with REML smoothing parameters
#' coxMPLests5 <- coxph_mpl_dc(surv, cova, control, )
#' mpl_beta_phi_zp5 <- coxMPLests5$mpl_beta_phi_zp
#' mpl_h0t5 <- coxMPLests5$mpl_h0t
#' mpl_h0Ti5 <- approx( X, mpl_h0t5, xout = seq(0, 5.4, 0.01),
#' method="constant", rule = 2, ties = mean)$y
#'
#' ##-- Real marginal baseline hazard for T
#' ht0b <- a * (seq(0, 5.4, 0.01) ^ (a - 1)) / (lambda ^ a)
#'
#'
#' ##-- Fitting cox ph hazard model for T using MPL and an correct copula
#' #with zero smoothing parameters
#' coxMPLests3 <- coxph_mpl_dc(surv, cova,
#' ordSp=4, binCount=100,
#' tau=0.8, copula=copula3,
#' pent='penalty_mspl', smpart=0, penc='penalty_mspl', smparc=0,
#' cat.smpar = 'No')
#' mpl_beta_phi_zp3 <- coxMPLests3$mpl_beta_phi_zp
#' mpl_h0t3 <- coxMPLests3$mpl_h0t
#' mpl_h0Ti3 <- approx( X, mpl_h0t3, xout = seq(0, 5.4, 0.01),
#' method="constant", rule = 2, ties = mean)$y
#'
#' ##-- Plot the true and estimated baseline hazards for T
#' t_up <- 3.5
#' y_uplim <- 2
#' Ti<-seq(0, 5.4, 0.01)[seq(0, 5.4, 0.01)<=t_up]
#' h0Ti<-ht0b[seq(0, 5.4, 0.01)<=t_up]
#' h0Ti5<-mpl_h0Ti5[seq(0, 5.4, 0.01)<=t_up]
#' h0Ti3<-mpl_h0Ti3[seq(0, 5.4, 0.01)<=t_up]
#'
#' plot( x = Ti, y = h0Ti5,
#' type="l", col="grey", lty=4, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim),
#' xlab='Time', ylab='Hazard')
#' lines(x = Ti, y = h0Ti,
#' col="green",
#' lty=1, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim)
#' )
#' lines(x = Ti, y = h0Ti3,
#' col="blue",
#' lty=4, lwd=3, cex.axis=1.6, cex.lab=1.6, ylim=c(0, y_uplim)
#' )
#'
#'}
#'
#'
#'
#' @export
coxph_mpl_dc<-function(surv, cova, control,...)
{
if(missing(control)){
control<-coxph_mpl_dc.control(...)}
ordSp<-control$ordSp
binCount<-control$binCount
tie<-control$tie
tau<-control$tau
copula<-control$copula
pent<-control$pent
smpart<-control$smpart
penc<-control$penc
smparc<-control$smparc
maxit2<-control$maxit2
maxit<-control$maxit
mid<-control$mid
asy<-control$asy
ac<-control$ac
cv<-control$cv
ac.theta<-control$ac.theta
ac.gamma<-control$ac.gamma
ac.Utheta<-control$ac.Utheta
ac.Ugamma<-control$ac.Ugamma
min.theta<-control$min.theta
min.gamma<-control$min.gamma
min.ht<-control$min.ht
min.hc<-control$min.hc
min.St<-control$min.St
min.Sc<-control$min.Sc
min.C<-control$min.C
min.dC<-control$min.dC
eps<-control$eps
tol.thga<-control$tol.thga
tol.bph<-control$tol.bph
tol.smpar<-control$tol.smpar
cat.smpar<-control$cat.smpar
if(missing(cova)){cova<-matrix(0, dim(surv)[1],1)}
if(missing(ordSp)){ordSp<-1}
if(missing(binCount)){binCount<-1}
if(missing(tie)){tie<-'No'}
if(missing(tau)){tau<-0}
if(missing(copula)){copula<-'independent'}
if(missing(pent)){pent<-'mat1'}
if(missing(penc)){penc<-'mat1'}
if(missing(smpart)){smpart<-0}
if(missing(smparc)){smparc<-0}
if(missing(maxit2)){maxit2<-50}
if(missing(maxit)){maxit<-5000}
if(missing(mid)){mid<-1}
if(missing(asy)){asy<-1}
if(missing(ac)){ac<-0}
if(missing(cv)){cv<-0}
if(missing(ac.theta)){ac.theta<-1e-5}
if(missing(ac.gamma)){ac.gamma<-1e-5}
if(missing(ac.Utheta)){ac.Utheta<--1e-2}
if(missing(ac.Ugamma)){ac.Ugamma<--1e-2}
if(missing(min.theta)){min.theta<-1e-7}
if(missing(min.gamma)){min.gamma<-1e-7}
if(missing(min.ht)){min.ht<-1e-7}
if(missing(min.hc)){min.hc<-1e-7}
if(missing(min.St)){min.St<-1e-7}
if(missing(min.Sc)){min.Sc<-1e-7}
if(missing(min.C)){min.C<-1e-7}
if(missing(min.dC)){min.dC<-1e-7}
if(missing(eps)){eps<-1e-5}
if(missing(tol.thga)){tol.thga<-1e-5}
if(missing(tol.bph)){tol.bph<-1e-5}
if(missing(cat.smpar)){cat.smpar<-'Yes'}
if(missing(tol.smpar)){tol.smpar<-1e-2}
smpart_init<-smpart
smparc_init<-smparc
if(smpart=='REML'){smpart<-0}
if(smparc=='REML'){smparc<-0}
{
if(copula=='independent'){alpha<-0}
else if(copula=='clayton'){alpha<-iTau(claytonCopula(100), tau)} #calculate the alpha value of Clayton copula give the tau value
else if(copula=='gumbel'){alpha<-iTau(gumbelCopula(100), tau)}
else if (copula=='frank'){alpha<- -iTau(frankCopula(100), tau)}
else (stop("Error! Define a proper copula function"))
}
del<-surv[,2]
eta<-surv[,3]
X<-surv[,1]
discrsurvdat<-discrBinNA(surv, binCount, tie)
groupsurvdat <- discrsurvdat$discretize
binwv<-groupsurvdat[,1] #the bin widths
ID<-discrsurvdat$ID #the bin ID of each subject
binedg<-discrsurvdat$binedg
n<-nrow(cova) #sample size
p<-ncol(cova)#number of covariates for each subject
numIntKnt<-length(binedg)-2
IntKnt<-binedg[2:(numIntKnt+1)]
bryKnt<-c(0, max(binedg))
numSp<-numIntKnt+ordSp
m<-numSp
tcova<-cova #regression covariates
cova <- cova - t(matrix(apply(cova,2,mean),p,n))
if(ordSp==1)
{
psix<-psi(ID);
Psix<-Psi(ID, binwv)
{
if(pent=='mat1') Rt<-mat1(psix, X, eps)
else if (pent=='mat2') Rt<-mat2(psix, X, eps)
else (stop("Error! Define a proper Penalty function for T"))
}
{
if(penc=='mat1') Rc<-mat1(psix, X, eps)
else if (penc=='mat2') Rc<-mat2(psix, X, eps)
else (stop("Error! Define a proper Penalty function for C"))
}
}
else
{
psix<-mSpline(X, knots = IntKnt, degree = ordSp-1, intercept = TRUE, derivs=0, Boundary.knots=bryKnt )
Psix<-iSpline(X, knots = IntKnt, degree = ordSp-1, intercept = TRUE, derivs=0, Boundary.knots=bryKnt )
{
if(pent=='penalty_mspl') Rt<-penalty_mspl(numSp, ordSp, IntKnt, bryKnt)
else (stop("Error! Define a proper Penalty function for T"))
}
{
if(penc=='penalty_mspl') Rc<-penalty_mspl(numSp, ordSp, IntKnt, bryKnt)
else (stop("Error! Define a proper Penalty function for C"))
}
}
if( identical( cova, matrix( 0, n, 1 ) ) )
{
beta0 <- 0
phi0 <- 0
}
else
{
coxpht<-coxph( Surv(surv[,1], surv[,2])~ tcova )
beta0<-coxpht$coef
beta0 <- rep( 0, dim(tcova)[2] )
coxphc<-coxph( Surv(surv[,1], surv[,3]) ~ tcova )
phi0<-coxphc$coef
phi0 <- rep( 0, dim(tcova)[2] )
}
beta_old<-beta0 # initial estimates of the regression coefficient of the proportional hazard for failure time
phi_old<-phi0 # initial estimates of the regression coefficient of the proportional hazard for dependent censoring time
theta0<-theta_initial(del, psix, Psix, beta0, tcova) #initial estimates of theta's assumed independent censoring
gamma0<-gamma_initial(eta, psix, Psix, phi0, tcova) #initial estimates of gamma's assumed independent censoring
theta_old<-theta0
gamma_old<-gamma0
theta_old[theta_old<min.theta]<-min.theta
gamma_old[gamma_old<min.gamma]<-min.gamma
h0t0<-baseHaz(theta_old, psix);
h0c0<-baseHaz(gamma_old, psix)
h0t_old<-h0t0 #intital estimates of the baseline hazards for failure time
h0c_old<-h0c0 #intital estimates of the baseline hazards for dependent censoring time
h0t_old[h0t_old<min.ht]<-min.ht
h0c_old[h0c_old<min.hc]<-min.hc
H0t_old<-baseCumm(theta_old, Psix)
H0c_old<-baseCumm(gamma_old, Psix)
S0t_old<-exp(-H0t_old)
S0c_old<-exp(-H0c_old)
regt_old<-cova%*%beta_old
regc_old<-cova%*%phi_old
eregt_old<-exp(regt_old)
eregc_old<-exp(regc_old)
Hcoxt_old<-H0t_old*eregt_old
Hcoxc_old<-H0c_old*eregc_old
Scoxt_old<-exp(-Hcoxt_old) #initial estimates of the survival function of failure based on Cox proportional hazard
Scoxc_old<-exp(-Hcoxc_old) #initial estimates of the survival function of dependent censoring based on Cox proportional hazard
Scoxt_old[Scoxt_old<min.St]<-min.St
Scoxc_old[Scoxc_old<min.Sc]<-min.Sc
{
if(copula=='independent')
{
C_old<-IndependentCopula(Scoxt_old, Scoxc_old) # independent copula function value
dC_old<-dC_ind(Scoxt_old, Scoxc_old)# the first and second derivate of the independent copula function
dC3_old<-dC3_ind(Scoxt_old, Scoxc_old)#the third derivative of the independent copula functiuon
}
else if(copula=='clayton')
{
C_old<-ClaytonCopula(Scoxt_old, Scoxc_old, alpha)
dC_old<-dC_clay(Scoxt_old, Scoxc_old, alpha)
dC3_old<-dC3_clay(Scoxt_old, Scoxc_old, alpha)
}
else if(copula=='gumbel')
{
C_old<-GumbelCopula(Scoxt_old, Scoxc_old, alpha)
dC_old<-dC_gumbel(Scoxt_old, Scoxc_old, alpha)
dC3_old<-dC3_gumbel(Scoxt_old, Scoxc_old, alpha)
}
else if (copula=='frank')
{
C_old<-FrankCopula(Scoxt_old, Scoxc_old, alpha)
dC_old<-dC_frank(Scoxt_old, Scoxc_old, alpha)
dC3_old<-dC3_frank(Scoxt_old, Scoxc_old, alpha)
}
else (stop("Error! Define a proper copula function"))
}
C_old[ C_old == Inf ] <- 1
C_old[ C_old == -Inf ] <- 0
dC_old[,1:3][dC_old[,1:3]<min.dC]<-min.dC
C_old[C_old<min.C]<-min.C
Rt_star <- rbind( matrix(0, 2*p, 2*p + 2*m), cbind( matrix(0, m, 2*p), Rt, matrix(0,m,m) ), matrix(0, m, 2*p+2*m) )
Rc_star <- rbind( matrix(0, 2*p, 2*p + 2*m), matrix(0,m,2*p+2*m), cbind( matrix(0, m, 2*p+m), Rc ) )
full.iter <- 0
dftold <- m
dfcold <- m
for(k2 in 1: maxit2)
{
pllik<-0 #penalized log likelihood function
pllik_old<-penlogreg_dep(del, eta, theta_old, gamma_old, h0t_old, h0c_old, regt_old, regc_old, Hcoxt_old, Hcoxc_old, C_old, dC_old, Rt, Rc, smpart, smparc)
pllik[1]<-pllik_old
for(k in 1:maxit)
{
coef_old<-c(beta_old, phi_old)
Ubeta<-U_beta(p, cova, Hcoxt_old, Scoxt_old, del, eta, dC_old, C_old) #score function respect to beta
Uphi<-U_phi(p, cova, Hcoxc_old, Scoxc_old, del, eta, dC_old, C_old) #score function respect to phi
Ucoef<-c(Ubeta, Uphi)
Hbb<-Hess_beta(p, C_old, dC_old, dC3_old, cova, Hcoxt_old, Scoxt_old, del, eta) #hessian matrix of beta
Hpp<-Hess_phi(p, C_old, dC_old, dC3_old, cova, Hcoxc_old, Scoxc_old, del, eta) #hessian matrix of phi
Hpb<-Hess_phi_beta(p, C_old, dC_old, dC3_old, cova, Hcoxt_old, Hcoxc_old, Scoxt_old, Scoxc_old, del, eta)
Hbp<-t(Hpb)
Hcoef<- -rbind(cbind(Hbb, Hbp), cbind(Hpb, Hpp))
if(is.positive.definite( round(Hcoef, digits = 4) )==TRUE){
Hcoef<-Hcoef
}else{
Ubet_i<-dplcox_dbeta_i( p, cova, Hcoxt_old, Scoxt_old, del, eta, dC_old, C_old )
Uph_i<-dplcox_dphi_i( p, cova, Hcoxc_old, Scoxc_old, del, eta, dC_old, C_old )
U_i<-cbind(Ubet_i, Uph_i)
Hcoef<-t(U_i)%*%U_i
}
diag(Hcoef)<-diag(Hcoef)+eps
inccoef<-solve(Hcoef)%*%Ucoef
coef_new<-coef_old+inccoef #updating beta using Newton algorithm
beta_new<-coef_new[1:p]
phi_new<-coef_new[(p+1):(2*p)]
regt_new<-cova%*%beta_new
regc_new<-cova%*%phi_new
eregt_new<-exp(regt_new)
eregc_new<-exp(regc_new)
#eregt_new[eregt_new==Inf]=max( eregt_new[eregt_new<Inf] )
#eregc_new[eregc_new==Inf]=max( eregc_new[eregc_new<Inf] )
Hcoxt_new1<-H0t_old*eregt_new
Hcoxc_new1<-H0c_old*eregc_new
Scoxt_new1<-exp(-Hcoxt_new1)
Scoxc_new1<-exp(-Hcoxc_new1)
Scoxt_new1[Scoxt_new1<min.St]<-min.St
Scoxc_new1[Scoxc_new1<min.Sc]<-min.Sc
{
if(copula=='independent')
{
C_new1<-IndependentCopula(Scoxt_new1, Scoxc_new1)
dC_new1<-dC_ind(Scoxt_new1, Scoxc_new1)
}
else if(copula=='clayton')
{
C_new1<-ClaytonCopula(Scoxt_new1, Scoxc_new1, alpha)
dC_new1<-dC_clay(Scoxt_new1, Scoxc_new1, alpha)
}
else if(copula=='gumbel')
{
C_new1<-GumbelCopula(Scoxt_new1, Scoxc_new1, alpha)
dC_new1<-dC_gumbel(Scoxt_new1, Scoxc_new1, alpha)
}
else if (copula=='frank')
{
C_new1<-FrankCopula(Scoxt_new1, Scoxc_new1, alpha)
dC_new1<-dC_frank(Scoxt_new1, Scoxc_new1, alpha)
}
else (stop("Error! Define a proper copula function"))
}
C_new1[ C_new1 == Inf ] <- 1
C_new1[ C_new1 == -Inf ] <- 0
dC_new1[,1:3][dC_new1[,1:3]<min.dC]<-min.dC
C_new1[C_new1<min.C]<-min.C
pllik_new1<-penlogreg_dep(del, eta, theta_old, gamma_old, h0t_old, h0c_old, regt_new, regc_new, Hcoxt_new1, Hcoxc_new1, C_new1, dC_new1, Rt, Rc, smpart, smparc)
Ubeta_new1<-U_beta(p, cova, Hcoxt_new1, Scoxt_new1, del, eta, dC_new1, C_new1) #score function respect to beta
Uphi_new1<-U_phi(p, cova, Hcoxc_new1, Scoxc_new1, del, eta, dC_new1, C_new1) #score function respect to phi
ome<-0.6
while( pllik_new1<pllik_old | min( c( Ubeta_new1/Ubeta, Uphi_new1/Uphi ) )<0 )
{
coef_new<-coef_old+ome*inccoef
beta_new<-coef_new[1:p]
phi_new<-coef_new[(p+1):(2*p)]
regt_new<-cova%*%beta_new
regc_new<-cova%*%phi_new
eregt_new<-exp(regt_new)
eregc_new<-exp(regc_new)
Hcoxt_new1<-H0t_old*eregt_new
Hcoxc_new1<-H0c_old*eregc_new
Scoxt_new1<-exp(-Hcoxt_new1)
Scoxc_new1<-exp(-Hcoxc_new1)
Scoxt_new1[Scoxt_new1<min.St]<-min.St
Scoxc_new1[Scoxc_new1<min.Sc]<-min.Sc
{
if(copula=='independent')
{
C_new1<-IndependentCopula(Scoxt_new1, Scoxc_new1)
dC_new1<-dC_ind(Scoxt_new1, Scoxc_new1)
}
else if(copula=='clayton')
{
C_new1<-ClaytonCopula(Scoxt_new1, Scoxc_new1, alpha)
dC_new1<-dC_clay(Scoxt_new1, Scoxc_new1, alpha)
}
else if(copula=='gumbel')
{
C_new1<-GumbelCopula(Scoxt_new1, Scoxc_new1, alpha)
dC_new1<-dC_gumbel(Scoxt_new1, Scoxc_new1, alpha)
}
else if (copula=='frank')
{
C_new1<-FrankCopula(Scoxt_new1, Scoxc_new1, alpha)
dC_new1<-dC_frank(Scoxt_new1, Scoxc_new1, alpha)
}
else (stop("Error! Define a proper copula function"))
}
C_new1[ C_new1 == Inf ] <- 1
C_new1[ C_new1 == -Inf ] <- 0
dC_new1[,1:3][dC_new1[,1:3]<min.dC]<-min.dC
C_new1[C_new1<min.C]<-min.C
pllik_new1<-penlogreg_dep(del, eta, theta_old, gamma_old, h0t_old, h0c_old, regt_new, regc_new, Hcoxt_new1, Hcoxc_new1,
C_new1, dC_new1, Rt, Rc, smpart, smparc)
Ubeta_new1<-U_beta(p, cova, Hcoxt_new1, Scoxt_new1, del, eta, dC_new1, C_new1) #score function respect to beta
Uphi_new1<-U_phi(p, cova, Hcoxc_new1, Scoxc_new1, del, eta, dC_new1, C_new1) #score function respect to phi
if(ome >= 1e-2) ome=ome*0.6
else if(ome < 1e-2 & ome >= 1e-5) ome=ome*5e-2
else if(ome < 1e-5 & ome >= 1e-20) ome=ome*1e-5
else
{
print("Iterations do not move coeficients")
break
}
}
denomt<-denomtreg(Rt, psix, Psix, theta_old, dC_new1, C_new1, Scoxt_new1, eregt_new, del, eta, smpart)
st<-theta_old/(denomt+1e-5) #step of function to update theta
Uth<-U_theta(psix, Psix, theta_old, h0t_old, eregt_new, Scoxt_new1, del, eta, dC_new1, C_new1, smpart, Rt) #score function respect to theta
theta_new<-theta_old+st*Uth#updating theta using MI algorithm
denomc<-denomcreg(Rc, psix, Psix, gamma_old, dC_new1, C_new1, Scoxc_new1, eregc_new, del, eta, smparc)
sc<-gamma_old/(denomc++1e-5)
Uga<-U_gamma(psix, Psix, gamma_old, h0c_old, eregc_new, Scoxc_new1, del, eta, dC_new1, C_new1, smparc, Rc)
gamma_new<-gamma_old+sc*Uga
theta_new[theta_new<min.theta]<-min.theta
gamma_new[gamma_new<min.gamma]<-min.gamma
h0t_new<-baseHaz(theta_new, psix)
h0c_new<-baseHaz(gamma_new, psix)
h0t_new[h0t_new<min.ht]<-min.ht
h0c_new[h0c_new<min.hc]<-min.hc
H0t_new<-baseCumm(theta_new,Psix)
H0c_new<-baseCumm(gamma_new,Psix)
S0t_new<-exp(-H0t_new)
S0c_new<-exp(-H0c_new)
Hcoxt_new2<-H0t_new*eregt_new
Hcoxc_new2<-H0c_new*eregc_new
Scoxt_new2<-exp(-Hcoxt_new2)
Scoxc_new2<-exp(-Hcoxc_new2)
Scoxt_new2[Scoxt_new2<min.St]<-min.St
Scoxc_new2[Scoxc_new2<min.Sc]<-min.Sc
{
if(copula=='independent')
{
C_new2<-IndependentCopula(Scoxt_new2, Scoxc_new2)
dC_new2<-dC_ind(Scoxt_new2, Scoxc_new2)
dC3_new2<-dC3_ind(Scoxt_new2, Scoxc_new2)
}
else if(copula=='clayton')
{
C_new2<-ClaytonCopula(Scoxt_new2, Scoxc_new2, alpha)
dC_new2<-dC_clay(Scoxt_new2, Scoxc_new2, alpha)
dC3_new2<-dC3_clay(Scoxt_new2, Scoxc_new2, alpha)
}
else if(copula=='gumbel')
{
C_new2<-GumbelCopula(Scoxt_new2, Scoxc_new2, alpha)
dC_new2<-dC_gumbel(Scoxt_new2, Scoxc_new2, alpha)
dC3_new2<-dC3_gumbel(Scoxt_new2, Scoxc_new2, alpha)
}
else if (copula=='frank')
{
C_new2<-FrankCopula(Scoxt_new2, Scoxc_new2, alpha)
dC_new2<-dC_frank(Scoxt_new2, Scoxc_new2, alpha)
dC3_new2<-dC3_frank(Scoxt_new2, Scoxc_new2, alpha)
}
else (stop("Error! Define a proper copula function"))
}
C_new2[ C_new2 == Inf ] <- 1
C_new2[ C_new2 == -Inf ] <- 0
dC_new2[,1:3][dC_new2[,1:3]<min.dC]<-min.dC
C_new2[C_new2<min.C]<-min.C
pllik_new2<-penlogreg_dep(del, eta, theta_new, gamma_new, h0t_new, h0c_new, regt_new, regc_new, Hcoxt_new2, Hcoxc_new2, C_new2, dC_new2, Rt, Rc, smpart, smparc)
incth<-theta_new-theta_old
incga<-gamma_new-gamma_old
ome<-0.6
while(pllik_new2<pllik_new1)
{
theta_new<-theta_old+ome*incth
gamma_new<-gamma_old+ome*incga
theta_new[theta_new<min.theta]<-min.theta
gamma_new[gamma_new<min.gamma]<-min.gamma
h0t_new<-baseHaz(theta_new, psix)
h0c_new<-baseHaz(gamma_new, psix)
h0t_new[h0t_new<min.ht]<-min.ht
h0c_new[h0c_new<min.hc]<-min.hc
H0t_new<-baseCumm(theta_new,Psix)
H0c_new<-baseCumm(gamma_new,Psix)
S0t_new<-exp(-H0t_new)
S0c_new<-exp(-H0c_new)
Hcoxt_new2<-H0t_new*eregt_new
Hcoxc_new2<-H0c_new*eregc_new
Scoxt_new2<-exp(-Hcoxt_new2)
Scoxc_new2<-exp(-Hcoxc_new2)
Scoxt_new2[Scoxt_new2<min.St]<-min.St
Scoxc_new2[Scoxc_new2<min.Sc]<-min.Sc
{
if(copula=='independent')
{
C_new2<-IndependentCopula(Scoxt_new2, Scoxc_new2)
dC_new2<-dC_ind(Scoxt_new2, Scoxc_new2)
dC3_new2<-dC3_ind(Scoxt_new2, Scoxc_new2)
}
else if(copula=='clayton')
{
C_new2<-ClaytonCopula(Scoxt_new2, Scoxc_new2, alpha)
dC_new2<-dC_clay(Scoxt_new2, Scoxc_new2, alpha)
dC3_new2<-dC3_clay(Scoxt_new2, Scoxc_new2, alpha)
}
else if(copula=='gumbel')
{
C_new2<-GumbelCopula(Scoxt_new2, Scoxc_new2, alpha)
dC_new2<-dC_gumbel(Scoxt_new2, Scoxc_new2, alpha)
dC3_new2<-dC3_gumbel(Scoxt_new2, Scoxc_new2, alpha)
}
else if (copula=='frank')
{
C_new2<-FrankCopula(Scoxt_new2, Scoxc_new2, alpha)
dC_new2<-dC_frank(Scoxt_new2, Scoxc_new2, alpha)
dC3_new2<-dC3_frank(Scoxt_new2, Scoxc_new2, alpha)
}
else (stop("Error! Define a proper copula function"))
}
C_new2[ C_new2 == Inf ] <- 1
C_new2[ C_new2 == -Inf ] <- 0
dC_new2[,1:3][dC_new2[,1:3]<min.dC]<-min.dC
C_new2[C_new2<min.C]<-min.C
pllik_new2<-penlogreg_dep(del, eta, theta_new, gamma_new, h0t_new, h0c_new, regt_new, regc_new, Hcoxt_new2, Hcoxc_new2, C_new2, dC_new2, Rt, Rc, smpart, smparc)
if(ome >= 1e-2) ome=ome*0.6
else if(ome < 1e-2 & ome >= 1e-5) ome=ome*5e-2
else if(ome < 1e-5 & ome >= 1e-30) ome=ome*1e-5
else
{
print("Iterations do not move baseline hazards")
break
}
}
h0_old<-c(theta_old, gamma_old)
h0_new<-c(theta_new, gamma_new)
pllik[k+1]<-pllik_new2
#print(k)
if (max(abs(h0_new-h0_old))<tol.thga & max(abs(coef_new-coef_old))<tol.bph)
{
break
}
else
{
theta_old<-theta_new
gamma_old<-gamma_new
h0t_old<-h0t_new
h0c_old<-h0c_new
beta_old<-beta_new
phi_old<-phi_new
H0t_old<-H0t_new
H0c_old<-H0c_new
S0t_old<-S0t_new
S0c_old<-S0c_new
regt_old<-regt_new
regc_old<-regc_new
eregt_old<-eregt_new
eregc_old<-eregc_new
Hcoxt_old<-Hcoxt_new2
Hcoxc_old<-Hcoxc_new2
Scoxt_old<-Scoxt_new2
Scoxc_old<-Scoxc_new2
C_old<-C_new2
dC_old<-dC_new2
dC3_old<-dC3_new2
pllik_old<-pllik_new2
}
}
fb<-apply(tcova, 2, mean)%*%beta_new
fph<-apply(tcova, 2, mean)%*%phi_new
theta_cop<-theta_new%*%exp(-fb)
gamma_cop<-gamma_new%*%exp(-fph)
theta_cop[theta_cop<min.theta]<-min.theta
gamma_cop[gamma_cop<min.gamma]<-min.gamma
h0t_cop<-baseHaz(theta_cop, psix)
h0c_cop<-baseHaz(gamma_cop, psix)
h0t_cop[h0t_cop<min.ht]<-min.ht
h0c_cop[h0c_cop<min.hc]<-min.hc
beta_cop<-beta_new
phi_cop<-phi_new
H0t_cop<-baseCumm(theta_cop,Psix)
H0c_cop<-baseCumm(gamma_cop,Psix)
S0t_cop<-exp(-H0t_cop) #the corresonding MPL estimates of S_{0T}(x)
S0c_cop<-exp(-H0c_cop) #the corresonding MPL estimates of S_{0C}(x)
regt<-tcova%*%beta_cop
regc<-tcova%*%phi_cop
eregt<-exp(regt)
eregc<-exp(regc)
hcoxt<-h0t_cop*eregt #the corresonding MPL estimates of h_{T}(x)
hcoxc<-h0c_cop*eregc #the corresonding MPL estimates of h_{C}(x)
Hcoxt<-H0t_cop*eregt #the corresonding MPL estimates of H_{T}(x)
Hcoxc<-H0c_cop*eregc #the corresonding MPL estimates of H_{C}(x)
Scoxt<-exp(-Hcoxt)
Scoxc<-exp(-Hcoxc)
Scoxt[Scoxt<min.St]<-min.St
Scoxc[Scoxc<min.Sc]<-min.Sc
{
if(copula=='independent')
{
dC3 <- dC3_ind(Scoxt, Scoxc)
dC <- dC_ind(Scoxt, Scoxc)
Co <- IndependentCopula(Scoxt, Scoxc)
}
else if(copula=='clayton')
{
dC3 <- dC3_clay(Scoxt, Scoxc, alpha)
dC <- dC_clay(Scoxt, Scoxc, alpha)
Co <- ClaytonCopula(Scoxt, Scoxc, alpha)
}
else if(copula=='gumbel')
{
dC3 <- dC3_gumbel(Scoxt, Scoxc, alpha)
dC <- dC_gumbel(Scoxt, Scoxc, alpha)
Co <- GumbelCopula(Scoxt, Scoxc, alpha)
}
else if (copula=='frank')
{
dC3 <- dC3_frank(Scoxt, Scoxc, alpha)
dC <- dC_frank(Scoxt, Scoxc, alpha)
Co <- FrankCopula(Scoxt, Scoxc, alpha)
}
else (stop("Error! Define a proper copula function"))
Co[ Co == Inf ] <- 1
Co[ Co == -Inf ] <- 0
dC[,1:3][dC[,1:3]<min.dC]<-min.dC
Co[Co<min.C]<-min.C
}
Ubeta_cop<-U_beta(p, cova, Hcoxt, Scoxt, del, eta, dC, Co) #score function respect to beta
Uphi_cop<-U_phi(p, cova, Hcoxc, Scoxc, del, eta, dC, Co) #score function respect to phi
Uth_cop<-U_theta(psix, Psix, theta_cop, h0t_cop, eregt, Scoxt, del, eta, dC, Co, smpart, Rt) #score function respect to theta
Uga_cop<-U_gamma(psix, Psix, gamma_cop, h0c_cop, eregc, Scoxc, del, eta, dC, Co, smparc, Rc)
if(asy==0) #without compute the asymptotic standard deviation of the MPL estimates
{
std_b_asy<-0;std_ph_asy<-0;std_th_asy<-0;std_ga_asy<-0;
cov_asy<-0
zp_b_ph<-0
std_ht0_asy<-0
std_hc0_asy<-0
std_Ht0_asy<-0
std_Hc0_asy<-0
std_St0_asy<-0
std_Sc0_asy<-0
}
else
{
asymstd<-asycoxphdepcencopmpl(Co, dC, dC3,
tcova,
psix, Psix,
theta_cop, gamma_cop, Uth_cop, Uga_cop,
h0t_cop, h0c_cop,
eregt, eregc,
Hcoxt, Scoxt, Hcoxc, Scoxc,
del, eta,
smpart, Rt, smparc, Rc,
eps, mid,
ac.theta, ac.gamma, ac.Utheta, ac.Ugamma)
std_b_asy<-asymstd$asympstd_beta #asymptotic std of beta
std_ph_asy<-asymstd$asympstd_phi #asymptotic std of phi
std_th_asy<-asymstd$asympstd_theta #asymptotic std of theta
std_ga_asy<-asymstd$asympstd_gamma #asymptotic std of gamma
cov_asy<-asymstd$asympcov_est #asymptotic covaraince of the estimates
z_b_asy<-beta_cop/std_b_asy #z score of beta
z_ph_asy<-phi_cop/std_ph_asy #z score of phi
p_b_asy<-2*pnorm(-abs(z_b_asy)) #p-value of beta
p_ph_asy<-2*pnorm(-abs(z_ph_asy)) #p-value of phi
zp_b_ph<-cbind(
c( beta_cop, phi_cop ),
c( std_b_asy, std_ph_asy ),
c( z_b_asy, z_ph_asy ),
c( p_b_asy, p_ph_asy )
) #regression ooefficient estimates with the corresponding standard deviations, z scores and p values
colnames(zp_b_ph)<-c("Est", "SD", "Z score", "P value")
rownames(zp_b_ph)<-rep(colnames(cova), 2)
std_ht0_asy<-sqrt( diag( psix%*%( cov_asy[(2*p+1):(2*p+m),(2*p+1):(2*p+m)] )%*%t( psix ) ) )
std_hc0_asy<-sqrt( diag( psix%*%( cov_asy[(2*p+m+1):(2*p+2*m),(2*p+m+1):(2*p+2*m)] )%*%t( psix ) ) )
std_Ht0_asy<-sqrt( diag( Psix%*%( cov_asy[(2*p+1):(2*p+m),(2*p+1):(2*p+m)] )%*%t( Psix ) ) )
std_Hc0_asy<-sqrt( diag( Psix%*%( cov_asy[(2*p+m+1):(2*p+2*m),(2*p+m+1):(2*p+2*m)] )%*%t( Psix ) ) )
std_St0_asy<-S0t_cop*std_Ht0_asy
std_Sc0_asy<-S0c_cop*std_Hc0_asy
}
if(ac==0){aic<-0}
else
{maxpl<-max(pllik)
aic<-2*2*(m+p)-2*maxpl}
if( identical(cova, matrix(0, n, 1)) ){noX <- TRUE}else{noX <- FALSE}
pos<-c( if(noX){rep(FALSE, 2*p)}else{rep(TRUE, 2*p)}, c(theta_cop>=ac.theta, gamma_cop>=ac.gamma) | c(Uth_cop>=ac.Utheta, Uga_cop>=ac.Ugamma) )
H<-hesscoxphcopmpl(p, Co, dC, dC3, tcova, h0t_cop, h0c_cop, eregt, eregc, Hcoxt, Scoxt, Hcoxc, Scoxc, psix, Psix, del, eta, 0, Rt, 0, Rc)
Hp<-hesscoxphcopmpl(p, Co, dC, dC3, tcova, h0t_cop, h0c_cop, eregt, eregc, Hcoxt, Scoxt, Hcoxc, Scoxc, psix, Psix, del, eta, smpart, Rt, smparc, Rc)
diag(Hp)<-diag(Hp)+eps
l<-penlogreg_dep(del, eta, theta_cop, gamma_cop, h0t_cop, h0c_cop, regt, regc, Hcoxt, Hcoxc, Co, dC, Rt, Rc, 0, 0)
if(cv==0){cvl<-0}
else
{
Hpinv<-matrix(0, (2*p+2*m), (2*p+2*m))
Hpinv[pos,pos]<-solve(Hp[pos, pos])
cvl<-l-sum(diag(Hpinv%*%H))
}
if(smpart_init=='REML')
{
HtRinv <- matrix(0, 2*p + 2*m, 2*p + 2*m)
sigmat2_old <- 1/( smpart)
sigmac2_old <- 1/( smparc)
diag(H)<-diag(H)+eps
HtRinv[pos, pos]<-solve( -H[pos, pos] + ( 1/sigmat2_old )*Rt_star[pos, pos] + ( 1/sigmac2_old )*Rc_star[pos, pos] )
dft<-m - ( sum( diag( HtRinv%*%Rt_star ) )/sigmat2_old )
sigmat2<-c( t( theta_cop )%*%Rt%*%theta_old/( dft ) )
smpart_old <- smpart
smpart <- 1/( sigmat2 )
}
else
{
sigmat2_old <- 1/( smpart)
sigmat2<-sigmat2_old
smpart_old <- smpart
smpart <- 1/( sigmat2 )
dft=m
}
if(smparc_init=='REML')
{
HcRinv <- matrix(0, 2*p + 2*m, 2*p + 2*m)
sigmat2_old <- 1/( smpart)
sigmac2_old <- 1/( smparc)
HcRinv[pos, pos] <- solve( -H[pos, pos] + ( 1/sigmat2_old )*Rt_star[pos, pos] + ( 1/sigmac2_old )*Rc_star[pos, pos] )
dfc <- m - ( sum( diag( HcRinv%*%Rc_star ) )/sigmac2_old )
sigmac2 <- c( t( gamma_cop )%*%Rc%*%gamma_old/( dfc ) )
smparc_old <- smparc
smparc <- 1/( sigmac2 )
}
else
{
sigmac2_old <- 1/( smparc)
sigmac2 <- sigmac2_old
smparc_old <- smparc
smparc <- 1/( sigmac2 )
dfc <- m
}
full.iter <- full.iter + k
if(cat.smpar=='Yes'){
cat(k2,"\tsmpart =",smpart,"\tsmparc =",smparc,"\tdft=",dft,"\tdfc=",dfc,"\titer =", k, "\n")
}
#print(k2)
if( abs( smpart - smpart_old )<( tol.smpar ) & abs( smparc - smparc_old ) < ( tol.smpar ) )
{
break
}
else{
dftold=dft
dfcold=dfc
smpart_old=smpart
smparc_old=smparc
sigmat2_old = sigmat2
sigmac2_old = sigmac2
}
}
iter <- c(k2, full.iter, k)
if(k2==1)
{
smpart <- smpart_old
smparc <- smparc_old
}
fit <- list(
mpl_theta=theta_cop, mpl_gamma=gamma_cop,
mpl_Utheta=Uth_cop, mpl_Ugamma=Uga_cop,
mpl_h0t=h0t_cop, mpl_h0c=h0c_cop,
mpl_H0t=H0t_cop, mpl_H0c=H0c_cop,
mpl_S0t=S0t_cop, mpl_S0c=S0c_cop,
mpl_beta=beta_cop, mpl_phi=phi_cop,
mpl_Ubeta=Ubeta_cop, mpl_Uphi=Uphi_cop,
mpl_ht=hcoxt, mpl_hc=hcoxc,
mpl_Ht=Hcoxt, mpl_Hc=Hcoxc,
mpl_St=Scoxt, mpl_Sc=Scoxc,
penloglik=pllik, iteration=iter,
smootht=smpart, smoothc=smparc, sigmat2=sigmat2, sigmac2=sigmac2,
mpl_cvl=cvl, mpl_aic=aic,
mpl_beta_sd=std_b_asy, mpl_phi_sd=std_ph_asy, mpl_h0t_sd=std_th_asy, mpl_h0c_sd=std_ga_asy,
mpl_ht0_sd=std_ht0_asy, mpl_hc0_sd=std_hc0_asy,
mpl_Ht0_sd=std_Ht0_asy, mpl_Hc0_sd=std_Hc0_asy,
mpl_St0_sd=std_St0_asy, mpl_Sc0_sd=std_Sc0_asy,
mpl_est_cov=cov_asy, mpl_beta_phi_zp=zp_b_ph,
binedg=binedg, binwv=binwv,ID=ID, psix=psix, Psix=Psix,
numIntKnt = numIntKnt, #length(binedg)-2
IntKnt = IntKnt, #binedg[2:(numIntKnt+1)]
bryKnt = bryKnt, #c(0, max(binedg))
numSp = numSp, #numIntKnt+ordSp
surv=surv,
cova=cova,
ordSp=ordSp,
binCount=binCount,
tie=tie,
tau=tau,
copula=copula,
pent=pent,
smpart_init=smpart_init,
penc=penc,
smparc_init=smparc_init,
maxit2=maxit2,
maxit=maxit,
mid=mid,
asy=asy,
ac=ac,
cv=cv,
ac.theta=ac.theta,
ac.gamma=ac.gamma,
ac.Utheta=ac.Utheta,
ac.Ugamma=ac.Ugamma,
min.theta=min.theta,
min.gamma=min.gamma,
min.ht=min.ht,
min.hc=min.hc,
min.St=min.St,
min.Sc=min.Sc,
min.C=min.C,
min.dC=min.dC,
eps=eps,
tol.thga=tol.thga,
tol.bph=tol.bph,
tol.smpar=tol.smpar,
cat.smpar=cat.smpar
)
class(fit) = "coxph_mpl_dc"
return( fit )
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.