SemiParBIVProbit: Semiparametric Bivariate Probit Modelling

Description

SemiParBIVProbit can be used to fit bivariate probit models where the linear predictors can be flexibly specified using parametric and regression spline components. The use of several bivariate copula distributions and of asymmetric link functions is supported. During the model fitting process, the possible presence of correlated error equations, endogeneity or non-random sample selection is accounted for. Regression spline bases are extracted from the package mgcv. Multi-dimensional smooths are available via the use of penalized thin plate regression splines (isotropic). The current implementation does not support scale invariant tensor product smooths.

Usage

SemiParBIVProbit(formula.eq1, formula.eq2, data = list(), weights = NULL,  
                 start.v = NULL, BivD = "N", nu = 3, PL = "P", eqPL = "both", 
                 valPL=c(0,0), fitPL="pLiksp", spPL=c(0.01,0.01), 
                 selection = FALSE, H.n = TRUE, gamma = 1, aut.sp = TRUE, 
                 fp = FALSE, pPen1 = NULL, pPen2 = NULL, rinit = 1, rmax = 100, 
                 fterm = sqrt(.Machine$double.eps), 
                 mterm = sqrt(.Machine$double.eps), 
                 iterlimsp = 50, pr.tolsp = 1e-6,    
                 control.sp = list(maxit = 50, tol = 1e-6, step.half = 25,
                                   rank.tol = sqrt(.Machine$double.eps)) )

Arguments

formula.eq1

A formula for equation 1. s terms are used to specify smooth smooth functions of predictors. SemiParBIVProbit supports the use shrinkage smoothers for variable selection purposes. See the ex

formula.eq2

A formula for equation 2.

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 SemiParBIVProbit

weights

Optional vector of prior weights to be used in fitting.

start.v

Starting values for the parameters of the two equations and association coefficient can be provided here.

BivD

Type of bivariate error distribution employed. Possible choices are "N", "C0", "C90", "C180", "C270", "J0", "J90", "J180", "J270", "G0", "G90", "G180", "G270", "F", "T" which stand for bivariate normal, Clayton, rotated Clayton (90 degrees),

Degrees of freedom for Student-t.

Link function to be employed. Possible choices are "P" (classic probit), "PP" (power probit), "RPP" (reciprocal power probit) and "SN" (skew normal).

eqPL

Equation(s) on which the asymmetric link function approach should be applied. Possible choices are "both", "first" and "second".

valPL

Initial values for link function shape parameters when an asymmetric link fucntion approach is employed.

fitPL

Fitting approach for shape parameters of the asymmetric link functions employed. Possible choices are "fixed", "unpLik", "pLik" and "pLiksp" which stand for fixed shape parameters, classic likelihood optimization, penalized likelihood with

spPL

Values for the smoothing parameters of the ridge penalties associated with the shape parameteres of the asymmetric link functions.

selection

If TRUE, then the numerical routine for bivariate probit modelling in the presence of non-random sample selection is employed.

H.n

For BivD = "N" only, if FALSE then the expected (rather than observed) information matrix is employed.

gamma

Inflation factor for the model degrees of freedom in the UBRE score. Smoother models can be obtained setting this parameter to a value greater than 1. Typically gamma = 1.4 achieves this.

aut.sp

If TRUE, then automatic multiple smoothing parameter selection is carried out. If FALSE, then smoothing parameters are set to the values obtained from the univariate fits.

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

pPen1,pPen2

Optional list specifying any penalties to be applied to the parametric model terms of equations 1 and 2.

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.

fterm

Positive scalar giving the tolerance at which the difference in objective function values in a step is considered close enough to zero to terminate the algorithm.

mterm

Positive scalar giving the tolerance at which the two-term Taylor-series approximation to the difference in objective function values in a step is considered close enough to zero to terminate the algorithm.

iterlimsp

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

pr.tolsp

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

control.sp

