Semiparametric Copula Bivariate Regression Models with Non-Random Sample Selection

Description

copulaSampleSel fit flexible copula bivariate sample selection models with several types of covariate effects, copula distributions and marginal distributions.

Usage

1
2
3
4
5
6
copulaSampleSel(formula, data = list(), weights = NULL, subset = NULL,  
                 BivD = "N", margins = c("probit", "N"), gamlssfit = FALSE,
                 fp = FALSE, infl.fac = 1, 
                 rinit = 1, rmax = 100, 
                 iterlimsp = 50, tolsp = 1e-07,
                 gc.l = FALSE, parscale, extra.regI = "t")

Arguments

formula

In the basic setup, this will be a list of two formulas, one for equation 1 and the other for equation 2. s terms are used to specify smooth smooth functions of predictors. For the case of more than two equations see the example below and the documentation of SemiParBIVProbit() for more details. Note that the first formula MUST refer to the selection equation.

data

An optional data frame, list or environment containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which copulaSampleSel is called.

weights

Optional vector of prior weights to be used in fitting.

subset

Optional vector specifying a subset of observations to be used in the fitting process.

margins

It indicates the distributions used for the two margins. The first is one of "probit", "logit", "cloglog" which refer to the link function of the first equation whose response is always assumed to be binary. The response for the second equation can be normal ("N"), log-normal ("LN"), Gumbel ("GU"), reverse Gumbel ("rGU"), logistic ("LO"), Weibull ("WEI"), inverse Gaussian ("iG"), gamma ("GA"), gamma with identity link for the location parameter ("GAi"), Dagum ("DAGUM"), Singh-Maddala ("SM"), beta ("BE"), Fisk ("FISK", also known as log-logistic distribution), Poisson ("PO"), zero truncated Poisson ("ZTP"), negative binomial - type I ("NBI"), negative binomial - type II ("NBII"), Poisson inverse Gaussian ("PIG").

gamlssfit

If gamlssfit = TRUE then a gamlss is also fitted for the outcome equation. This is may used for obtaining better calibrated starting values, for instance.

BivD

Type of bivariate error distribution employed. Possible choices are "N", "C0", "C90", "C180", "C270", "J0", "J90", "J180", "J270", "G0", "G90", "G180", "G270", "F", "AMH", "FGM" which stand for bivariate normal, Clayton, rotated Clayton (90 degrees), survival Clayton, rotated Clayton (270 degrees), Joe, rotated Joe (90 degrees), survival Joe, rotated Joe (270 degrees), Gumbel, rotated Gumbel (90 degrees), survival Gumbel, rotated Gumbel (270 degrees), Frank, Ali-Mikhail-Haq and Farlie-Gumbel-Morgenstern.

fp

If TRUE then a fully parametric model with unpenalised regression splines if fitted. See the example below.

infl.fac

Inflation factor for the model degrees of freedom in the approximate AIC. Smoother models can be obtained setting this parameter to a value greater than 1.

rinit

Starting trust region radius. The trust region radius is adjusted as the algorithm proceeds. See the documentation of trust for further details.

rmax

Maximum allowed trust region radius. This may be set very large. If set small, the algorithm traces a steepest descent path.

iterlimsp

A positive integer specifying the maximum number of loops to be performed before the smoothing parameter estimation step is terminated.

tolsp

Tolerance to use in judging convergence of the algorithm when automatic smoothing parameter estimation is used.

gc.l

This is relevant when working with big datasets. If TRUE then the garbage collector is called more often than it is usually done. This keeps the memory footprint down but it will slow down the routine.

parscale

The algorithm will operate as if optimizing objfun(x / parscale, ...) where parscale is a scalar. If missing then no rescaling is done. See the documentation of trust for more details.

extra.regI

If "t" then regularization as from trust is applied to the information matrix if needed. If different from "t" then extra regularization is applied via the options "pC" (pivoted Choleski - this will only work when the information matrix is semi-positive or positive definite) and "sED" (symmetric eigen-decomposition).

Details

