R code for fitting a copula to censored data

The following R code fits a bivariate (Archimedean or elliptical) copula to data where one of the variables contains censored observations. The censored observations can be left, right or interval censored.

copula-censored-data

Two-stage parametric ML method
The copula is fitted using the two-stage parametric ML approach (also referred to as the Inference Functions for Margins [IFM] method). This method fits a copula in two steps:

  1. Estimate the parameters of the marginals
  2. Fix the marginal parameters to the values estimated in first step, and subsequently estimate the copula parameters.

The two-stage parametric ML method is described by Joe in his 2014 book Dependence modeling with copulas (pp. 16-17).

Suggestions and/or questions? Please contact Stefan Gelissen (email: info at datall-analyse.nl).
See this page for an overview of all of Stefan’s
R code blog posts.

R code (last update: 2017-05-17)

library(copula)
library(MASS)
library(survival)
library(SPREDA)

##---Define functions
#All of the following functions are needed for fitting a copula
#to censored data by means of the two-stage parametric ML estimation method.
#The two-stage parametric ML estimation method is also referred to as the 
#method of Inference Functions for Margins (IFM).

#Function for computing the log-likelihood
loglikCensored <- function (param, copula, u, which.censored=NULL,
                            which.censored2=NULL,
                            event=rep(1, length(x))){
  
  #-param is a vector of parameter values.
  #-copula is a copula object. At the moment, only Archimedean and elliptical
  # copulas are accepted.
  #-u is a n x d matrix, where n are the number of observations (i.e, sample size),
  # and d the number of dimensions. In case of interval censored observations,
  # the dimensions of this matrix become n x (d+1).
  # The values in the u-matrix are all defined in [0, 1].
  #-which.censored is a number specifying in which column of the u-matrix
  # the censored observations can be found.
  #-which.censored2 is an optional number that is specified only in case of interval
  # censoring. More specifically, in case of interval censoring, this number
  # indicates in which column of the u-matrix the upper value of the censoring
  # interval can be found.
  #-event is the censoring indicator (0=right censored, 1=uncensored,
  # 2=left censored, 3=interval censored).
  
  if (copula@dimension>2) {
    stop("Only bivariate copulas are allowed")
  }
  if (is.null(attr(copula@parameters, "fixed"))) {
    copula@parameters <- param
    cop.param <- getTheta(copula, freeOnly=FALSE, attr=TRUE)
  }
  else {
    copula@parameters[!attr(copula@parameters, "fixed")] <- param
    cop.param <- getTheta(copula, freeOnly=TRUE, attr=TRUE)}

  lower <- attr(cop.param, "param.lowbnd")
  upper <- attr(cop.param, "param.upbnd")
  admissible <- !(any(is.na(cop.param) | cop.param > upper | 
                        cop.param < lower))
  
  if (admissible) {
    #compute the log-likelihood
    lik <- ifelse(
      #right censoring
      event==0, dRight(u, copula, which.censored, which.censored2), 
      #uncensored observations
      ifelse(
        event==1, dUncensored(u, copula, which.censored, which.censored2),
        #left censoring
        ifelse(
          event==2, dLeft(u, copula, which.censored, which.censored2),
          #interval censoring
          dInterval(u, copula, which.censored, which.censored2) ))) 
    
    sum(log(lik))}
  else NaN
}

#Function for computing the conditional probability of right censored observations
dRight <- function (u, copula, which.censored, which.censored2) {
  
  ndim <- copula@dimension
  
  if(is.vector(u)){
    u <- matrix(u, nrow=1)
  }
  
  if(is.null(which.censored2)) {
    uu <- u[, -which.censored] #uncensored variable
    uc <- u[, which.censored] #censored variable
  }
  else {
    uu <- u[, -c(which.censored, which.censored2)] #uncensored variable
    uc <- u[, which.censored] #censored variable
  }
  
  #compute the survival probabilities
  su <- 1-uu; sc <- 1-uc
  
  #conditional probability (using the survival copula)
  1 - matrix(cCopula(cbind(uu, 1-sc), copula), ncol=ndim)[,ndim]
}

