R code for fitting a model to unbalanced longitudinal data with a copula

In my previous blog post I showed how to fit a model to longitudinal data with a copula.

Remember that the focus in my previous post was on balanced longitudinal data. However, the following R code demonstrates how to fit a copula when dealing with unbalanced longitudinal data.

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

library(copula)
library(SPREDA)


#R functions for fitting a copula to unbalanced longitudinal data
#can be downloaded from: 
source("http://www.datall-analyse.nl/R/copula_unbalanced.R")



#EXAMPLE 1: Gaussian copula

##Generate some artificial data using a n-dimensional Gaussian copula

#n-dimensional Gaussian copula
ndim <- 20 #number of dimensions
myCop <- normalCopula(param=.6, dim=ndim, dispstr="ex")

#Marginals:
#all n marginals follow a Gumbel distribution with location=0 and scale=.1
myMvdc <- mvdc(copula=myCop, rep("myGumbel", ndim),
               rep(list(list(location=0, scale=.1)), ndim))

#The Gumbel distribution is skewed to the left
plot(1, type="n", xlab="x", ylab="f(x)", xlim=c(-.5, .2), ylim=c(0, 4))
curve(dmyGumbel(x, location=0, scale=.1), add=TRUE)

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

#Inspect the data from the first 5 cases, and observe that there is considerable
#individual heterogeneity (that is, values within the same case tend
#to be more alike than values across cases) and serial dependence (that is,
#values at different time points within the same case are correlated).
matplot(t(rs[1:5,]), type="b", pch=1, lty=1, ylab="y(t)", xlab="t-index")

#Next, transform the generated data by employing a nonlinear model.
#The employed nonlinear model is an asymptotic function through the origin.
plot(1, type="n", xlab="x", ylab="f(x)", xlim=c(0, 6000), ylim=c(-1.4, 0))
curve(-1.4*(1-exp(-exp(-7.5)*x)), add=TRUE) #nonlinear function

#Transform the data
t <- rep(seq(250, 250*ndim, 250), times=ns) #time points
phi1 <- -1.4 #asymptote
x1 <- rep(c(-2, 0, 2), each=(ns/3)*ndim) #covariate values
phi2 <- -7.5 + .7*x1 #rate of change depending on covariate x1
ys <- phi1*(1-exp(-exp(phi2)*t)) + as.vector(t(rs)) #response
dataL <- data.frame(id=rep(1:ns, each=ndim), ys , time=t, x1) #combine data

#Visualize the transformed data
matplot(t(matrix(as.vector(dataL[,2]), ncol=ndim, byrow=TRUE)), type="b", pch=1,
        lty=1, col=rep(c("red", "darkgreen", "blue"), each=ns/3),
        ylab="y(t)", xlab="t-index")
legend("topright", lty=c(1, 1, 1), pch=c(1,1,1),
       col=c("red", "darkgreen", "blue"),
       legend=c("x1 = -2", "x1 = 0", "x1 = 2"),
       bty="n", y.intersp=1.2, cex=.7)

#Inspect 5 cases with value x1=2, and notice the individual heterogeneity
#and serial dependence
dataLselect <- dataL[which(x1==2),]
matplot(t(matrix(as.vector(dataLselect[,2]), ncol=ndim, byrow=TRUE))[,1:5],
        type="b", pch=1, lty=1, ylab="y(t)", xlab="t-index")


#Create an unbalanced data set
#x=-2, keep all observations
dataL1 <-dataL[which(dataL$x1==-2), ]

#x=0, discard observations with time > 3750
dataL2 <-dataL[which(dataL$x1==0), ]
dataL2 <- dataL2[-which(dataL2$time>3750),]

#x=2, discard observations with time > 2500
dataL3 <-dataL[which(dataL$x1==2), ]
dataL3 <- dataL3[-which(dataL3$time>2500),]

#combine the data
dataLu <- rbind(dataL1, dataL2, dataL3)
dataLuWide <- reshape(data=dataLu, v.names="ys", timevar="time",
                      idvar="id", drop="x1", direction="wide")[,-1]

