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.
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:
- Estimate the parameters of the marginals
- 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]