#Function for computing the conditional probability of left censored observations
dLeft <- function (u, copula, which.censored, which.censored2) {
  
  ndim <- copula@dimension
  
  if(is.vector(u)){
    u <- matrix(u, nrow=1)
  }
  
  if(is.null(which.censored2)) {
    uu <- u[, -which.censored] #uncensored variable
    uc <- u[, which.censored] #censored variable
  }
  else {
    uu <- u[, -c(which.censored, which.censored2)] #uncensored variable
    uc <- u[, which.censored] #censored variable
  }
  
  #conditional probability
  matrix(cCopula(cbind(uu, uc), copula), ncol=ndim)[,ndim]
}

#Function for computing the conditional probability of interval censored observations
dInterval <- function (u, copula, which.censored, which.censored2) {
  
  ndim <- copula@dimension
  
  if(is.vector(u)){
    u <- matrix(u, nrow=1)
  }
  
  if(is.null(which.censored2)) {
    NA
  }
  else {
    #assign some arbitrary value (within [0, 1]) to observations having NA
    u[which(is.na(u[, which.censored2])), which.censored2] <- 0
    
    uu <- u[, -c(which.censored,  which.censored2)] #uncensored variable
    ucl <- u[, which.censored] #lower bound censoring interval
    ucu <- u[, which.censored2] #upper bound censoring interval
    
    #conditional probability
    matrix(cCopula(cbind(uu, ucu), copula), ncol=ndim)[,ndim] -
      matrix(cCopula(cbind(uu, ucl), copula), ncol=ndim)[,ndim]
  }
}

#Function for computing the joint density of uncensored observations
dUncensored <- function (u, copula, which.censored, which.censored2) {
  
  if(is.vector(u)){
    u <- matrix(u, nrow=1)
  }
  
  if(is.null(which.censored2)) {
    dCopula(u, copula)
  }
  else {
    u <- u[, -which.censored2]
    dCopula(u, copula)
  }
}

#Functions for the cdf, pdf, quantile function, and random number generation
#of a Weibull distribution having a location scale parametrization.
pmyWeibull <- function(x, location=0, scale=1, threshold=0) {
  xt <- (log(x-threshold)-location)/scale
  psev(xt)
}

dmyWeibull <- function(x, location=0, scale=1, threshold=0) {
  xt <- (log(x-threshold)-location)/scale
  dsev(xt)/(scale*x)
}

qmyWeibull <- function(p, location=0, scale=1, threshold=0) {
  exp(location + qsev(p)*scale) + threshold
}

rmyWeibull <- function(n, location=0, scale=1, threshold=0) {
  exp(location + rsev(n)*scale) + threshold
}

##---End of function definitions





##EXAMPLE 1: right censored observations

##Generate some artificial data using a bivariate Gaussian copula

#Bivariate Gaussian copula
myCop <- normalCopula(param=-.7, dim=2)

#Marginals:
#the first marginal is a normal distribution, and
#the second marginal is a Weibull distribution
myMvdc <- mvdc(copula=myCop, c("norm", "myWeibull"),
               list(list(mean=5, sd=.8), list(location=2.4, scale=1.3)))

#Visualize the copula
persp(myMvdc, dMvdc, xlim=c(0, 8), ylim=c(0, 110))

#Draw sample (this sample will subsequently be used as the artificial data)
ns <- 200 #sample size
xys <- rMvdc(ns, myMvdc)

#Inspect the generated data
plot(xys[,1], xys[,2], xlab="x", ylab="y")

#Apply random censoring to the observations of variable y. The censored
#observations will subsequently be treated as right censored cases.
censored <- rmyWeibull(ns, location=3.4, scale=.8)
events <- ifelse(xys[,2] > censored, 0, 1)
xys[,2] <- pmin(xys[,2], censored)

#Plot the censored data
plot(xys[,1], xys[,2], xlab="x", ylab="y", col=events+1)
legend("topright", pch=c(1, 1), col=c(1, 2),
       legend=c("censored", "uncensored"),
       bty="n", y.intersp=1.2, cex=.7)