List containing iteration control constants for the automatic smoothing parameter selection method. These are: maxit: maximum number of iterations of the magic algorithm; tol: toler

Value

The function returns an object of class SemiParBIVProbit as described in SemiParBIVProbitObject.

WARNINGS

In the current context, convergence failure may occur when the number of observations is low and the data have low information content. When convergence failure is associated with an infinite cycling between the two steps detailed above, several practical solutions are available. For instance, one might either (i) lower the total number of parameters to estimate by reducing the dimension of the regression spline bases, (ii) set the smoothing parameters to the values obtained from the univariate fits (aut.sp=FALSE), or (iii) set the smoothing parameters to the values obtained from the non-converged algorithm. The default option is (iii).

Details

The bivariate model has a probit link for each of the two equations, and models the association between the responses through the correlation parameter $\rho$ of a standardised bivariate normal distribution, or that of a bivariate copula distribution, $\theta$. In a semiparametric bivariate probit model the linear predictors are flexibly specified using parametric components and smooth functions of covariates. Replacing the smooth components with their regression spline expressions yields a fully parametric bivariate probit model. In principle, classic maximum likelihood estimation can be employed. However, to avoid overfitting, penalized likelihood maximization has to be employed instead. Here the use of penalty matrices allows for the suppression of that part of smooth term complexity which has no support from the data. The trade-off between smoothness and fitness is controlled by smoothing parameters associated with the penalty matrices. Smoothing parameters are chosen to minimize the approximate Un-Biased Risk Estimator (UBRE). Automatic smoothing parameter selection is integrated using a performance-oriented iteration approach (Gu, 1992; Wood, 2004). At each iteration, (i) the penalized weighted least squares problem is solved and (ii) the smoothing parameters of that problem estimated by approximate UBRE. Steps (i) and (ii) are iterated until convergence. Details of the underlying fitting methods are given in Radice, Marra and Wojtys (submitted). Releases previous to 3.2-7 were based on the algorithms detailed in Marra and Radice (2011, 2013). Asymmetric link functions are available; see Marra and Radice (submitted) for details.

References

Gu C. (1992), Cross validating non-Gaussian data. Journal of Computational and Graphical Statistics, 1(2), 169-179. Marra G. and Radice R. (2011), Estimation of a Semiparametric Recursive Bivariate Probit in the Presence of Endogeneity. Canadian Journal of Statistics, 39(2), 259-279. Marra G. and Radice R. (2013), A Penalized Likelihood Estimation Approach to Semiparametric Sample Selection Binary Response Modeling. Electronic Journal of Statistics, 7, 1432-1455. Marra G. and Radice R. (submitted), On the Use of Asymmetric Link Functions in Copula Regression Spline Models for Bivariate Binary Responses. Marra G., Radice R. and Missiroli S. (2014), Testing the Hypothesis of Absence of Unobserved Confounding in Semiparametric Bivariate Probit Models. Computational Statistics, 29(3-4), 715-741. Marra G., Radice R. and Filippou P. (submitted), Testing the Hypothesis of Exogeneity in Regression Spline Bivariate Probit Models. Radice R., Marra G. and M. Wojtys (submitted), Copula Regression Spline Models for Binary Outcomes. Wood S.N. (2004), Stable and efficient multiple smoothing parameter estimation for generalized additive models. Journal of the American Statistical Association, 99(467), 673-686.

Examples

Run this code

library(SemiParBIVProbit)

############
## EXAMPLE 1
############
## Generate data
## Correlation between the two equations 0.5 - Sample size 400 

set.seed(0)

n <- 400

Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u     <- rmvnorm(n, rep(0,2), Sigma)

x1 <- round(runif(n)); x2 <- runif(n); x3 <- runif(n)

f1   <- function(x) cos(pi*2*x) + sin(pi*x)
f2   <- function(x) x+exp(-30*(x-0.5)^2)   