cNames <- colnames(dataLuWide)
orderCols <- order(nchar(cNames), cNames)
dataLuWide <- dataLuWide[, orderCols]

#Visualize the unbalanced data
matplot(t(dataLuWide), type="b", pch=1,
        lty=1, col=rep(c("red", "darkgreen", "blue"), each=ns/3),
        ylab="y(t)", xlab="t-index")
legend("topright", lty=c(1, 1, 1), pch=c(1,1,1),
       col=c("red", "darkgreen", "blue"),
       legend=c("x1 = -2", "x1 = 0", "x1 = 2"),
       bty="n", y.intersp=1.2, cex=.7)


##Fit a n-dimensional 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 of the IFM: fit the marginals, and obtain the marginal parameters
#For obtaining the marginal parameters, we will fit a nonlinear model to the data

#Use the getInitial function for SSasympOrig() for obtaining start values
#for phi1 and beta0 (note: phi1=Asym and beta0=lrc)
(start <- getInitial(ys ~ SSasympOrig(time, Asym, lrc), data = dataLu))

#Define the log-likelihood function for the nonlinear model
llik <- function (theta, y, t, cov) {
  #parameters
  phi1 <- theta[1]
  beta0 <- theta[2]
  beta1 <- theta[3]
  sigma <- exp(theta[4])
  
  phi2 <- beta0 + beta1*cov
  mu <- phi1*(1-exp(-exp(phi2)*t))
  
  #log-likelihood
  sum(log( dsev((y-mu) / sigma) / sigma ))
}

#For obtaining a start value for sigma (that is, the scale of the Gumbel
#distribution), set beta1=0, and inspect the resulting log-likelihood function
sigmas <- seq(.001, 5, length.out=100)
startSigma <- sapply(1:length(sigmas),
                     function (i) llik(theta=c(start[1], start[2], 0, log(sigmas[i])),
                                       y=dataLu$ys, t=dataLu$time, cov=dataLu$x1))

#Plot the resulting log-likelihood function, and notice that the maximum of
#the log-likelihood is at about .5. Hence, a plausible start value for sigma is .5
plot(sigmas, startSigma, xlab="sigma", ylab="log-likelihood", type="l")

#Now that we have plausible start values, fit the nonlinear model by
#maximizing the log-likehood function.
#Note that the log transformation of sigma enforces the scale parameter
#of the Gumbel distribution to be positive.
nlmod <- optim(c(start[1], start[2], 0, log(.5)), llik,
               control=list(fnscale=-1),
               y=dataLu$ys, t=dataLu$time, cov=dataLu$x1)
nlmod #results
nlmod$par #MLEs (the true values were -1.4, -7.5, .7, log(.1)=-2.3, respectively)
exp(nlmod$par[4]) #the MLE for scale parameter

#Based on the estimated parameters, subsequently compute the cumulative
#probabilities for the marginals
phi1hat <- nlmod$par[1]; beta0hat <- nlmod$par[2]
beta1hat <- nlmod$par[3]; sigmahat <- exp(nlmod$par[4])
phi2hat <- cbind(1, dataLu$x1)%*%c(beta0hat, beta1hat)
muhat <- phi1hat*(1-exp(-exp(phi2hat)*dataLu$time))

udat <- pmyGumbel(dataLu$ys, location=muhat, scale=sigmahat)
udatDf <- cbind(dataLu, udat)
udatWide <- reshape(data=udatDf, v.names="udat", timevar="time",
                    idvar="id", drop=c("x1", "ys"), direction="wide")[,-1]

cNames <- colnames(udatWide)
orderCols <- order(nchar(cNames), cNames)
udatWide <- as.matrix(udatWide[, orderCols])


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

#Starting value for the copula parameter
corMat <- sin(cor(udatWide, method="kendall", use="complete.obs")*pi/2)
(tauEst <- mean(corMat[upper.tri(corMat)]))


#Fit a Gaussian copula with an exchangeable dispersion matrix. Note that is
#also possible to fit other dispersion matrices (that is, AR(1) and Toeplitz).

