`copulaSampleSel`

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

1 2 3 4 5 |

`formula` |
In the basic setup, this will be a list of two formulas, one for equation 1 and the other for equation 2. |

`data` |
An optional data frame, list or environment containing the variables in the model. If not found in |

`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"), normal where sigma2 corresponds to the standard deviation instead of the variance ("N2"), 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"). |

`dof` |
If |

`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", "T" 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,
Farlie-Gumbel-Morgenstern and Student-t with fixed |

`fp` |
If |

`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 |

`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 |

`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 |

`extra.regI` |
If "t" then regularization as from |

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.

The function returns an object of class `copulaSampleSel`

as described in `copulaSampleSelObject`

.

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.

Convergence failure may sometimes occur. This is not necessarily a bad thing as it may indicate specific problems
with a fitted model. In such a situation, the user may use some extra regularisation (see `extra.regI`

) and/or
rescaling (see `parscale`

). However, the user should also consider
re-specifying/simplifying 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). When using
smooth functions, if the covariate's values are too sparse then convergence may be affected by this.
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.

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

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.

`copulaReg`

, `SemiParBIVProbit`

, `adjCov`

, `VuongClarke`

, `plot.SemiParBIVProbit`

, `SemiParBIVProbit-package`

, `copulaSampleSelObject`

, `conv.check`

, `summary.copulaSampleSel`

, `predict.SemiParBIVProbit`

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 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 | ```
## Not run:
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)
## IMPUTE MISSING VALUES
n.m <- 10
res <- imputeSS(out, n.m)
bet <- NA
for(i in 1:n.m){
dataSim$yo[dataSim$ys == 0] <- res[[i]]
outg <- gamlss(list(yo ~ bi + s(x1)), data = dataSim)
bet[i] <- coef(outg)["bi"]
print(i)
}
mean(bet)
##
## 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)
## IMPUTE MISSING VALUES
n.m <- 10
res <- imputeSS(outC0, n.m)
#
#
#######################################################
## 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)
n.m <- 10
res <- imputeSS(out, n.m)
for(i in 1:n.m){
dataSim$yo[dataSim$ys == 0] <- res[[i]]
outg <- gamlss(list(yo ~ bi + s(x1)), margin = "GA", data = dataSim)
bet[i] <- coef(outg)["bi"]
print(i)
}
mean(bet)
#
#
## End(Not run)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.