The underlying algorithm is based on an extension of the procedure used for SemiParBIVProbit(). For more details see ?SemiParBIVProbit.

This function works as SemiParSampleSel() in SemiParSampleSel and has been included in SemiParBIVProbit (which already included sample selection models for binary outcomes) for the user's convenience (given several requests). copulaSampleSel() allows for several continuous/discrete distributions and link functions for the selection equation. SemiParSampleSel() allows for a probit link, normal, gamma and several discrete distributions.

If there are factors in the model, before fitting, the user has to ensure that the numbers of factor variables' levels in the selected sample are the same as those in the complete dataset. Even if a model could be fitted in such a situation, the model may produce fits which are not coherent with the nature of the correction sought. For more details see ?SemiParBIVProbit.

There are many continuous/discrete distributions and copula functions to choose from and we plan to include more options. Get in touch if you are interested in a particular distribution.

Value

The function returns an object of class copulaSampleSel as described in copulaSampleSelObject.

WARNINGS

Convergence failure may sometimes occur. Convergence can be checked using conv.check which provides some information about the score and information matrix associated with the fitted model. The former should be close to 0 and the latter positive definite. SemiParBIVProbit() will produce some warnings if there is a convergence issue.

In such a situation, the user may use some extra regularisation (see extra.regI) and/or rescaling (see parscale). Using gamlssfit = TRUE is typically more effective than the first two options as this may provide better calibrated starting values as compared to those obtained from the default starting value procedure. The default option is, however, gamlssfit = FALSE only because it tends to be computationally cheaper and because the default starting value procedure has typically been found to do a satisfactory job in most cases. (The results obtained when using gamlssfit = FALSE and gamlssfit = TRUE could also be compared to check if starting values make any difference.)

The above suggestions may help, especially the latter option. However, the user should also consider re-specifying the model and/or using a diferrent dependence structure and/or checking that the chosen marginal distributions are adequate. In our experience, we found that convergence failure typically occurs when the model has been misspecified and/or the sample size/number of selected observations is low compared to the complexity of the model. Examples of misspecification include using a Clayton copula rotated by 90 degrees when a positive association between the margins is present instead, using marginal distributions that do not fit the responses or that contain many parameters, and employing a copula which does not accommodate the type and/or strength of the dependence between the margins (e.g., using AMH when the association between the margins is strong). It is also worth bearing in mind that the use of a three parameter marginal distribution requires the data to be more informative than a situation in which a two parameter distribution is used instead.

Extra attention is required when specifying the dependence parameter as a function of covariates. This is because in these situations the dependence parameter mainly models the association between the unobserved confounders in the two equations. Therefore, this option would make sense when it is believed that the strength of the association between the unobservables in the two equations varies based on some grouping factor or across geographical areas, for instance.

Author(s)

Maintainer: Giampiero Marra giampiero.marra@ucl.ac.uk

References

Marra G. and Radice R. (2013), Estimation of a Regression Spline Sample Selection Model. Computational Statistics and Data Analysis, 61, 158-173.

Marra G. and Wyszynski K. (in press), Semi-Parametric Copula Sample Selection Models for Count Responses. Computational Statistics and Data Analysis.

Wojtys M. and Marra G. (submitted). Copula based generalized additive models with non-random sample selection.

See Also

copulaReg, SemiParBIVProbit, adjCov, VuongClarke, plot.SemiParBIVProbit, SemiParBIVProbit-package, copulaSampleSelObject, conv.check, summary.copulaSampleSel, predict.SemiParBIVProbit

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
library(SemiParBIVProbit)

######################################################################
## Generate data
## Correlation between the two equations and covariate correlation 0.5 
## Sample size 2000 
######################################################################

set.seed(0)

n <- 2000

rh <- 0.5      

sigmau <- matrix(c(1, rh, rh, 1), 2, 2)
u      <- rMVN(n, rep(0,2), sigmau)

sigmac <- matrix(rh, 3, 3); diag(sigmac) <- 1
cov    <- rMVN(n, rep(0,3), sigmac)
cov    <- pnorm(cov)