y1 <- ifelse(-1.55 + 2*x1    + f1(x2) + u[,1] > 0, 1, 0)
y2 <- ifelse(-0.25 - 1.25*x1 + f2(x2) + u[,2] > 0, 1, 0)

dataSim <- data.frame(y1, y2, x1, x2, x3)

#
#

## CLASSIC BIVARIATE PROBIT

out  <- SemiParBIVProbit(y1 ~ x1 + x2 + x3, 
                         y2 ~ x1 + x2 + x3, 
                         data = dataSim)
summary(out)
AIC(out)
BIC(out)

## SEMIPARAMETRIC BIVARIATE PROBIT

## "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  <- SemiParBIVProbit(y1 ~ x1 + s(x2, bs = "tp", k = 10, m = 2) + s(x3), 
                         y2 ~ x1 + s(x2) + s(x3), 
                         data = dataSim)
sem.checks(out)
summary(out)
AIC(out)


## estimated smooth function plots - red lines are true curves

x2 <- sort(x2)
f1.x2 <- f1(x2)[order(x2)] - mean(f1(x2))
f2.x2 <- f2(x2)[order(x2)] - mean(f2(x2))
f3.x3 <- rep(0, length(x3))

par(mfrow=c(2,2),mar=c(4.5,4.5,2,2))
plot(out, eq = 1, select = 1, scale = 0); lines(x2, f1.x2, col = "red")
plot(out, eq = 1, select = 2, scale = 0); lines(x3, f3.x3, col = "red")
plot(out, eq = 2, select = 1, scale = 0); lines(x2, f2.x2, col = "red")
plot(out, eq = 2, select = 2, scale = 0); lines(x3, f3.x3, col = "red")

#
## same plots but CIs 'with intercept' 

par(mfrow=c(2,2),mar=c(4.5,4.5,2,2))
plot(out, eq = 1, select = 1, seWithMean = TRUE, scale = 0); lines(x2, f1.x2, col = "red")
plot(out, eq = 1, select = 2, seWithMean = TRUE, scale = 0); lines(x3, f3.x3, col = "red")
plot(out, eq = 2, select = 1, seWithMean = TRUE, scale = 0); lines(x2, f2.x2, col = "red")
plot(out, eq = 2, select = 2, seWithMean = TRUE, scale = 0); lines(x3, f3.x3, col = "red")


## p-values suggest to drop x3 from both equations, with a stronger 
## evidence for eq. 2. This can be also achieved using shrinkage smoothers

outSS <- SemiParBIVProbit(y1 ~ x1 + s(x2, bs = "cs") + s(x3, bs = "cs"), 
                          y2 ~ x1 + s(x2, bs = "cs") + s(x3, bs = "cs"), 
                          data = dataSim)

par(mfrow = c(2,2), mar = c(4.5,4.5,2,2))
plot(outSS, eq = 1, select = 1, scale = 0)
plot(outSS, eq = 1, select = 2, ylim = c(-0.1,0.1))
plot(outSS, eq = 2, select = 1, scale = 0)
plot(outSS, eq = 2, select = 2, ylim = c(-0.1,0.1))

#
#

############
## EXAMPLE 2
############
## Generate data with one endogenous variable and exclusion restriction

set.seed(0)

n <- 300

Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u     <- rmvnorm(n, rep(0,2), Sigma)

cov   <- rmvnorm(n, rep(0,2), Sigma, method = "svd")
cov   <- pnorm(cov)
x1 <- round(cov[,1]); x2 <- cov[,2]

f1   <- function(x) cos(pi*2*x) + sin(pi*x)
f2   <- function(x) x+exp(-30*(x-0.5)^2)   

y1 <- ifelse(-1.55 + 2*x1    + f1(x2) + u[,1] > 0, 1, 0)
y2 <- ifelse(-0.25 - 1.25*y1 + f2(x2) + u[,2] > 0, 1, 0)

dataSim <- data.frame(y1, y2, x1, x2)