myCopUn <- normalCopula(param=.5, dim=ncol(udatWide), dispstr="ex")

fitGaussian <- optim(tauEst, loglikUnbNormal, method="BFGS",
                     hessian=TRUE, control=list(fnscale=-1),
                     copula=myCopUn, u=udatWide)

#Did the optimization converge? 
fitGaussian$convergence # 0 indicates successful convergence
#MLE for the copula parameter (remember: the true value for the parameter was .6)
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)


#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.
#You may use a resampling method (such as a parametric bootstrap method)
#for estimating standard errors that include the variability of marginal MLEs.





#EXAMPLE 2: Gaussian copula with a Toeplitz dispersion matrix

##Generate some artificial data using a n-dimensional Gaussian copula

#n-dimensional Gaussian copula
ndim <- 20 #number of dimensions
myCop <- normalCopula(param=c(.6, .3, rep(0, 17)), dim=ndim, dispstr="toep")
#Elements in this Toeplitz dispersion matrix with a horizontal distance of 3
#or more elements from the diagonal are zero. This yields the following matrix:
toeplitz(c(1, c(.6, .3, rep(0, 17))))

#Marginals:
#all n marginals follow a Gumbel distribution with location=0 and scale=.7
myMvdc <- mvdc(copula=myCop, rep("myGumbel", ndim),
               rep(list(list(location=0, scale=.1)), ndim))

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

#Transform the data
t <- rep(seq(250, 250*ndim, 250), times=ns) #time points
phi1 <- -1.4 #asymptote
x1 <- rep(c(-2, 0, 2), each=(ns/3)*ndim) #covariate values
phi2 <- -7.5 + .7*x1 #rate of change depending on covariate x1
ys <- phi1*(1-exp(-exp(phi2)*t)) + as.vector(t(rs)) #response
dataL <- data.frame(id=rep(1:ns, each=ndim), ys , time=t, x1) #combine data

#Visualize the transformed data
matplot(t(matrix(as.vector(dataL[,2]), ncol=ndim, byrow=TRUE)), type="b", pch=1,
        lty=1, col=rep(c("red", "darkgreen", "blue"), each=ns/3),
        ylab="y(t)", xlab="t-index")
legend("topright", lty=c(1, 1, 1), pch=c(1,1,1),
       col=c("red", "darkgreen", "blue"),
       legend=c("x1 = -2", "x1 = 0", "x1 = 2"),
       bty="n", y.intersp=1.2, cex=.7)

#Create unbalanced data by removing 300 observations at random
nrs <- sample(1:nrow(dataL), 300)
dataLu <- dataL[-nrs, ]

#Visualize the unbalanced data
dataLuWide <- reshape(data=dataLu, v.names="ys", timevar="time",
                      idvar="id", drop="x1", direction="wide")[,-1]

cNames <- colnames(dataLuWide)
orderCols <- order(nchar(cNames), cNames)
dataLuWide <- dataLuWide[, orderCols]

matplot(t(dataLuWide), type="b", pch=1,
        lty=1, col=rep(c("red", "darkgreen", "blue"), each=ns/3),
        ylab="y(t)", xlab="t-index")
legend("topright", lty=c(1, 1, 1), pch=c(1,1,1),
       col=c("red", "darkgreen", "blue"),
       legend=c("x1 = -2", "x1 = 0", "x1 = 2"),
       bty="n", y.intersp=1.2, cex=.7)


##Fit a n-dimensional 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 of the IFM: fit the marginals, and obtain the marginal parameters
#For obtaining the marginal parameters, we will fit a nonlinear model to the data

#Use the getInitial function for SSasympOrig() for obtaining start values
#for phi1 and beta0 (note: phi1=Asym and beta0=lrc)
(start <- getInitial(ys ~ SSasympOrig(time, Asym, lrc), data = dataLu))