bi <- round(cov[,1]); x1 <- cov[,2]; x2 <- cov[,3]
  
f11 <- function(x) -0.7*(4*x + 2.5*x^2 + 0.7*sin(5*x) + cos(7.5*x))
f12 <- function(x) -0.4*( -0.3 - 1.6*x + sin(5*x))  
f21 <- function(x) 0.6*(exp(x) + sin(2.9*x)) 

ys <-  0.58 + 2.5*bi + f11(x1) + f12(x2) + u[, 1] > 0
y  <- -0.68 - 1.5*bi + f21(x1) +           u[, 2]
yo <- y*(ys > 0)
  
dataSim <- data.frame(ys, yo, bi, x1, x2)

## CLASSIC SAMPLE SELECTION MODEL
## the first equation MUST be the selection equation

resp.check(yo[ys > 0], "N")

out <- copulaSampleSel(list(ys ~ bi + x1 + x2, 
                            yo ~ bi + x1), 
                       data = dataSim)
conv.check(out)
post.check(out)
summary(out)

AIC(out)
BIC(out)


## SEMIPARAMETRIC SAMPLE SELECTION MODEL

## "cr" cubic regression spline basis      - "cs" shrinkage version of "cr"
## "tp" thin plate regression spline basis - "ts" shrinkage version of "tp"
## for smooths of one variable, "cr/cs" and "tp/ts" achieve similar results 
## k is the basis dimension - default is 10
## m is the order of the penalty for the specific term - default is 2

out <- copulaSampleSel(list(ys ~ bi + s(x1, bs = "tp", k = 10, m = 2) + s(x2), 
                            yo ~ bi + s(x1)), 
                       data = dataSim)
conv.check(out) 
post.check(out)
AIC(out)

## compare the two summary outputs
## the second output produces a summary of the results obtained when only 
## the outcome equation is fitted, i.e. selection bias is not accounted for

summary(out)
summary(out$gam2)

## estimated smooth function plots
## the red line is the true curve
## the blue line is the naive curve not accounting for selection bias

x1.s <- sort(x1[dataSim$ys>0])
f21.x1 <- f21(x1.s)[order(x1.s)] - mean(f21(x1.s))

plot(out, eq = 2, ylim = c(-1, 0.8)); lines(x1.s, f21.x1, col = "red")
par(new = TRUE)
plot(out$gam2, se = FALSE, lty = 3, lwd = 2, ylim = c(-1, 0.8), 
     ylab = "", rug = FALSE)


## Not run:  

## SEMIPARAMETRIC SAMPLE SELECTION MODEL with association 
## and dispersion parameters 
## depending on covariates as well

eq.mu.1   <- ys ~ bi + s(x1) + s(x2)
eq.mu.2   <- yo ~ bi + s(x1)
eq.sigma2 <-    ~ bi
eq.theta  <-    ~ bi + x1

fl <- list(eq.mu.1, eq.mu.2, eq.sigma2, eq.theta)

out <- copulaSampleSel(fl, data = dataSim)
conv.check(out)   
post.check(out)
summary(out)
out$sigma2
out$theta

jc.probs(out, 0, 0.3, intervals = TRUE)[1:4,]

outC0 <- copulaSampleSel(fl, data = dataSim, BivD = "C0")
conv.check(outC0)
post.check(outC0)
AIC(out, outC0)
BIC(out, outC0)

#
#

#######################################################
## example using Gumbel copula and normal-gamma margins
#######################################################

set.seed(1)

y  <- rgamma(n, shape = 1/2^2, exp(-0.68 - 1.5*bi + f21(x1) + u[, 2])*2^2)
yo <- y*(ys > 0)
  
dataSim <- data.frame(ys, yo, bi, x1, x2)


out <- copulaSampleSel(list(ys ~ bi + s(x1) + s(x2), 
                            yo ~ bi + s(x1)), 
                        data = dataSim, BivD = "G0", 
                        margins = c("probit", "GA"))
conv.check(out)
post.check(out)
summary(out)


#
#

## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.