#Censoring distribution: 0=right censored, 1=uncensored
table(events)



##Fit a bivariate Gaussian copula to the artificial data
#by means of the two-stage parametric ML estimation method, also referred
#to as the method of Inference Functions for Margins (IFM).

#Stage 1: fit the marginals, and obtain the marginal parameters
(normMLE <- fitdistr(xys[,1], densfun="normal"))
normMean1 <- normMLE$estimate[1]
normSD1 <- normMLE$estimate[2]

(weibullMLE <- survreg(Surv(time=xys[,2], events) ~ 1, dist="weibull"))
weibullLocation2 <- coef(weibullMLE)
weibullScale2 <- weibullMLE$scale

#Compute the cumulative probabilities for the marginals
udat <- cbind(pnorm(xys[,1], mean=normMean1, sd=normSD1),
              pmyWeibull(xys[,2], location=weibullLocation2, scale=weibullScale2))


#Stage 2: fix the parameters of the marginals and estimate the copula parameter

#Starting value for the copula parameter
(tauEst <- sin(cor(xys[events==1, 1], xys[events==1, 2], method="kendall")*pi/2))

fitGaussian <- optim(tauEst, loglikCensored, method="BFGS",
                     hessian=TRUE,
                     control=list(fnscale=-1), copula=myCop,
                     u=udat, which.censored=2, event=events)

#Did the optimization converge? 
fitGaussian$convergence # 0 indicates successful convergence
#MLE for the copula parameter
fitGaussian$par
#Standard error for MLE copula parameter
SEparms <- sqrt(diag(solve(-1*fitGaussian$hessian)))
#Construct 95% confidence interval for the copula parameter
(ci <- fitGaussian$par + c(1,-1)*qnorm(.05/2)*SEparms)

#Visually inspect the log-likelihood function used for estimating
#the copula parameter
taus <- seq(-1, 0, length.out=100)
liks <- sapply(taus, function(k) loglikCensored(param=k, copula=myCop, u=udat,
                                                which.censored=2, event=events))
plot(taus, liks, type="l", xlab="Copula parameter", ylab="log-likelihood")
#Include the MLE for the copula parameter
abline(v=fitGaussian$par, col="red", lty=2)
#Include the confidence interval for the copula parameter
abline(v=ci, col="darkgreen", lty=2)
#Include the true value for the copula parameter
abline(v=myCop@parameters, col="blue")
#Add a legend
legend("bottomright", lty=c(2, 2, 1), col=c("red", "darkgreen", "blue"),
       legend=c("MLE", "95% confidence interval", "true value"),
       bty="n", y.intersp=1.2, cex=.7)


#When using the two-stage parametric ML method, the standard error of the MLE
#for the copula parameter is usually underestimated (i.e., too small).
#The reason for this underestimation is that the marginal MLEs were fixed at
#the second stage of the fitting process. As a consequence, when estimating
#the copula parameter at this second stage, variation in the MLEs of the
#marginals is not taken into account.
#However, a resampling method (such as a nonparametric bootstrap method)
#will yield standard errors that include the variability of marginal MLEs.

#Function for performing the nonparametric bootstrap method
bootIFM <- function (x, event) {
  
  #bootstrap sample
  bsample <- sample(1:nrow(x), replace=TRUE)
  bx <- x[bsample, ]
  bevent <- event[bsample]
  
  #stage 1: fit the marginals
  normMLEb <- fitdistr(bx[,1], densfun="normal")
  normMean1b <- normMLEb$estimate[1]
  normSD1b <- normMLEb$estimate[2]
  
  weibullMLEb <- survreg(Surv(time=bx[,2], bevent) ~ 1, dist="weibull")
  weibullLocation2b <- coef(weibullMLEb)
  weibullScale2b <- weibullMLEb$scale
  
  #Cumulative probabilities for marginals
  udatb <- cbind(pnorm(bx[,1], mean=normMean1b, sd=normSD1b),
                 pmyWeibull(bx[,2], location=weibullLocation2b, scale=weibullScale2b))
  
  #stage 2: fit the copula
  fitCop <- optim(fitGaussian$par, loglikCensored, method="BFGS",
                  control=list(fnscale=-1), copula=myCop,
                  u=udatb, which.censored=2, event=bevent)
  
  #return MLEs for the bootstrap sample
  bootResult <- rep(NA, 5)
  bootResult <- c(normMean1b, normSD1b, weibullLocation2b, weibullScale2b,
                  fitCop$par)
  names(bootResult) <- c("normMean1", "normSD1",
                         "weibullLocation2", "weibullScale2", "theta")
  bootResult
}