#Define the log-likelihood function for the nonlinear model
llik <- function (theta, y, t, cov) {
  #parameters
  phi1 <- theta[1]
  beta0 <- theta[2]
  beta1 <- theta[3]
  sigma <- exp(theta[4])
  
  phi2 <- beta0 + beta1*cov
  mu <- phi1*(1-exp(-exp(phi2)*t))
  
  #log-likelihood
  sum(log( dsev((y-mu) / sigma) / sigma ))
}

#For obtaining a start value for sigma (that is, the scale of the Gumbel
#distribution), set beta1=0, and inspect the resulting log-likelihood function
sigmas <- seq(.001, 5, length.out=100)
startSigma <- sapply(1:length(sigmas),
                     function (i) llik(theta=c(start[1], start[2], 0, log(sigmas[i])),
                                       y=dataLu$ys, t=dataLu$time, cov=dataLu$x1))

#Plot the resulting log-likelihood function, and notice that the maximum of
#the log-likelihood is at about .5. Hence, a plausible start value for sigma is .5
plot(sigmas, startSigma, xlab="sigma", ylab="log-likelihood", type="l")

#Now that we have plausible start values, fit the nonlinear model by
#maximizing the log-likehood function.
#Note that the log transformation of sigma enforces the scale parameter
#of the Gumbel distribution to be positive.
nlmod <- optim(c(start[1], start[2], 0, log(.5)), llik,
               control=list(fnscale=-1),
               y=dataLu$ys, t=dataLu$time, cov=dataLu$x1)
nlmod #results
nlmod$par #MLEs (the true values were -1.4, -7.5, .7, log(.1)=-2.3, respectively)
exp(nlmod$par[4]) #the MLE for scale parameter

#Based on the estimated parameters, subsequently compute the cumulative
#probabilities for the marginals
phi1hat <- nlmod$par[1]; beta0hat <- nlmod$par[2]
beta1hat <- nlmod$par[3]; sigmahat <- exp(nlmod$par[4])
phi2hat <- cbind(1, dataLu$x1)%*%c(beta0hat, beta1hat)
muhat <- phi1hat*(1-exp(-exp(phi2hat)*dataLu$time))

udat <- pmyGumbel(dataLu$ys, location=muhat, scale=sigmahat)
udatDf <- cbind(dataLu, udat)
udatWide <- reshape(data=udatDf, v.names="udat", timevar="time",
                    idvar="id", drop=c("x1", "ys"), direction="wide")[,-1]

cNames <- colnames(udatWide)
orderCols <- order(nchar(cNames), cNames)
udatWide <- as.matrix(udatWide[, orderCols])


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

#Starting value for the copula parameters
corMat <- sin(cor(udatWide, method="kendall", use="pairwise.complete.obs")*pi/2)
round(corMat, 2) #have a closer look at the correlation matrix
upper1 <- sapply(1:(nrow(corMat)-1), function (i) corMat[i, i+1])
upper2 <- sapply(1:(nrow(corMat)-2), function (i) corMat[i, i+2])
(tauEst <- c(mean(upper1), mean(upper2)))


#Fit a Gaussian copula with a Toeplitz dispersion matrix.

#Note: we will fix elements of the Toeplitz matrix with a horizontal distance
#of 3 or more elements from the diagonal to zero.
myCopUn <- normalCopula(param=fixParam(c(.5, .5, rep(0, 17)),
                                       c(FALSE, FALSE, rep(TRUE, 17))),
                        dim=ncol(udatWide), dispstr="toep")

fitGaussian <- optim(c(tauEst), loglikUnbNormal, method="BFGS",
                     hessian=TRUE, control=list(fnscale=-1),
                     copula=myCopUn, u=udatWide)

#Did the optimization converge? 
fitGaussian$convergence # 0 indicates successful convergence
#MLE for the copula parameters
#(remember: the true values for the parameters were .6 and .3, respectively)
fitGaussian$par
#Standard errors for MLE copula parameters
SEparms <- sqrt(diag(solve(-1*fitGaussian$hessian)))
#Construct 95% confidence interval for the copula parameters
fitGaussian$par[1] + c(1,-1)*qnorm(.05/2)*SEparms[1] #parameter 1 (true value: .6)
fitGaussian$par[2] + c(1,-1)*qnorm(.05/2)*SEparms[2] #parameter 2 (true value: .3)

