MGL-regression.RmdThis example implements the MGL and MGL-EV copula regression models of Zhengxiao Li et al. (2021, “MGL copulas”).
library(rMGLReg)
library(fitdistrplus)
library(splines)
library(snpar)
data("danishmulti")
dt <- data.table::data.table(danishmulti)
dtnew <- dt[Building>0&Contents>0]
y1 <- dtnew$Building
y2 <- dtnew$Contents
y <- cbind(y1, y2)
# empirical cdf
u1 <- snpar::kde(y[,1], kernel = "gaus",
xgrid = y[,1],
h = 0.2)$Fhat
u2 <- snpar::kde(y[,2], kernel = "gaus",
xgrid = y[,2],
h = 0.2)$Fhat
U <- cbind(u1, u2) # two-dim
library(ggplot2)
library(latex2exp)
Usample <- U
XY <- y
newtheme <- theme_bw() + theme(axis.line = element_line(colour = "black"),
axis.text.x = element_text(margin = margin(t = 0.25, unit = "cm")),
axis.text.y = element_text(margin = margin(r = 0.25, unit = "cm"),
size = 10,
vjust = 0.5,
hjust = 0.5),
axis.ticks.length = unit(-0.1, "cm"),
plot.title = element_text(hjust = 0.5),
legend.direction = 'vertical',
legend.position = c(.95, .99),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.text = element_text(size = 10))
dtplot <- data.frame(U1 = Usample[,1], U2 = Usample[,2],
logY1 = log(XY[,1]), logY2 = log(XY[,2]))
p1 <- ggplot(data = dtplot, mapping = aes(y = logY1,
x = logY2)) +
newtheme +
geom_point() +
labs(title = "",
x = TeX("$log Y_1$"),
y = TeX("$log Y_2$")) +
scale_x_continuous(limits = c(-7.5, 5),
breaks = seq(-7.5, 5, by = 2.5)) +
scale_y_continuous(limits = c(-5, 5),
breaks = seq(-5, 5, by = 2.5))
p1
p2 <- ggplot(data = dtplot, mapping = aes(y = U2,
x = U1)) +
newtheme +
geom_point() +
labs(title = "",
x = TeX("$u_1$"),
y = TeX("$u_2$")) +
scale_x_continuous(limits = c(0, 1),
breaks = seq(0, 1, by = 0.2)) +
scale_y_continuous(limits = c(0, 1),
breaks = seq(0, 1, by = 0.2))
p2
library(patchwork)
p0 <- p1 + p2 + plot_layout(ncol = 2)
p0
dtnew[, year := as.numeric(substr(Date, start = 1, stop = 4))]
dtnew[, x1 := (year - mean(year))/(sd(year))]
dtnew[, day := difftime(as.Date(Date), as.Date("1980-01-03"), units="days")]
X <- splines::ns(dtnew$year, knots = quantile(dtnew$year, c(0.5)), intercept = T)
# survival MGL copual regression
m.MGL180 <- MGL.reg(U = U, copula = "MGL180",
X = X, method = "Nelder-Mead",
initpar = c(-0.32, 0.001, 0.001)
)
# survival MGL-EV copual regression
m.MGLEV180 <- MGL.reg(U = U, copula = "MGL-EV180",
method = "Nelder-Mead",
X = X,
initpar = c(-0.32, 0.001, 0.001)
)
# gumbel regression
m.gumbel <- MGL.reg(U = U, copula = "Gumbel",
method = "Nelder-Mead",
X = X,
initpar = c(-0.32, 0.001, 0.001)
)
# MGB2 regression
m.MGB2 <- MGL.reg(U = U, copula = "MGB2",
method = "Nelder-Mead",
X = X,
initpar = c(0.1, -0.5, -0.5, -0.5, -0.5)
)
estimates.copula <- rbind(c(m.MGL180$estimates, NA, NA),
c(m.MGLEV180$estimates, NA, NA),
c(m.gumbel$estimates, NA, NA),
c(m.MGB2$estimates)
)
sd.copula <- rbind(c(m.MGL180$se, NA, NA),
c(m.MGLEV180$se, NA, NA),
c(m.gumbel$se, NA, NA),
c(m.MGB2$se)
)
ll.copula <- rbind((m.MGL180)$loglike, (m.MGLEV180)$loglike, (m.gumbel)$loglike, (m.MGB2)$loglike)
AIC.copula <- rbind((m.MGL180)$AIC, (m.MGLEV180)$AIC, (m.gumbel)$AIC, (m.MGB2)$AIC)
BIC.copula <- rbind((m.MGL180)$BIC, (m.MGLEV180)$BIC, (m.gumbel)$BIC, (m.MGB2)$BIC)
table1 <- cbind(estimates.copula, sd.copula, ll.copula, AIC.copula, BIC.copula)
row.names(table1) <- c('MGL180', 'MGLEV180', 'Gumbel', 'MGB2')
colnames(table1) <- c("estimates", "estimates","estimates", "estimates","estimates",
"se", "se","se","se","se",
'loglike',
'AIC', 'BIC')
knitr::kable(t(table1), digits = 3)
# Survival MGL regression model
par.est <- m.MGL180$estimates
agepred <- seq(from = 1980, to = 1990)
Xcopula <- splines::ns(agepred, knots = quantile(dtnew$year, c(0.5)), intercept = T)
delta.est <- exp(Xcopula%*%par.est)
cov.est <- -solve(m.MGL180$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}
matplot(delta.mat, col = 'gray', type = 'l', ylim = c(0, 2), main = 'Survival MGL')
lines(delta.est, col = 'red', lwd = 2)
# Surivial MGL-EV regression
par.est <- m.MGLEV180$estimates
delta.est <- exp(Xcopula%*%par.est)
cov.est <- -solve(m.MGLEV180$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}
matplot(delta.mat, col = 'gray', type = 'l', ylim = c(0, 2), main = 'Survival MGL-EV')
lines(delta.est, col = 'red', lwd = 2)
# Gumbel regression
par.est <- m.gumbel$estimates
delta.est <- exp(Xcopula%*%par.est) + 1
cov.est <- -solve(m.gumbel$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,]) + 1
}
matplot(delta.mat, col = 'gray', type = 'l', ylim = c(0, 2), main = 'Gumbel')
lines(delta.est, col = 'red', lwd = 2)
# MGB2 regression
par.est <- m.MGB2$estimates
delta.est <- exp(Xcopula%*%par.est[1:ncol(X)])
cov.est <- -solve(m.MGB2$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 5)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
pars.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
beta.sim <- pars.sim[,1:ncol(X)]
for(i in 1:100){
delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}
matplot(delta.mat, col = 'gray', type = 'l', ylim = c(0, 4), main = 'MGB2')
lines(delta.est, col = 'red', lwd = 2)