#Number of bootstrap samples
nboot <- 1000
#Perform the bootstrap method for the two-stage parametric ML method.
#Note that this bootstrap can be computationally expensive.
bootMLE <- replicate(nboot, bootIFM(x=xys, event=events))


#Compare the original standard error for the copula parameter with that
#of the bootstrap method.

#First we extract the MLEs estimated for the original sample
originalMLE <- c(normMean1, normSD1, weibullLocation2, weibullScale2,
                 fitGaussian$par)
names(originalMLE) <- c("normMean1", "normSD1",
                        "weibullLocation2", "weibullScale2", "theta")

#Compute the standard errors for the bootstrap
(SEparmsBoot <- sqrt(diag(((bootMLE-originalMLE)%*%t((bootMLE-originalMLE)))/nboot)))

SEparms #standard error of the copula parameter for the original sample
SEparmsBoot[5] #standard error of the copula parameter for the bootstrap
#Notice that the standard error obtained with the bootstrap is larger
#than that of the original sample.


#Compute a bias-corrected percentile confidence interval for the copula
#parameter (based on the bootstrap results)
q0 <- mean(bootMLE[5,] <= originalMLE[5])
alpha <- .05
zl.bias2 <- round(nboot*pnorm(2*qnorm(q0)+qnorm(alpha/2)))
zu.bias2 <- round(nboot*pnorm(2*qnorm(q0)+qnorm(1-alpha/2)))
#95% confidence interval
(ciBoot <- c(sort(bootMLE[5,])[zl.bias2], sort(bootMLE[5,])[zu.bias2]))

#Visually Inspect (as above) the log-likelihood function used for estimating
#the copula parameter
plot(taus, liks, type="l", xlab="Copula parameter", ylab="log-likelihood")
abline(v=fitGaussian$par, col="red", lty=2)
#Include the 95% bootstap confidence interval
abline(v=ciBoot, col="darkgreen", lty=2)
#Include the true copula parameter
abline(v=myCop@parameters, col="blue")
#Add a legend
legend("bottomright", lty=c(2, 2, 1), col=c("red", "darkgreen", "blue"),
       legend=c("MLE", "95% confidence interval", "true value"),
       bty="n", y.intersp=1.2, cex=.7)





##EXAMPLE 2: arbitrary censoring

##Generate some artificial data using a bivariate Gumbel copula

#Bivariate Gumbel copula
myCop <- gumbelCopula(param=3, dim=2)

#Marginals:
#the first marginal is a normal distribution, and
#the second marginal is a Weibull distribution
myMvdc <- mvdc(copula=myCop, c("norm", "myWeibull"),
               list(list(mean=5, sd=.8), list(location=2.4, scale=1.3)))

#Visualize the copula
persp(myMvdc, dMvdc, xlim=c(0, 8), ylim=c(0, 110))

#Draw sample (this sample will subsequently be used as the artificial data)
ns <- 200 #sample size
xys <- rMvdc(ns, myMvdc)

#Inspect data
plot(xys[,1], xys[,2], xlab="x", ylab="y")


#Apply censoring: mixture of left, right, and interval censoring

#Generate intervals for censoring
intCens <- seq(0,70,10)

#Assign observations to censoring intervals
assignCens <- findInterval(xys[,2], intCens)

#Observations with >100 are right censored and get value NA for time2
intCensRC <- c(intCens, NA)

#Lower bounds of censoring interval
L <- sapply(assignCens, function (k) intCensRC[k])

#Upper bounds of censoring interval
U <- sapply(assignCens, function (k) intCensRC[k+1])