#Remember that when using the two-stage parametric ML method, the standard errors
#of the MLE for the copula parameters are usually underestimated (i.e., too small).





##EXAMPLE 3: ARCHIMEDEAN COPULAS

##Generate some artificial data using a n-dimensional Gumbel copula

#n-dimensional Gumbel copula
ndim <- 20 #number of dimensions
myCop <- gumbelCopula(param=3, dim=ndim)

#Marginals:
#all n marginals follow a Gumbel distribution with location=0 and scale=.7
myMvdc <- mvdc(copula=myCop, rep("myGumbel", ndim),
               rep(list(list(location=0, scale=.1)), ndim))

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

#Inspect the data from the first 5 cases, and observe that there is considerable
#individual heterogeneity (that is, values within the same case tend
#to be more alike than values across cases) and serial dependence (that is,
#values at different time points within the same case are correlated).
matplot(t(rs[1:5,]), type="b", pch=1, lty=1, ylab="y(t)", xlab="t-index")

#Next, transform the generated data by employing a nonlinear model.
#The employed nonlinear model is an asymptotic function through the origin.
t <- rep(seq(250, 250*ndim, 250), times=ns) #time points
phi1 <- -1.4 #asymptote
x1 <- rep(c(-2, 0, 2), each=(ns/3)*ndim) #covariate values
phi2 <- -7.5 + .7*x1 #rate of change depending on covariate x1
ys <- phi1*(1-exp(-exp(phi2)*t)) + as.vector(t(rs)) #response
dataL <- data.frame(id=rep(1:ns, each=ndim), ys , time=t, x1) #combine data

#Visualize the transformed data
matplot(t(matrix(as.vector(dataL[,2]), ncol=ndim, byrow=TRUE)), type="b", pch=1,
        lty=1, col=rep(c("red", "darkgreen", "blue"), each=ns/3),
        ylab="y(t)", xlab="t-index")
legend("topright", lty=c(1, 1, 1), pch=c(1,1,1),
       col=c("red", "darkgreen", "blue"),
       legend=c("x1 = -2", "x1 = 0", "x1 = 2"),
       bty="n", y.intersp=1.2, cex=.7)

#Inspect 5 cases with value x1=2, and notice the individual heterogeneity
#and serial dependence
dataLselect <- dataL[which(x1==2),]
matplot(t(matrix(as.vector(dataLselect[,2]), ncol=ndim, byrow=TRUE))[,1:5],
        type="b", pch=1, lty=1, ylab="y(t)", xlab="t-index")

#Create unbalanced data
#x=-2, keep all observations
dataL1 <-dataL[which(dataL$x1==-2), ]

#x=0, discard observations with time > 3750
dataL2 <-dataL[which(dataL$x1==0), ]
dataL2 <- dataL2[-which(dataL2$time>3750),]

#x=2, discard observations with time > 2500
dataL3 <-dataL[which(dataL$x1==2), ]
dataL3 <- dataL3[-which(dataL3$time>2500),]

#combine the data
dataLu <- rbind(dataL1, dataL2, dataL3)
dataLuWide <- reshape(data=dataLu, v.names="ys", timevar="time",
                      idvar="id", drop="x1", direction="wide")[,-1]

cNames <- colnames(dataLuWide)
orderCols <- order(nchar(cNames), cNames)
dataLuWide <- dataLuWide[, orderCols]

#Visualize the unbalanced data
matplot(t(dataLuWide), type="b", pch=1,
        lty=1, col=rep(c("red", "darkgreen", "blue"), each=ns/3),
        ylab="y(t)", xlab="t-index")
legend("topright", lty=c(1, 1, 1), pch=c(1,1,1),
       col=c("red", "darkgreen", "blue"),
       legend=c("x1 = -2", "x1 = 0", "x1 = 2"),
       bty="n", y.intersp=1.2, cex=.7)


