biclaytoncop.RdEstimate the correlation parameter of the (bivariate) Clayton copula distribution by maximum likelihood estimation.
biclaytoncop(lapar = "loglink", iapar = NULL, imethod = 1,
parallel = FALSE, zero = NULL)Details at CommonVGAMffArguments.
See Links for more link function choices.
Details at CommonVGAMffArguments.
If parallel = TRUE then the constraint is also applied
to the intercept.
The cumulative distribution function is $$P(u_1, u_2;\alpha) = (u_1^{-\alpha} + u_2^{-\alpha}-1)^{-1/\alpha}$$ for \(0 \leq \alpha \). Here, \(\alpha\) is the association parameter. The support of the function is the interior of the unit square; however, values of 0 and/or 1 are not allowed (currently). The marginal distributions are the standard uniform distributions. When \(\alpha = 0\) the random variables are independent.
This VGAM family function can handle multiple responses, for example, a six-column matrix where the first 2 columns is the first out of three responses, the next 2 columns being the next response, etc.
An object of class "vglmff"
(see vglmff-class).
The object is used by modelling functions
such as vglm
and vgam.
Clayton, D. (1982). A model for association in bivariate survival data. Journal of the Royal Statistical Society, Series B, Methodological, 44, 414–422.
Schepsmeier, U. and Stober, J. (2014). Derivatives and Fisher information of bivariate copulas. Statistical Papers 55, 525–542.
The response matrix must have a multiple of two-columns. Currently, the fitted value is a matrix with the same number of columns and values equal to 0.5. This is because each marginal distribution corresponds to a standard uniform distribution.
This VGAM family function is fragile; each response must be in the interior of the unit square.
ymat <- rbiclaytoncop(n = (nn <- 1000), apar = exp(2))
bdata <- data.frame(y1 = ymat[, 1], y2 = ymat[, 2],
y3 = ymat[, 1], y4 = ymat[, 2], x2 = runif(nn))
summary(bdata)
#> y1 y2 y3 y4
#> Min. :0.000828 Min. :0.0008981 Min. :0.000828 Min. :0.0008981
#> 1st Qu.:0.254812 1st Qu.:0.2605758 1st Qu.:0.254812 1st Qu.:0.2605758
#> Median :0.508578 Median :0.5055172 Median :0.508578 Median :0.5055172
#> Mean :0.507874 Mean :0.5064419 Mean :0.507874 Mean :0.5064419
#> 3rd Qu.:0.764826 3rd Qu.:0.7455065 3rd Qu.:0.764826 3rd Qu.:0.7455065
#> Max. :0.999361 Max. :0.9995598 Max. :0.999361 Max. :0.9995598
#> x2
#> Min. :0.002145
#> 1st Qu.:0.270375
#> Median :0.526698
#> Mean :0.519788
#> 3rd Qu.:0.780546
#> Max. :0.999487
if (FALSE) plot(ymat, col = "blue") # \dontrun{}
fit1 <-
vglm(cbind(y1, y2, y3, y4) ~ 1, # 2 responses, e.g., (y1,y2) is the 1st
biclaytoncop, data = bdata,
trace = TRUE, crit = "coef") # Sometimes a good idea
#> Iteration 1: coefficients = 1.993075, 1.991763
#> Iteration 2: coefficients = 1.9935812, 1.9935912
#> Iteration 3: coefficients = 1.9935769, 1.9935769
#> Iteration 4: coefficients = 1.993577, 1.993577
coef(fit1, matrix = TRUE)
#> loglink(apar1) loglink(apar2)
#> (Intercept) 1.993577 1.993577
Coef(fit1)
#> apar1 apar2
#> 7.341748 7.341748
head(fitted(fit1))
#> y1 y2 y3 y4
#> 1 0.5 0.5 0.5 0.5
#> 2 0.5 0.5 0.5 0.5
#> 3 0.5 0.5 0.5 0.5
#> 4 0.5 0.5 0.5 0.5
#> 5 0.5 0.5 0.5 0.5
#> 6 0.5 0.5 0.5 0.5
summary(fit1)
#>
#> Call:
#> vglm(formula = cbind(y1, y2, y3, y4) ~ 1, family = biclaytoncop,
#> data = bdata, trace = TRUE, crit = "coef")
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 1.99358 0.03099 64.34 <2e-16 ***
#> (Intercept):2 1.99358 0.03099 64.34 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Names of linear predictors: loglink(apar1), loglink(apar2)
#>
#> Log-likelihood: 2439.828 on 1998 degrees of freedom
#>
#> Number of Fisher scoring iterations: 4
#>
#> No Hauck-Donner effect found in any of the estimates
#>
# Another example; apar is a function of x2
bdata <- transform(bdata, apar = exp(-0.5 + x2))
ymat <- rbiclaytoncop(n = nn, apar = with(bdata, apar))
bdata <- transform(bdata, y5 = ymat[, 1], y6 = ymat[, 2])
fit2 <- vgam(cbind(y5, y6) ~ s(x2), data = bdata,
biclaytoncop(lapar = "loglink"), trace = TRUE)
#> VGAM s.vam loop 1 : loglikelihood = 186.30248
#> VGAM s.vam loop 2 : loglikelihood = 186.37729
#> VGAM s.vam loop 3 : loglikelihood = 186.37027
#> VGAM s.vam loop 4 : loglikelihood = 186.36697
#> VGAM s.vam loop 5 : loglikelihood = 186.36587
#> VGAM s.vam loop 6 : loglikelihood = 186.36553
#> VGAM s.vam loop 7 : loglikelihood = 186.36542
#> VGAM s.vam loop 8 : loglikelihood = 186.36539
#> VGAM s.vam loop 9 : loglikelihood = 186.36537
#> VGAM s.vam loop 10 : loglikelihood = 186.36537
if (FALSE) plot(fit2, lcol = "blue", scol = "orange", se = TRUE) # \dontrun{}