#Indicator for type of censoring
Cens <- ifelse(is.na(U), 0, 3)

#Censoring in the interval [0, 10] is identical to left censoring at t=10
Cens <- ifelse(L==0, 2, Cens) #set censoring indicator to 2 (=left censoring)
#For these left censored observations, set L to 100, and U to NA
L <- ifelse(Cens==2, U, L)
U <- ifelse(Cens==2, NA, U)

#Inspect data
head(cbind(L, U, xys[,2], Cens), n=15)

#Censoring distribution: 0=right censored, 2=left censored, 3=interval censored
table(Cens)



##Fit a bivariate Gumbel copula to the artificial data
#by means of the two-stage parametric ML estimation method

#Stage 1: fit the marginals, and obtain the marginal parameters
(normMLE <- fitdistr(xys[,1], densfun="normal"))
normMean1 <- normMLE$estimate[1]
normSD1 <- normMLE$estimate[2]

(weibullMLE <- survreg(Surv(time=L, time2=U, event=Cens, type="interval") ~ 1,
                       dist="weibull"))
weibullLocation2 <- coef(weibullMLE)
weibullScale2 <- weibullMLE$scale

#Cumulative probabilities for marginals
udat <- cbind(pnorm(xys[,1], mean=normMean1, sd=normSD1),
              pmyWeibull(L, location=weibullLocation2, scale=weibullScale2),
              pmyWeibull(U, location=weibullLocation2, scale=weibullScale2))


#Stage 2: fix the parameters of the marginals and estimate the copula parameter

#Starting value for the copula parameter
corEst <- cor(xys[Cens==3,1], U[Cens==3], method="kendall")
(tauEst <- 1 / (1-corEst))

fitGumbel <- optim(tauEst, loglikCensored, method="BFGS",
                   hessian=TRUE,
                   control=list(fnscale=-1), copula=myCop,
                   u=udat, which.censored=2, which.censored2=3,
                   event=Cens)

#Starting value for the copula parameter
fitGumbel$convergence # 0 indicates successful convergence
#MLE for the copula parameter
fitGumbel$par
#Standard error for MLE copula parameter: remember that for the two-stage
#method this standard error is usually too small
SEparms <- sqrt(diag(solve(-1*fitGumbel$hessian)))
#Construct 95% confidence interval for the copula parameter
(ci <- fitGumbel$par + c(1,-1)*qnorm(.05/2)*SEparms)

#Visually inspect the log-likelihood function used for estimating
#the copula parameter
taus <- seq(1, 5, length.out=100)
liks <- sapply(taus, function(k) loglikCensored(param=k, copula=myCop, u=udat,
                                                which.censored=2, which.censored2=3,
                                                event=Cens))
plot(taus, liks, type="l", xlab="Copula parameter", ylab="log-likelihood")
#Include the MLE for the copula parameter
abline(v=fitGumbel$par, col="red", lty=2)
#Include the confidence interval for the copula parameter
abline(v=ci, col="darkgreen", lty=2)
#Include the true value for the copula parameter
abline(v=myCop@parameters, col="blue")
#Add a legend
legend("bottomleft", lty=c(2, 2, 1), col=c("red", "darkgreen", "blue"),
       legend=c("MLE", "95% confidence interval", "true value"),
       bty="n", y.intersp=1.2, cex=.7)

In the next example we will use the full-likelihood for fitting
a bivariate copula to data having right censored observations.

##---Define functions
#All of the following functions are needed for fitting a bivariate copula
#to right censored data by means of the full-likelihood.

#Function for computing the log-likelihood
loglikMvdcCensored <- function (param, mvdc, x, which.censored=NULL,
                                event=rep(1, length(x))){
  
  #-param is a vector of parameter values.
  #-mvdc is a mvdc object. At the moment, only Archimedean and elliptical
  # copulas are accepted.
  #-x is a n x d matrix, where n are the number of observations (i.e, sample size),
  # and d the number of dimensions.
  #-which.censored is a number specifying in which column of the x-matrix
  # the censored observations can be found.
  #-event is the censoring indicator (0=right censored, 1=uncensored).
  
  #set the parameters for the copula
  mvdc <- copula:::setMvdcPar(mvdc, param, noCheck=TRUE)
  
  #compute the log-likelihood
  lik <- ifelse(
    #right censoring
    event==0, dMvdcRight(x, mvdc, which.censored), 
    #uncensored observations
    dMvdc(x, mvdc))
  
  sum(log(lik))
}