##Fit a n-dimensional Gumbel 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 of the IFM: fit the marginals, and obtain the marginal parameters
#For obtaining the marginal parameters, we will fit a nonlinear model to the data

#Use the getInitial function for SSasympOrig() for obtaining start values
#for phi1 and beta0 (note: phi1=Asym and beta0=lrc)
(start <- getInitial(ys ~ SSasympOrig(time, Asym, lrc), data = dataLu))

#Define the log-likelihood function for the nonlinear model
llik <- function (theta, y, t, cov) {
  #parameters
  phi1 <- theta[1]
  beta0 <- theta[2]
  beta1 <- theta[3]
  sigma <- exp(theta[4])
  
  phi2 <- beta0 + beta1*cov
  mu <- phi1*(1-exp(-exp(phi2)*t))
  
  #log-likelihood
  sum(log( dsev((y-mu) / sigma) / sigma ))
}

#For obtaining a start value for sigma (that is, the scale of the Gumbel
#distribution), set beta1=0, and inspect the resulting log-likelihood function
sigmas <- seq(.001, 5, length.out=100)
startSigma <- sapply(1:length(sigmas),
                     function (i) llik(theta=c(start[1], start[2], 0, log(sigmas[i])),
                                       y=dataLu$ys, t=dataLu$time, cov=dataLu$x1))

#Plot the resulting log-likelihood function, and notice that the maximum of
#the log-likelihood is at about .5. Hence, a plausible start value for sigma is .5
plot(sigmas, startSigma, xlab="sigma", ylab="log-likelihood", type="l")

#Now that we have plausible start values, fit the nonlinear model by
#maximizing the log-likehood function.
#Note that the log transformation of sigma enforces the scale parameter
#of the Gumbel distribution to be positive.
nlmod <- optim(c(start[1], start[2], 0, log(.5)), llik,
               control=list(fnscale=-1),
               y=dataLu$ys, t=dataLu$time, cov=dataLu$x1)
nlmod #results
nlmod$par #MLEs (the true values were -1.4, -7.5, .7, log(.1)=-2.3, respectively)
exp(nlmod$par[4]) #the MLE for scale parameter

#Based on the estimated parameters, subsequently compute the cumulative
#probabilities for the marginals
phi1hat <- nlmod$par[1]; beta0hat <- nlmod$par[2]
beta1hat <- nlmod$par[3]; sigmahat <- exp(nlmod$par[4])
phi2hat <- cbind(1, dataLu$x1)%*%c(beta0hat, beta1hat)
muhat <- phi1hat*(1-exp(-exp(phi2hat)*dataLu$time))

udat <- pmyGumbel(dataLu$ys, location=muhat, scale=sigmahat)
udatDf <- cbind(dataLu, udat)
udatWide <- reshape(data=udatDf, v.names="udat", timevar="time",
                    idvar="id", drop=c("x1", "ys"), direction="wide")[,-1]

cNames <- colnames(udatWide)
orderCols <- order(nchar(cNames), cNames)
udatWide <- as.matrix(udatWide[, orderCols])


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

#Starting value for the copula parameter
corMat <- cor(udatWide, method="kendall", use="complete.obs")
tauEst <- mean(corMat[upper.tri(corMat)])
(thetaGumbel <- 1 / (1-tauEst))


#Fit a Gumbel copula
myCopUn <- gumbelCopula(param=2, dim=ncol(udatWide))

fitGumbel <- optim(thetaGumbel, loglikUnbArchi, method="BFGS",
                   hessian=TRUE, control=list(fnscale=-1),
                   copula=myCopUn, u=udatWide)

#Did the optimization converge? 
fitGumbel$convergence # 0 indicates successful convergence
#MLE for the copula parameter (remember: the true value for the parameter was 3)
fitGumbel$par
#Standard error for MLE copula parameter
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)

#Remember that 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).





##EXAMPLE 4: t-COPULA

##Generate some artificial data using a n-dimensional t-copula