#

## Testing the hypothesis of absence of endogeneity... 

LM.bpm(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2), dataSim)

# p-value suggests presence of endogeneity, hence fit a bivariate model


## CLASSIC RECURSIVE BIVARIATE PROBIT

out <- SemiParBIVProbit(y1 ~ x1 + x2, 
                        y2 ~ y1 + x2, 
                        data = dataSim)
summary(out)

## SEMIPARAMETRIC RECURSIVE BIVARIATE PROBIT

out <- SemiParBIVProbit(y1 ~ x1 + s(x2), 
                        y2 ~ y1 + s(x2), 
                        data = dataSim)
summary(out)

#

## Testing the hypothesis of absence of endogeneity post estimation... 

gt.bpm(out)

#
## average treatment effect with CIs

AT(out, eq = 2, nm.bin = "y1") 


## try a Clayton copula model... 

outC <- SemiParBIVProbit(y1 ~ x1 + s(x2), 
                         y2 ~ y1 + s(x2), 
                         data = dataSim, BivD = "C0")
summary(outC)
AT(outC, eq = 2, nm.bin = "y1") 

## try a Joe copula model... 

outJ <- SemiParBIVProbit(y1 ~ x1 + s(x2), 
                         y2 ~ y1 + s(x2), 
                         data = dataSim, BivD = "J0")
summary(outJ)
AT(outJ, eq = 2, nm.bin = "y1") 


VuongClarke.bcm(out, outJ)

#
## recursive bivariate probit modelling with unpenalized splines 
## can be achieved as follows

outFP <- SemiParBIVProbit(y1 ~ x1 + s(x2, bs = "cr", k = 5), 
                          y2 ~ y1 + s(x2, bs = "cr", k = 6), fp = TRUE, 
                          data = dataSim)
summary(outFP)

#
#

############
## EXAMPLE 3
############
## Generate data with a non-random sample selection mechanism and exclusion restriction

set.seed(0)

n <- 1100

Sigma <- matrix(0.5, 2, 2); diag(Sigma) <- 1
u     <- rmvnorm(n, rep(0,2), Sigma)

SigmaC <- matrix(0.5, 3, 3); diag(SigmaC) <- 1
cov    <- rmvnorm(n, rep(0,3), SigmaC, method = "svd")
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] > 0
yo <- y*(ys > 0)
  
dataSim <- data.frame(y, ys, yo, bi, x1, x2)

## Testing the hypothesis of absence of non-random sample selection... 

LM.bpm(ys ~ bi + s(x1) + s(x2), yo ~ bi + s(x1), dataSim, selection = TRUE)

# p-value suggests presence of sample selection, hence fit a bivariate model

#
## SEMIPARAMETRIC SAMPLE SELECTION BIVARIATE PROBIT
## the first equation MUST be the selection equation

out <- SemiParBIVProbit(ys ~ bi + s(x1) + s(x2), 
                        yo ~ bi + s(x1), 
                        data = dataSim, selection = TRUE)
gt.bpm(out)                        

## compare the two summary outputs
## the second output produces a summary of the results obtained when
## selection bias is not accounted for

summary(out)
summary(out$gam2)

## corrected predicted probability that 'yo' is equal to 1

est.prev(out)

## 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))

x11()
plot(out, eq = 2, ylim = c(-1.65,0.95)); lines(x1.s, f21.x1, col="red")
par(new = TRUE)
plot(out$gam2, se = FALSE, col = "blue", ylim = c(-1.65,0.95), ylab = "", rug = FALSE)

#
#
## try a Clayton copula model... 

outC <- SemiParBIVProbit(ys ~ bi + s(x1) + s(x2), 
                         yo ~ bi + s(x1), 
                         data = dataSim, selection = TRUE, BivD = "C0")
sem.checks(outC)
summary(outC)
est.prev(outC)

#
#

Run the code above in your browser using DataLab