#Function for computing the density of right censored observations
dMvdcRight <- function (x, mvdc, which.censored) {
  if (is.vector(x)) 
    x <- matrix(x, nrow = 1)
  
  i <- which.censored #censored variable
  j <- ifelse(i==1, 2, 1) #uncensored variable
  
  u <- x
  
  for (k in 1:2) {
    cdf.expr <- copula:::asCall(paste0("p", mvdc@margins[k]),
                                mvdc@paramMargins[[k]])
    u[, k] <- eval(cdf.expr, list(x = x[, k]))}
  
  pdf.expr <- copula:::asCall(paste0("d", mvdc@margins[j]),
                              mvdc@paramMargins[[j]])
  fx <- eval(pdf.expr, list(x = x[, j]))
  
  #conditional probability (using the survival copula)
  (1 - matrix(cCopula(cbind(u[, j], u[, i]), mvdc@copula), ncol=2)[,2]) * fx  
}

#Functions for the cdf, pdf, quantile function, and random number generation
#of a Weibull distribution having a location scale parametrization.
pmyWeibull <- function(x, location=0, scale=1, threshold=0) {
  xt <- (log(x-threshold)-location)/scale
  psev(xt)
}

dmyWeibull <- function(x, location=0, scale=1, threshold=0) {
  xt <- (log(x-threshold)-location)/scale
  dsev(xt)/(scale*x)
}

qmyWeibull <- function(p, location=0, scale=1, threshold=0) {
  exp(location + qsev(p)*scale) + threshold
}

rmyWeibull <- function(n, location=0, scale=1, threshold=0) {
  exp(location + rsev(n)*scale) + threshold
}

##---End of function definitions



##EXAMPLE 3: full-likelihood for right censored observations

##Generate some artificial data using a Gumbel copula

#Bivariate Gumbel copula
myCop <- gumbelCopula(param=3, dim=2)

#Marginals:
#the first marginal is a normal distribution, and
#the second marginal is a Weibull distribution
myMvdc <- mvdc(copula=myCop, c("norm", "myWeibull"),
               list(list(mean=5, sd=.8), list(location=2.4, scale=1.3)))

#Draw sample (this sample will subsequently be used as the artificial data)
ns <- 200 #sample size
xys <- rMvdc(ns, myMvdc)

#Inspect data
plot(xys[,1], xys[,2], xlab="x", ylab="y")

#Apply random censoring to the observations of variable y. The censored
#observations will subsequently be treated as right censored cases.
censored <- rmyWeibull(ns, location=3.4, scale=.8)
events <- ifelse(xys[,2] > censored, 0, 1)
xys[,2] <- pmin(xys[,2], censored)

#Censoring distribution: 0=right censored, 1=uncensored
table(events) #censoring



##Fit a bivariate Gaussian copula to the artificial data

#Obtain start values for fitting the copula
normMLE <- fitdistr(xys[,1], densfun="normal")
normMean1 <- normMLE$estimate[1]
normSD1 <- normMLE$estimate[2]

weibullMLE <- survreg(Surv(time=xys[,2], events) ~ 1, dist="weibull")
weibullLocation2 <- coef(weibullMLE)
weibullScale2 <- weibullMLE$scale

#Starting value for the copula parameter
corEst <- cor(xys[events==1,1], xys[events==1,2], method="kendall")
(tauEst <- 1 / (1-corEst))

#Fit a Gumbel copula to the artificial data using the full-likelihood
myGumbel <- mvdc(copula=gumbelCopula(param=2, dim=2),
                 c("norm", "myWeibull"),
                 paramMargins=list(list(mean=0, sd=1),
                                   list(location=0, scale=1)))