#n-dimensional t-copula
ndim <- 20 #number of dimensions
myCop <- tCopula(param=.6, dim=ndim, df=4, dispstr="ex")

#Marginals:
#all n marginals follow a Gumbel distribution with location=0 and scale=.7
myMvdc <- mvdc(copula=myCop, rep("myGumbel", ndim),
               rep(list(list(location=0, scale=.1)), ndim))

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

#Next, transform the generated data by employing a nonlinear model.
t <- rep(seq(250, 250*ndim, 250), times=ns) #time points
phi1 <- -1.4 #asymptote
x1 <- rep(c(-2, 0, 2), each=(ns/3)*ndim) #covariate values
phi2 <- -7.5 + .7*x1 #rate of change depending on covariate x1
ys <- phi1*(1-exp(-exp(phi2)*t)) + as.vector(t(rs)) #response
dataL <- data.frame(id=rep(1:ns, each=ndim), ys , time=t, x1) #combine data

#Visualize the transformed data
matplot(t(matrix(as.vector(dataL[,2]), ncol=ndim, byrow=TRUE)), type="b", pch=1,
        lty=1, col=rep(c("red", "darkgreen", "blue"), each=ns/3),
        ylab="y(t)", xlab="t-index")
legend("topright", lty=c(1, 1, 1), pch=c(1,1,1),
       col=c("red", "darkgreen", "blue"),
       legend=c("x1 = -2", "x1 = 0", "x1 = 2"),
       bty="n", y.intersp=1.2, cex=.7)

#Inspect 5 cases with value x1=2, and notice the individual heterogeneity
#and serial dependence
dataLselect <- dataL[which(x1==2),]
matplot(t(matrix(as.vector(dataLselect[,2]), ncol=ndim, byrow=TRUE))[,1:5],
        type="b", pch=1, lty=1, ylab="y(t)", xlab="t-index")

#Create unbalanced data
#x=-2, keep all observations
dataL1 <-dataL[which(dataL$x1==-2), ]

#x=0, discard observations with time > 3750
dataL2 <-dataL[which(dataL$x1==0), ]
dataL2 <- dataL2[-which(dataL2$time>3750),]

#x=2, discard observations with time > 2500
dataL3 <-dataL[which(dataL$x1==2), ]
dataL3 <- dataL3[-which(dataL3$time>2500),]

#combine the data
dataLu <- rbind(dataL1, dataL2, dataL3)
dataLuWide <- reshape(data=dataLu, v.names="ys", timevar="time",
                      idvar="id", drop="x1", direction="wide")[,-1]

cNames <- colnames(dataLuWide)
orderCols <- order(nchar(cNames), cNames)
dataLuWide <- dataLuWide[, orderCols]

#Visualize the unbalanced data
matplot(t(dataLuWide), type="b", pch=1,
        lty=1, col=rep(c("red", "darkgreen", "blue"), each=ns/3),
        ylab="y(t)", xlab="t-index")
legend("topright", lty=c(1, 1, 1), pch=c(1,1,1),
       col=c("red", "darkgreen", "blue"),
       legend=c("x1 = -2", "x1 = 0", "x1 = 2"),
       bty="n", y.intersp=1.2, cex=.7)


##Fit a n-dimensional t-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 of the IFM: fit the marginals, and obtain the marginal parameters
#For obtaining the marginal parameters, we will fit a nonlinear model to the data

#Use the getInitial function for SSasympOrig() for obtaining start values
#for phi1 and beta0 (note: phi1=Asym and beta0=lrc)
(start <- getInitial(ys ~ SSasympOrig(time, Asym, lrc), data = dataLu))

#Define the log-likelihood function for the nonlinear model
llik <- function (theta, y, t, cov) {
  #parameters
  phi1 <- theta[1]
  beta0 <- theta[2]
  beta1 <- theta[3]
  sigma <- exp(theta[4])
  
  phi2 <- beta0 + beta1*cov
  mu <- phi1*(1-exp(-exp(phi2)*t))
  
  #log-likelihood
  sum(log( dsev((y-mu) / sigma) / sigma ))
}

