This file contains the source code of an exemplary application of the D-vine copula based quantile regression approach implemented in the R-package vinereg and presented in Kraus and Czado (2017): D-vine copula based quantile regression, Computational Statistics and Data Analysis, 110, 1-18. Please, feel free to address questions to daniel.kraus@tum.de.
library(vinereg)
library(ggplot2)
library(dplyr)
library(tidyr)
set.seed(5)
We consider the data set abalone from the UCI Machine Learning Repository (https://archive.ics.uci.edu/ml/datasets/abalone) and focus on the female
sub-population. In a first application we only consider continuous variables and
fit models to predict the quantiles of weight (whole) given the predictors
length, diameter, and height.
data(abalone, package = "PivotalR")
abalone_f <- abalone %>%
dplyr::filter(sex == "F") %>% # select female abalones
dplyr::select(-id, -sex) %>% # remove id and sex variables
dplyr::filter(height < max(height)) # remove height outlier
pairs(abalone_f, pch = ".")
We consider the female subset and fit a parametric regression D-vine for the response weight given the covariates len, diameter and height (ignoring the discreteness of some of the variables). The D-vine based model is selected sequentially by maximizing the conditional log-likelihood of the response given the covariates. Covariates are only added if they increase the (possibly AIC- or BIC-corrected) conditional log-likelihood.
We use the function vinereg() to fit a regression D-vine for predicting the
response weight given the covariates length, diameter, and height. The
argument family_set determines how the pair-copulas are estimated. We will
only use one-parameter families and the t copula here. The
selcrit argument specifies the penalty type for the conditional
log-likelihood criterion for variable selection.
fit_vine_par <- vinereg(
whole ~ length + diameter + height,
data = abalone_f,
family_set = c("onepar", "t"),
selcrit = "aic"
)
The result has a field order that shows the selected covariates and their
ranking order in the D-vine.
fit_vine_par$order
## [1] "length" "height" "diameter"
The field vine contains the fitted D-vine, where the first variable
corresponds to the response. The object is of class "vinecop_dist" so we can
use rvineocpulib's functionality to summarize the model
summary(fit_vine_par$vine)
## # A data.frame: 6 x 11
## tree edge conditioned conditioning var_types family rotation parameters
## 1 1 1, 2 c,c gumbel 180 5.2
## 1 2 2, 4 c,c gumbel 180 2.4
## 1 3 4, 3 c,c gumbel 180 2.4
## 2 1 1, 4 2 c,c t 0 0.45, 21.87
## 2 2 2, 3 4 c,c t 0 0.91, 4.68
## 3 1 1, 3 4, 2 c,c t 0 0.31, 7.40
## df tau loglik
## 1 0.81 1629
## 1 0.58 672
## 1 0.59 709
## 2 0.30 147
## 2 0.73 1193
## 2 0.20 72
We can also plot the contours of the fitted pair-copulas.
contour(fit_vine_par$vine)
In order to visualize the predicted influence of the covariates, we plot the estimated quantiles arising from the D-vine model at levels 0.1, 0.5 and 0.9 against each of the covariates.
# quantile levels
alpha_vec <- c(0.1, 0.5, 0.9)
We call the fitted() function on fit_vine_par to extract the fitted values
for multiple quantile levels. This is equivalent to predicting the quantile at
the training data. The latter function is more useful for out-of-sample
predictions.
pred_vine_par <- fitted(fit_vine_par, alpha = alpha_vec)
# equivalent to:
# predict(fit_vine_par, newdata = abalone.f, alpha = alpha_vec)
head(pred_vine_par)
## 0.1 0.5 0.9
## 1 0.6717860 0.7667223 0.8716322
## 2 0.6968823 0.7903849 0.9014626
## 3 0.6794480 0.7843661 0.8914415
## 4 0.7754026 0.8800504 0.9976871
## 5 0.5950430 0.6996904 0.8273507
## 6 0.6755811 0.7729168 0.8860188
To examine the effect of the individual variables, we will plot the predicted quantiles against each of the variables. To visualize the relationship more clearly, we add a smoothed line for each quantile level. This gives an estimate of the expected effect of a variable (taking expectation with respect to all other variables).
plot_effects(fit_vine_par)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
The fitted quantile curves suggest a non-linear effect of all three variables.
This can be compared to linear quantile regression:
pred_lqr <- pred_vine_par
for (a in seq_along(alpha_vec)) {
my.rq <- quantreg::rq(
whole ~ length + diameter + height,
tau = alpha_vec[a],
data = abalone_f
)
pred_lqr[, a] <- quantreg::predict.rq(my.rq)
}
plot_marginal_effects <- function(covs, preds) {
cbind(covs, preds) %>%
tidyr::gather(alpha, prediction, -seq_len(NCOL(covs))) %>%
dplyr::mutate(prediction = as.numeric(prediction)) %>%
tidyr::gather(variable, value, -(alpha:prediction)) %>%
ggplot(aes(value, prediction, color = alpha)) +
geom_point(alpha = 0.15) +
geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE) +
facet_wrap(~ variable, scale = "free_x") +
ylab(quote(q(y* "|" * x[1] * ",...," * x[p]))) +
xlab(quote(x[k])) +
theme(legend.position = "bottom")
}
plot_marginal_effects(abalone_f[, 1:3], pred_lqr)
We also want to check whether these results change, when we estimate the pair-copulas nonparametrically.
fit_vine_np <- vinereg(
whole ~ length + diameter + height,
data = abalone_f,
family_set = "nonpar",
selcrit = "aic"
)
fit_vine_np
## D-vine regression model: whole | length, height, diameter
## nobs = 1306, edf = 95.83, cll = 1207.07, caic = -2222.48, cbic = -1726.57
contour(fit_vine_np$vine)
Now only the length and height variables are selected as predictors. Let's have a look at the marginal effects.
plot_effects(fit_vine_np, var = c("diameter", "height", "length"))
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
The effects look similar to the parametric one, but slightly more wiggly. Note that even the diameter was not selected as a covariate, it's marginal effect is captured by the model. It just does not provide additional information when height and length are already accounted for.
To deal with discrete variables, we use the methodology of Schallhorn et al. (2017). For estmiating nonparametric pair-copulas with discrete variable(s), jittering is used (Nagler, 2017).
We let vinereg() know that a variable is discrete by declaring it ordered.
abalone_f$rings <- as.ordered(abalone_f$rings)
fit_disc <- vinereg(rings ~ ., data = abalone_f, selcrit = "aic")
fit_disc
## D-vine regression model: rings | shell, shucked, whole, viscera, height, diameter, length
## nobs = 1306, edf = 12, cll = -2763.64, caic = 5551.29, cbic = 5613.39
plot_effects(fit_disc)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Kraus and Czado (2017), D-vine copula based quantile regression, Computational Statistics and Data Analysis, 110, 1-18
Nagler (2017), A generic approach to nonparametric function estimation with mixed data, Statistics & Probability Letters, 137:326–330
Schallhorn, Kraus, Nagler and Czado (2017), D-vine quantile regression with discrete variables, arXiv preprint