start <- c(normMean1, normSD1, weibullLocation2, weibullScale2, tauEst)
start

fitGumbel <- optim(start, loglikMvdcCensored, method="L-BFGS-B",
                   lower=c(-Inf, 0, -Inf, 0, 0),
                   upper=c(Inf, Inf, Inf, Inf, Inf),
                   hessian=TRUE,
                   control=list(fnscale=-1), mvdc=myGumbel,
                   x=xys, which.censored=2, event=events)

#MLEs (true values for the parameters are 5, .8, 2.4, 1.3, 3, respectively)
fitGumbel$par
#Converged?
fitGumbel$convergence

#Standard errors for MLE parameters
SEparms <- sqrt(diag(solve(-1*fitGumbel$hessian)))

#Construct some confidence intervals:
#95% confidence interval for the Weibull location parameter (true value is 2.4)
fitGumbel$par[3] + c(1,-1)*qnorm(.05/2)*SEparms[3]
#95% confidence interval for the copula parameter (true value is 3)
fitGumbel$par[5] + c(1,-1)*qnorm(.05/2)*SEparms[5]




##Generate some new artificial data, but this time using a Gaussian copula

#Bivariate Gaussian copula
myCop <- normalCopula(param=-.7, dim=2)

#Marginals:
#the first marginal is a normal distribution, and
#the second marginal is a Weibull distribution
myMvdc <- mvdc(copula=myCop, c("norm", "myWeibull"),
               list(list(mean=5, sd=.8), list(location=2.4, scale=1.3)))

#Draw sample (this sample will subsequently be used as the artificial data)
ns <- 200 #sample size
xys <- rMvdc(ns, myMvdc)

#Inspect data
plot(xys[,1], xys[,2], xlab="x", ylab="y")

#Apply random censoring to the observations of variable y. The censored
#observations will subsequently be treated as right censored cases.
censored <- rmyWeibull(ns, location=3.4, scale=.8)
events <- ifelse(xys[,2] > censored, 0, 1)
xys[,2] <- pmin(xys[,2], censored)

#Censoring distribution: 0=right censored, 1=uncensored
table(events) #censoring



##Fit a bivariate Gaussian copula to the artificial data

#Obtain start values for fitting the copula
normMLE <- fitdistr(xys[,1], densfun="normal")
normMean1 <- normMLE$estimate[1]
normSD1 <- normMLE$estimate[2]

weibullMLE <- survreg(Surv(time=xys[,2], events) ~ 1, dist="weibull")
weibullLocation2 <- coef(weibullMLE)
weibullScale2 <- weibullMLE$scale

#Starting value for the copula parameter
(tauEst <- sin(cor(xys[events==1, 1], xys[events==1, 2], method="kendall")*pi/2))


#Fit a Gaussian copula to the artificial data using the full-likelihood
myGaussian <- mvdc(copula=normalCopula(param=-1, dim=2),
                   c("norm", "myWeibull"),
                   paramMargins=list(list(mean=0, sd=1),
                                     list(location=0, scale=1)))

start <- c(normMean1, normSD1, weibullLocation2, weibullScale2, tauEst)
start

fitGaussian <- optim(start, loglikMvdcCensored, method="L-BFGS-B",
                     lower=c(-Inf, 0, -Inf, 0, -1),
                     upper=c(Inf, Inf, Inf, Inf, 1),
                     hessian=TRUE,
                     control=list(fnscale=-1), mvdc=myGaussian,
                     x=xys, which.censored=2, event=events)

#MLEs (true values for the parameters are 5, .8, 2.4, 1.3, -.7, respectively)
fitGaussian$par
#converged?
fitGaussian$convergence

#Standard errors for MLE parameters
SEparms <- sqrt(diag(solve(-1*fitGaussian$hessian)))

#Construct some confidence intervals:
#95% confidence interval for the Weibull location parameter (true value is 2.4)
fitGaussian$par[3] + c(1,-1)*qnorm(.05/2)*SEparms[3]
##95% confidence interval for the copula parameter (true value is -.7)
fitGaussian$par[5] + c(1,-1)*qnorm(.05/2)*SEparms[5]