#For obtaining a start value for sigma (that is, the scale of the Gumbel
#distribution), set beta1=0, and inspect the resulting log-likelihood function
sigmas <- seq(.001, 5, length.out=100)
startSigma <- sapply(1:length(sigmas),
                     function (i) llik(theta=c(start[1], start[2], 0, log(sigmas[i])),
                                       y=dataLu$ys, t=dataLu$time, cov=dataLu$x1))

#Plot the resulting log-likelihood function, and notice that the maximum of
#the log-likelihood is at about .5. Hence, a plausible start value for sigma is .5
plot(sigmas, startSigma, xlab="sigma", ylab="log-likelihood", type="l")

#Now that we have plausible start values, fit the nonlinear model by
#maximizing the log-likehood function.
#Note that the log transformation of sigma enforces the scale parameter
#of the Gumbel distribution to be positive.
nlmod <- optim(c(start[1], start[2], 0, log(.5)), llik,
               control=list(fnscale=-1),
               y=dataLu$ys, t=dataLu$time, cov=dataLu$x1)
nlmod #results
nlmod$par #MLEs (the true values were -1.4, -7.5, .7, log(.1)=-2.3, respectively)
exp(nlmod$par[4]) #the MLE for scale parameter

#Based on the estimated parameters, subsequently compute the cumulative
#probabilities for the marginals
phi1hat <- nlmod$par[1]; beta0hat <- nlmod$par[2]
beta1hat <- nlmod$par[3]; sigmahat <- exp(nlmod$par[4])
phi2hat <- cbind(1, dataLu$x1)%*%c(beta0hat, beta1hat)
muhat <- phi1hat*(1-exp(-exp(phi2hat)*dataLu$time))

udat <- pmyGumbel(dataLu$ys, location=muhat, scale=sigmahat)
udatDf <- cbind(dataLu, udat)
udatWide <- reshape(data=udatDf, v.names="udat", timevar="time",
                    idvar="id", drop=c("x1", "ys"), direction="wide")[,-1]

cNames <- colnames(udatWide)
orderCols <- order(nchar(cNames), cNames)
udatWide <- as.matrix(udatWide[, orderCols])


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

#Starting value for the copula parameter
corMat <- sin(cor(udatWide, method="kendall", use="pairwise.complete.obs")*pi/2)
(tauEst <- mean(corMat[upper.tri(corMat)]))


#Fit a t-copula with an exchangeable dispersion matrix. Note that is also possible
#to fit other dispersion matrices (that is, AR(1) and Toeplitz).

myCopUn <- tCopula(param=.5, df=4, df.fixed=TRUE, dim=ncol(udatWide), dispstr="ex")

#choose df=3 as a naïve start value for the degrees of freedom
fitT <- optim(c(tauEst), loglikUnbT, method="BFGS",
              hessian=TRUE, control=list(fnscale=-1),
              copula=myCopUn, u=udatWide)

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

#Remember that 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).


#Fit the t-copula again, but this time also estimate the degrees of freedom

myCopUn2 <- tCopula(param=.5, df=4, df.fixed=FALSE, dim=ncol(udatWide), dispstr="ex")

#choose 3 as a naïve start value for the degrees of freedom
fitT <- optim(c(tauEst, 3), loglikUnbT, method="BFGS",
              hessian=TRUE, control=list(fnscale=-1),
              copula=myCopUn2, u=udatWide)

#Did the optimization converge? 
fitT$convergence # 0 indicates successful convergence
#MLE for the copula parameters
fitT$par
#Standard errors for MLE copula parameters
SEparms <- sqrt(diag(solve(-1*fitT$hessian)))
#Construct 95% confidence intervals for the copula parameters
fitT$par[1] + c(1,-1)*qnorm(.05/2)*SEparms[1] #true value: .6
fitT$par[2] + c(1,-1)*qnorm(.05/2)*SEparms[2] #true value: 4

#Remember that when using the two-stage parametric ML method, the standard errors
#of the MLE for the copula parameters are usually underestimated (i.e., too small).