Package ‘SemiParBIVProbit’ January 15, 2015 Version 3.2-13.2 Author Giampiero Marra <giampiero.marra@ucl.ac.uk> and Rosalba Radice <r.radice@bbk.ac.uk> Maintainer Giampiero Marra <giampiero.marra@ucl.ac.uk> Title Semiparametric Bivariate Probit Modelling Description Routine for fitting bivariate models for binary responses with semiparametric predictors (including linear and nonlinear effects) in the presence of correlated error equations, endogeneity, non-random sample selection or partial observability. Depends R (>= 3.1.1), mgcv, mvtnorm Imports magic, CDVine, VineCopula, VGAM, survey, trust, Matrix, sn LazyLoad yes License GPL (>= 2) URL http://www.ucl.ac.uk/statistics/people/giampieromarra Repository CRAN Date/Publication 2015-01-15 12:45:10 NeedsCompilation no R topics documented: SemiParBIVProbit-package . adjCov . . . . . . . . . . . . adjCovSD . . . . . . . . . . AT . . . . . . . . . . . . . . bprobgHs . . . . . . . . . . bprobgHsPL . . . . . . . . . bprobgHsPO . . . . . . . . bprobgHsSS . . . . . . . . . copgHs . . . . . . . . . . . est.prev . . . . . . . . . . . gt.bpm . . . . . . . . . . . . LM.bpm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 4 5 6 8 8 9 9 10 10 11 12 2 SemiParBIVProbit-package logLik.SemiParBIVProbit . . . . mb . . . . . . . . . . . . . . . . . pen . . . . . . . . . . . . . . . . . penPL . . . . . . . . . . . . . . . plot.SemiParBIVProbit . . . . . . predict.SemiParBIVProbit . . . . print.AT . . . . . . . . . . . . . . print.est.prev . . . . . . . . . . . print.mb . . . . . . . . . . . . . . print.SemiParBIVProbit . . . . . . print.summary.SemiParBIVProbit S.m . . . . . . . . . . . . . . . . S.mPL . . . . . . . . . . . . . . . sem.checks . . . . . . . . . . . . SemiParBIVProbit . . . . . . . . SemiParBIVProbit.fit . . . . . . . SemiParBIVProbit.fit.post . . . . SemiParBIVProbit.fit1 . . . . . . SemiParBIVProbitObject . . . . . startSS . . . . . . . . . . . . . . . summary.SemiParBIVProbit . . . VuongClarke.bcm . . . . . . . . . working.comp . . . . . . . . . . . working.comp1 . . . . . . . . . . working.comp3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 15 16 16 17 18 19 20 20 21 22 23 23 23 24 34 35 35 35 37 38 41 42 43 43 44 SemiParBIVProbit-package Semiparametric Bivariate Probit Modelling Description SemiParBIVProbit provides a function for fitting bivariate probit models with semiparametric predictors, including linear and nonlinear effects. Several bivariate copula distributions are supported. Smoothness selection is achieved automatically and interval calculations are based on a Bayesian approach. Details SemiParBIVProbit provides a function for fitting flexible bivariate probit models, in the presence of correlated error equations, endogeneity, non-random sample selection or partial observability. The underlying representation and estimation of the model is based on a penalized likelihood-based regression spline approach, with automatic smoothness selection. Several bivariate copula distributions are available. The numerical routine carries out function minimization using a trust region algorithm from the package trust in combination with an adaptation of a low level smoothness estimation fitting procedure from the package mgcv. SemiParBIVProbit-package 3 SemiParBIVProbit supports the use of many smoothers as extracted from mgcv. Estimation is by penalized maximum likelihood with automatic smoothness estimation achieved by using the approximate Un-Biased Risk Estimator (UBRE), which can also be viewed as an approximate AIC. Confidence intervals for smooth components and nonlinear functions of the model parameters are derived using a Bayesian approach. Approximate p-values for testing individual smooth terms for equality to the zero function are also provided. Functions plot.SemiParBIVProbit and summary.SemiParBIVProbit extract such information from a fitted SemiParBIVProbit object. Model/variable selection is also possible via the use of shrinakge smoothers or information criteria. Tools for testing the hypothesis of uncorrelated error equations/absence of unobserved confounding/absence of endogeneity/absence of non-random sample selection are available (see gt.bpm and LM.bpm). For recursive bivariate and sample selection models AT and est.prev calculate the average effect of an endogenous covariate and corrected prevalence. mb provides the nonparametric (worst-case) Manski’s bound which is useful to check whether the average effect from a recursive model is included within the possibilites of the bound. Models with asymmetric link functions are also implemented. However, in our experience, this approach is not likely to lead to significant improvements as compared to using probit links. Author(s) Giampiero Marra (University College London, Department of Statistical Science) and Rosalba Radice (Birkbeck, University of London, Department of Economics, Mathematics and Statistics) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> Thanks to Bear Braumoeller (Department of Political Science, The Ohio State University) for suggesting the implementation of bivariate models with partial observability. References Marra G. (2013), On P-values for Semiparametric Bivariate Probit Models. Statistical Methodology, 10(1), 23-28. 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. (2015), Flexible Bivariate Binary Models for Estimating the Efficacy of Phototherapy for Newborns with Jaundice. International Journal of Statistics and Probability. 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), 715741. Marra G., Radice R. and Filippou P. (submitted), Regression Spline Bivariate Probit Models: A Practical Approach to Testing for Exogeneity. McGovern M.E., Barnighausen T., Marra G. and Radice R. (2015), On the Assumption of Joint Normality in Selection Models: A Copula Based Approach Applied to Estimating HIV Prevalence. Epidemiology. Radice R., Marra G. and M. Wojtys (submitted), Copula Regression Spline Models for Binary Outcomes. 4 adjCov See Also SemiParBIVProbit adjCov Adjustment for the covariance matrix from a fitted SemiParBIVProbit model Description adjCov can be used to adjust the covariance matrix of a fitted SemiParBIVProbit object. Usage adjCov(object, id) Arguments object A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). id Cluster identifier. Details This adjustment can be made when dealing with clustered data and the cluster structure is neglected when fitting the model. The basic idea is that the model is fitted as though observations were independent, and subsequently adjust the covariance matrix of the parameter estimates. Using the terminology of Liang and Zeger (1986), this would correspond to using an independence structure within the context of generalized estimating equations. The parameter estimators are still consistent but are inefficient as compared to a model which accounts for the correct cluster dependence structure. The covariance matrix of the independence estimators can be adjusted as described Liang and Zeger (1986, Section 2). Value This function returns a fitted SemiParBIVProbit object which is identical to that supplied in adjCov but with adjusted covariance matrix. WARNINGS This correction may not be appropriate for models that require the use of penalties in fitting. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> adjCovSD 5 References Liang K.-Y. and Zeger S. (1986), Longitudinal Data Analysis Using Generalized Linear Models. Biometrika, 73(1), 13-22. See Also SemiParBIVProbit-package, SemiParBIVProbit, summary.SemiParBIVProbit adjCovSD Adjustment for the covariance matrix from a SemiParBIVProbit model fitted to complex survey data. Description adjCovSD can be used to adjust the covariance matrix of a fitted SemiParBIVProbit object. Usage adjCovSD(object, design) Arguments object A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). design A svydesign object as produced by svydesign() from the survey package. Details This function has been extracted from the survey package and adapted to the class of bivariate models. It computes the sandwich variance estimator for a bivariate model fitted to data from a complex sample survey (Lumley, 2004). Value This function returns a fitted SemiParBIVProbit object which is identical to that supplied in adjCovSD but with adjusted covariance matrix. WARNINGS This correction may not be appropriate for models that require the use of penalties in fitting. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> 6 AT References Lumley T. (2004), Analysis of Complex Survey Samples. Journal of Statistical Software, 9(8), 1-19. See Also SemiParBIVProbit-package, SemiParBIVProbit, summary.SemiParBIVProbit AT Average treatment effect of a binary endogenous variable Description AT can be used to calculate the sample average treatment effect of a binary endogenous predictor/treatment, with corresponding confidence/credible intervals. Usage AT(x, eq, nm.bin = "", E = TRUE, treat = TRUE, delta = FALSE, prob.lev = 0.05, s.meth = "svd", n.sim = 500, naive = FALSE) Arguments x A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). eq Equation containing the binary endogenous predictor of interest. nm.bin Name of the binary endogenous variable. E If TRUE then AT calculates the sample ATE. If FALSE then it calculates the sample AT for the treated individuals only. treat If TRUE then AT calculates the AT using the treated only. If FALSE then it calculates the effect on the control group. This only makes sense if used jointly with E = FALSE. delta If TRUE then an approximate delta method is used for confidence interval calculation, otherwise Bayesian posterior simulation (the most reliable option, despite a bit slower) is employed. Note that for models involving asymmetric link functions only posterior simulation is allowed for. prob.lev Probability of the left and right tails of the AT distribution used for interval calculations. s.meth Matrix decomposition used to determine the matrix root of the covariance matrix. This is used when delta = FALSE. See the documentation of the mvtnorm package for further details. n.sim Number of simulated coefficient vectors from the posterior distribution of the estimated model parameters. This is used when delta = FALSE. It may be increased if more precision is required. AT 7 naive It indicates whether the AT is calculated from a univariate probit model. This option has been introduced to compare adjusted (for unobserved confounding) and unadjusted estimates. Note that, although the unadjusted estimates do not account for unobserved confounding, they do control for observed confounding. Details AT measures the sample average difference in outcomes under treatment (the binary predictor or treatment assumes value 1) and under control (the binary treatment assumes value 0). Note that this function is only suitable for the calculation of the average treatment effect in recursive bivariate binary models. Posterior simulation and delta method can be used to obtain confidence/credible intervals. The former is typically more reliable at small sample sizes. See the references below for details. Value res It returns three values: lower confidence interval limit, estimated AT and upper confidence interval limit. prob.lev Probability level used. est.ATb If delta = FALSE then it returns a vector containing simulated values of the average treatment effect. This is used to calculate intervals. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References 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. Radice R., Marra G. and M. Wojtys (submitted), Copula Regression Spline Models for Binary Outcomes. See Also SemiParBIVProbit-package, SemiParBIVProbit, summary.SemiParBIVProbit Examples ## see examples for SemiParBIVProbit 8 bprobgHsPL bprobgHs Internal Function Description It provides the log-likelihood, gradient and observed or expected information matrix for penalized or unpenalized maximum likelihood optimization. Possible bivariate distributions are 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 and Student-t with fixed degrees of freedom. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References 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. Radice R., Marra G. and M. Wojtys (submitted), Copula Regression Spline Models for Binary Outcomes. bprobgHsPL Internal Function Description It provides the log-likelihood, gradient and observed or expected information matrix for penalized or unpenalized maximum likelihood optimization, when using asymmetric link functions. Possible bivariate distributions are 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 and Student-t with fixed degrees of freedom. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References Marra G. and Radice R. (2015), Flexible Bivariate Binary Models for Estimating the Efficacy of Phototherapy for Newborns with Jaundice. International Journal of Statistics and Probability. bprobgHsPO bprobgHsPO 9 Internal Function Description It provides the log-likelihood, gradient and observed or expected information matrix for penalized or unpenalized maximum likelihood optimization for the case of bivariate models with partial observability. Possible bivariate distributions are 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 and Student-t with fixed degrees of freedom. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> bprobgHsSS Internal Function Description It provides the log-likelihood, gradient and observed or expected information matrix for penalized or unpenalized maximum likelihood optimization, for the non-random sample selection model case. Possible bivariate distributions are 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 and Student-t with fixed degrees of freedom. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References 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. McGovern M.E., Barnighausen T., Marra G. and Radice R. (2015), On the Assumption of Joint Normality in Selection Models: A Copula Based Approach Applied to Estimating HIV Prevalence. Epidemiology. 10 est.prev copgHs Internal Function Description This function evaluates the first and second derivatives with respect to the margins and association parameter of several bivariate copulas. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> est.prev Estimated overall prevalence from sample selection model Description est.prev can be used to calculate the overall estimated prevalence from a sample selection model, with corresponding confidence intervals obtained using the delta method. Usage est.prev(x, sig.lev = 0.05, sw = NULL, naive = FALSE) Arguments x A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). sig.lev Significance level. sw Survey weights. naive If FALSE then the prevalence is calculated using the (naive/classic imputation) probit model. This option has been introduced to compare adjusted (for nonrandom sample selection) and unadjusted estimates. Details est.prev estimates the overall prevalence of a disease (e.g., HIV) when there are missing values that are not at random. Value res It returns three values: lower confidence interval limit, estimated prevalence and upper confidence interval limit. sig.lev Significance level used. gt.bpm 11 Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References McGovern M.E., Barnighausen T., Marra G. and Radice R. (2015), On the Assumption of Joint Normality in Selection Models: A Copula Based Approach Applied to Estimating HIV Prevalence. Epidemiology. See Also SemiParBIVProbit-package, SemiParBIVProbit, summary.SemiParBIVProbit Examples ## see examples for SemiParBIVProbit gt.bpm Gradient test Description gt.bpm can be used to test the hypothesis of absence of endogeneity, correlated model equations/errors or non-random sample selection. Usage gt.bpm(object) Arguments object A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). Details The gradient test was first proposed by Terrell (2002) and it is based on classic likelihood theory. See Marra et al. (submitted) for full details. Value It returns a numeric p-value corresponding to the null hypothesis that ρ is equal to 0. WARNINGS This test’s implementation is only valid for models with normal errors and probit links. 12 LM.bpm Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References Marra G., Radice R. and Filippou P. (submitted), Regression Spline Bivariate Probit Models: A Practical Approach to Testing for Exogeneity. Terrell G. (2002), The Gradient Statistic. Computing Science and Statistics, 34, 206-215. See Also SemiParBIVProbit Examples ## see examples for SemiParBIVProbit LM.bpm Lagrange Multiplier Test (Score Test) Description Before fitting a bivariate probit model, LM.bpm can be used to test the hypothesis of absence of endogeneity, correlated model equations/errors or non-random sample selection. Usage LM.bpm(formula, data = list(), weights = NULL, subset = NULL, Model, hess = TRUE, infl.fac = 1, pPen1 = NULL, pPen2 = NULL) Arguments formula 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. Note that if Model = "BSS" then 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). 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. Model It indicates the type of model to be used in the analysis. Possible values are "B" (bivariate model) and "BSS" (bivariate model with sample selection). hess If FALSE then the expected (rather than observed) information matrix is employed. LM.bpm 13 infl.fac 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. pPen1, pPen2 Optional list specifying any penalties to be applied to the parametric model terms of equations 1 and 2. Details This Lagrange multiplier test (also known as score test) is used here for testing the null hypothesis that ρ is equal to 0 (i.e. no endogeneity, non-random sample selection or correlated model equations/errors, depending on the model being fitted). Its main advantage is that it does not require an estimate of the model parameter vector under the alternative hypothesis. Asymptotically, it takes a Chi-squared distribution with one degree of freedom. Full details can be found in Marra et al. (2014) and Marra et al. (submitted). Value It returns a numeric p-value corresponding to the null hypothesis that ρ is equal to 0. WARNINGS This test’s implementation is only valid for models with normal errors and probit links. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References Marra G., Radice R. and Filippou P. (submitted), Regression Spline Bivariate Probit Models: A Practical Approach to Testing for Exogeneity. 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), 715741. See Also SemiParBIVProbit Examples ## see examples for SemiParBIVProbit 14 logLik.SemiParBIVProbit logLik.SemiParBIVProbit Extract the log likelihood for a fitted SemiParBIVProbit Description It extracts the log-likelihood for a fitted SemiParBIVProbit model. Usage ## S3 method for class 'SemiParBIVProbit' logLik(object, ...) Arguments object A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). ... Un-used for this function. Details Modification of the classic logLik which accounts for the estimated degrees of freedom used in SemiParBIVProbit objects. This function is provided so that information criteria work correctly with SemiParBIVProbit objects by using the correct degrees of freedom. Value Standard logLik object. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> See Also AIC, BIC Examples ## see examples for SemiParBIVProbit mb 15 mb Nonparametric (worst-case and IV) Manski’s bounds Description mb can be used to calculate the (worst-case and IV) Manski’s bounds and confidence interval covering the true effect of interest with a fixed probability. Usage mb(treat, outc, IV = NULL, Model, B = 100, sig.lev = 0.05) Arguments treat outc IV Model B sig.lev Binary treatment/selection variable. Binary outcome variable. An instrumental binary variable can be used if available. Possible values are "B" (model with endogenous variable) and "BSS" (model with non-random sample selection). Number of bootstrap replicates. This is used to obtain some components needed for confidence interval calculations. Significance level. Details Based on Manski (1990), this function returns the nonparametric lower and upper (worst-case) Manski’s bounds for the average treatment effect (ATE) when Model = "B" or prevalence when Model = "BSS". When an IV is employed the function returns IV Manski bounds. For comparison, it also returns the estimated effect assuming random assignment (i.e., the treatment received or selection relies on the assumption of ignorable observed and unobserved selection). Note that this is different from what provided by AT or est.prev when naive = FALSE as observed confounders are accounted for and the assumption here is of ignorable unobserved selection. A confidence interval covering the true ATE/prevalence with a fixed probability is also provided. This is based on the approach described in Imbens and Manski (2004). NOTE that this interval is typically very close (if not identical) to the lower and upper bounds. The ATE can be at most 1 (or 100 in percentage) and the worst-case Manski’s bounds have width 1. This means that 0 is always included within the possibilites of these bounds. Nevertheless, this may be useful to check whether the effect from a bivariate recursive model is included within the possibilites of the bounds. When estimating the prevalance the worst-case Manski’s bounds have width equal to the nonresponse probability, which provides a measure of the uncertainty about the prevalence caused by non-response. Again, this may be useful to check whether the prevalence from a bivariate nonrandom sample selection model is included within the possibilites of the bounds. See SemiParBIVProbit for some examples. 16 penPL Value LB, UP CI av.p Lower and upper bounds for the true effect of interest. Confidence interval covering the true effect of interest with a fixed probability. Estimated effect of interest assuming random assignment. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References Manski C.F. (1990), Nonparametric Bounds on Treatment Effects. American Economic Review, Papers and Proceedings, 80(2), 319-323. Imbens G.W. and Manski C.F (2004), Confidence Intervals for Partially Identified Parameters. Econometrica, 72(6), 1845-1857. See Also SemiParBIVProbit Examples ## see examples for SemiParBIVProbit pen Internal Function Description It provides an overall penalty matrix in a format suitable for estimation conditional on smoothing parameters. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> penPL Internal Function Description It provides an overall penalty matrix in a format suitable for estimation conditional on smoothing parameters, when asymmetric links are used. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> plot.SemiParBIVProbit 17 plot.SemiParBIVProbit SemiParBIVProbit plotting Description It takes a fitted SemiParBIVProbit object produced by SemiParBIVProbit() and plots the component smooth functions that make it up on the scale of the linear predictor. This function is based on plot.gam() in mgcv. Please see the documentation of plot.gam() for full details. Usage ## S3 method for class 'SemiParBIVProbit' plot(x, eq, ...) Arguments x eq ... A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). The equation from which smooth terms should be considered for printing. Other graphics parameters to pass on to plotting commands, as described for plot.gam in mgcv. Details This function produces plot showing the smooth terms of a fitted semiparametric bivariate probit model. For plots of 1-D smooths, the x axis of each plot is labelled using the name of the regressor, while the y axis is labelled as s(regr, edf) where regr is the regressor name, and edf the effective degrees of freedom of the smooth. As for 2-D smooths, perspective plots are produced with the x axes labelled with the first and second variable names and the y axis is labelled as s(var1, var2, edf), which indicates the variables of which the term is a function and the edf for the term. If seWithMean = TRUE then the intervals include the uncertainty about the overall mean. That is, although each smooth is shown centred, the confidence intervals are obtained as if every other term in the model was constrained to have average 0 (average taken over the covariate values) except for the smooth being plotted. The theoretical arguments and simulation study of Marra and Wood (2012) suggests that seWithMean = TRUE results in intervals with close to nominal frequentist coverage probabilities. Value The function generates plots. WARNING The function can not deal with smooths of more than 2 variables. This function is currently not suitable when plotting on the response scale if asymmetric links are employed. 18 predict.SemiParBIVProbit Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References Marra G. and Wood S.N. (2012), Coverage Properties of Confidence Intervals for Generalized Additive Model Components. Scandinavian Journal of Statistics, 39(1), 53-74. See Also AT, est.prev, SemiParBIVProbit, summary.SemiParBIVProbit, predict.SemiParBIVProbit Examples ## see examples for SemiParBIVProbit predict.SemiParBIVProbit Prediction from a fitted SemiParBIVProbit model Description It takes a fitted SemiParBIVProbit object produced by SemiParBIVProbit() and, for each equation, produces predictions for a new set of values of the model covariates or the original values used for the model fit. Standard errors of predictions can be produced. These are based on the posterior distribution of the model coefficients. This function is based on predict.gam() in mgcv. Please see the documentation of predict.gam() for full details. Usage ## S3 method for class 'SemiParBIVProbit' predict(object, eq, ...) Arguments object A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). eq The equation to be considered for prediction. ... Other arguments as in predict.gam() in mgcv. WARNING This function is currently not suitable when plotting on the response scale if asymmetric links are employed. print.AT 19 Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> See Also SemiParBIVProbit, plot.SemiParBIVProbit, summary.SemiParBIVProbit print.AT Print an AT object Description The print method for an AT object. Usage ## S3 method for class 'AT' print(x, ...) Arguments x AT object produced by AT(). ... Other arguments. Details print.AT prints the lower confidence interval limit, estimated AT and upper confidence interval limit. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> See Also AT 20 print.mb print.est.prev Print an est.prev object Description The print method for an est.prev object. Usage ## S3 method for class 'est.prev' print(x, ...) Arguments x est.prev object produced by est.prev(). ... Other arguments. Details print.est.prev prints the lower confidence interval limit, estimated prevalence and upper confidence interval limit. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> See Also est.prev print.mb Print an mb object Description The print method for an mb object. print.SemiParBIVProbit 21 Usage ## S3 method for class 'mb' print(x, ...) Arguments x mb object produced by mb(). ... Other arguments. Details print.mb prints the lower and upper bounds, confidence interval, and effect assuming random assignment. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> See Also mb print.SemiParBIVProbit Print a SemiParBIVProbit object Description The print method for a SemiParBIVProbit object. Usage ## S3 method for class 'SemiParBIVProbit' print(x, ...) Arguments x SemiParBIVProbit object produced by SemiParBIVProbit(). ... Other arguments. 22 print.summary.SemiParBIVProbit Details print.SemiParBIVProbit prints out the family, model equations, total number of observations, estimated association coefficient and total effective degrees of freedom for the penalized or unpenalized model. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> See Also SemiParBIVProbit print.summary.SemiParBIVProbit Print a summary.SemiParBIVProbit object Description The print method for a summary.SemiParBIVProbit object. Usage ## S3 method for class 'summary.SemiParBIVProbit' print(x, digits = max(3, getOption("digits") - 3), signif.stars = getOption("show.signif.stars"), ...) Arguments x summary.SemiParBIVProbit object produced by summary.SemiParBIVProbit(). digits Number of digits printed in output. signif.stars By default significance stars are printed alongside output. ... Other arguments. Details print.summary.SemiParBIVProbit prints model term summaries. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> S.m 23 See Also summary.SemiParBIVProbit S.m Internal Function Description It provides penalty matrices in a format suitable for automatic smoothness estimation. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> S.mPL Internal Function Description It provides penalty matrices in a format suitable for automatic smoothness estimation, when an asymmetric link approach is used. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> sem.checks Some convergence diagnostics Description It takes a fitted SemiParBIVProbit object produced by SemiParBIVProbit() and produces some diagnostic information about the fitting procedure. Usage sem.checks(x) Arguments x SemiParBIVProbit object produced by SemiParBIVProbit(). 24 SemiParBIVProbit Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> See Also SemiParBIVProbit 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. Several bivariate copula distributions are supported. During the model fitting process, the possible presence of correlated error equations, endogeneity, non-random sample selection or partial observability 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). Usage SemiParBIVProbit(formula, data = list(), weights = NULL, subset = NULL, start.v = NULL, Model = "B", method = "trwlF", BivD = "N", nu = 3, fp = FALSE, PL = "P", eqPL = "both", valPL = c(0,0), fitPL = "pLiksp", spPL = c(0.01,0.01), hess = TRUE, infl.fac = 1, pPen1 = NULL, pPen2 = NULL, rinit = 1, rmax = 100, iterlimsp = 50, pr.tolsp = 1e-6, gc.l = FALSE, awlm = FALSE, parscale, extra.regI = "t") Arguments formula 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. SemiParBIVProbit supports the use shrinkage smoothers for variable selection purposes. See the examples below and the documentation of mgcv for further details on formula specifications. Note that if Model = "BSS" then 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 SemiParBIVProbit 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. start.v Starting values for the parameters of the two equations and association coefficient can be provided here. SemiParBIVProbit 25 Model It indicates the type of model to be used in the analysis. Possible values are "B" (bivariate model), "BSS" (bivariate model with non-random sample selection) and "BPO" (bivariate model with partial observability). method Estimation method used. It can be either "trwl" or "trwlF". The latter option is much faster and more stable when smoothing parameter estimation is employed, and is the default. 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), 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, and Student-t (with fixed degrees of freedom nu). Note that for empirical applications options N, C (and all rotated versions), F and T may suffice. This is becuase the Clayton, Joe and Gumbel copulas are very similar; for instace, C90, J270 and G270 are very close. Nevertheless, there might be situations in which the use of Joe or Gumbel may result in more stable computation as compared to Clayton. nu Degrees of freedom for Student-t. fp If TRUE then a fully parametric model with unpenalised regression splines if fitted. See the example below. PL Link function to be employed. Possible choices are "P" (classic probit - default option), "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 fixed smoothing parameters and penalized likelihood with automatic smoothing parameter estimation. spPL Values for the smoothing parameters of the ridge penalties associated with the shape parameteres of the asymmetric link functions. hess If FALSE then the expected/Fisher (rather than observed) information matrix is employed. infl.fac 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. 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. 26 SemiParBIVProbit 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 estimation is used. gc.l This is relevant when working with big datasets. If TRUE then the garbage collector is called more often than when R decides to. This keeps the memory footprint down but it will slow down the routine. awlm If TRUE then the association parameter of the two equations is used in the construction of the working linear model needed for smoothing parameter estimation. This option makes sense only when method = "trwl". parscale The algorithm will operate as if optimizing objfun(x / parscale, ...). 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 only makes sense when hess = FALSE) and "sED" (symmetric eigen-decomposition). Details The bivariate models considered in this package have probit links for the two model equations, and model the association between the responses through the correlation parameter ρ of a standardised bivariate normal distribution, or that of a bivariate copula distribution, θ. 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 is 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), which can also be viewed as an approximate AIC. Automatic smoothing parameter estimation 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). Value The function returns an object of class SemiParBIVProbit as described in SemiParBIVProbitObject. WARNINGS In our experience, convergence failure may occur, for instance, when the number of observations is low, there are important mistakes in the model specification (i.e., using C90 when the model equations are positively associated), or there are many smooth components in the model as compared to the number of observations. Convergence failure may also mean that an infinite cycling SemiParBIVProbit 27 between steps (i) and (ii) occurs. In this case, the smoothing parameters are set to the values obtained from the non-converged algorithm (sem.checks will give a warning). In such cases, we recommend re-specifying the model, using some extra regularisation (see extra.regI) and/or using some rescaling (see parscale). The use of asymmetric links is not likely to lead to significant improvements as compared to using probit links. We recommend using this option at the end of the empirical analysis as a sensitivity analysis tool, and not during the model building process. Also, the theoretical and numerical complexity introduced by asymmetric links means that convergence failures are more likely to occur and that data requirements are typically high. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References Key References: Gu C. (1992), Cross validating non-Gaussian data. Journal of Computational and Graphical Statistics, 1(2), 169-179. Marra G. (2013), On P-values for Semiparametric Bivariate Probit Models. Statistical Methodology, 10(1), 23-28. 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. (2015), Flexible Bivariate Binary Models for Estimating the Efficacy of Phototherapy for Newborns with Jaundice. International Journal of Statistics and Probability. 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), 715741. Marra G., Radice R. and Filippou P. (submitted), Regression Spline Bivariate Probit Models: A Practical Approach to Testing for Exogeneity. McGovern M.E., Barnighausen T., Marra G. and Radice R. (2015), On the Assumption of Joint Normality in Selection Models: A Copula Based Approach Applied to Estimating HIV Prevalence. Epidemiology. Radice R., Marra G. and M. Wojtys (submitted), Copula Regression Spline Models for Binary Outcomes. Marra G. and Wood S.N. (2011), Practical Variable Selection for Generalized Additive Models. Computational Statistics and Data Analysis, 55(7), 2372-2387. Marra G. and Wood S.N. (2012), Coverage Properties of Confidence Intervals for Generalized Additive Model Components. Scandinavian Journal of Statistics, 39(1), 53-74. Poirier D.J. (1980), Partial Observability in Bivariate Probit Models. Journal of Econometrics, 12, 209-217. 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. 28 SemiParBIVProbit See Also AT, adjCov, est.prev, gt.bpm, LM.bpm, VuongClarke.bcm, plot.SemiParBIVProbit, SemiParBIVProbit-package, SemiParBIVProbitObject, sem.checks, summary.SemiParBIVProbit, predict.SemiParBIVProbit Examples 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 f2 <- function(x) cos(pi*2*x) + sin(pi*x) <- 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(list(y1 ~ x1 + x2 + x3, y2 ~ x1 + x2 + x3), data = dataSim) sem.checks(out) summary(out) AIC(out) BIC(out) ## SEMIPARAMETRIC BIVARIATE PROBIT ## ## ## ## ## "cr" cubic regression spline basis - "cs" shrinkage version "tp" thin plate regression spline basis - "ts" shrinkage version for smooths of one variable, "cr/cs" and "tp/ts" achieve similar k is the basis dimension - default is 10 m is the order of the penalty for the specific term - default is of "cr" of "tp" results 2 SemiParBIVProbit 29 ## For COPULA models use BivD argument out <- SemiParBIVProbit(list(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, cplot = TRUE) AIC(out) ## estimated smooth function plots - red lines are true curves x2 <f1.x2 f2.x2 f3.x3 sort(x2) <- f1(x2)[order(x2)] - mean(f1(x2)) <- f2(x2)[order(x2)] - mean(f2(x2)) <- rep(0, length(x3)) par(mfrow=c(2,2),mar=c(4.5,4.5,2,2)) plot(out, eq = 1, select = 1, scale = plot(out, eq = 1, select = 2, scale = plot(out, eq = 2, select = 1, scale = plot(out, eq = 2, select = 2, scale = 0); 0); 0); 0); lines(x2, lines(x3, lines(x2, lines(x3, f1.x2, f3.x3, f2.x2, f3.x3, col col col col = = = = "red") "red") "red") "red") # ## same plots but CIs 'with intercept' plot(out, plot(out, plot(out, plot(out, eq eq eq eq = = = = 1, 1, 2, 2, select select select select = = = = 1, 2, 1, 2, seWithMean seWithMean seWithMean seWithMean = = = = TRUE, TRUE, TRUE, TRUE, scale scale scale scale = = = = 0); 0); 0); 0); lines(x2, lines(x3, lines(x2, lines(x3, f1.x2, f3.x3, f2.x2, f3.x3, col col col col ## 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(list(y1 ~ x1 + s(x2, bs = "cs") + s(x3, bs = "cs"), y2 ~ x1 + s(x2, bs = "cs") + s(x3, bs = "cs")), data = dataSim) sem.checks(outSS) plot(outSS, eq plot(outSS, eq plot(outSS, eq plot(outSS, eq graphics.off() = = = = 1, 1, 2, 2, select select select select = = = = 1, 2, 1, 2, scale = 0) ylim = c(-0.1,0.1)) scale = 0) ylim = c(-0.1,0.1)) # # ############ ## EXAMPLE 2 ############ ## Generate data with one endogenous variable = = = = "red") "red") "red") "red") 30 SemiParBIVProbit ## and exclusion restriction ## Not run: 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 f2 <- function(x) cos(pi*2*x) + sin(pi*x) <- 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(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), dataSim, Model = "B") # p-value suggests presence of endogeneity, hence fit a bivariate model ## CLASSIC RECURSIVE BIVARIATE PROBIT out <- SemiParBIVProbit(list(y1 ~ x1 + x2, y2 ~ y1 + x2), data = dataSim) sem.checks(out) summary(out) ## SEMIPARAMETRIC RECURSIVE BIVARIATE PROBIT out <- SemiParBIVProbit(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), data = dataSim) sem.checks(out) summary(out) # ## Testing the hypothesis of absence of endogeneity post estimation... gt.bpm(out) SemiParBIVProbit # ## average treatment effect with CI mb(y1, y2, Model = "B") AT(out, eq = 2, nm.bin = "y1") AT(out, eq = 2, nm.bin = "y1", naive = TRUE) ## try a Clayton copula model... outC <- SemiParBIVProbit(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), data = dataSim, BivD = "C0") sem.checks(outC) summary(outC) AT(outC, eq = 2, nm.bin = "y1") ## try a Joe copula model... outJ <- SemiParBIVProbit(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), data = dataSim, BivD = "J0") sem.checks(outJ) summary(outJ, cplot = TRUE) 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(list(y1 ~ x1 + s(x2, bs = "cr", k = 5), y2 ~ y1 + s(x2, bs = "cr", k = 6)), fp = TRUE, data = dataSim) sem.checks(outFP) summary(outFP) # ## SEMIPARAMETRIC RECURSIVE BIVARIATE PROBIT ## WITH ASYMMETRIC LINK FUNCTIONS out <- SemiParBIVProbit(list(y1 ~ x1 + s(x2), y2 ~ y1 + s(x2)), data = dataSim, PL = "SN", hess = FALSE) sem.checks(out) summary(out) # same output as model with probit links ############ ## EXAMPLE 3 31 32 SemiParBIVProbit ############ ## 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(list(ys ~ bi + s(x1) + s(x2), yo ~ bi + s(x1)), dataSim, Model = "BSS") # 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(list(ys ~ bi + s(x1) + s(x2), yo ~ bi + s(x1)), data = dataSim, Model = "BSS") sem.checks(out) 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 mb(ys, yo, Model = "BSS") est.prev(out) SemiParBIVProbit 33 est.prev(out, naive = TRUE) ## 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.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(list(ys ~ bi + s(x1) + s(x2), yo ~ bi + s(x1)), data = dataSim, Model = "BSS", BivD = "C0") sem.checks(outC) summary(outC, cplot = TRUE) est.prev(outC) # ############ ## EXAMPLE 4 ############ ## Generate data with partial observability set.seed(0) n <- 10000 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) y1 <- ifelse(-1.55 + 2*x1 + x2 + u[,1] > 0, 1, 0) y2 <- ifelse( 0.45 - x3 + u[,2] > 0, 1, 0) y <- y1*y2 dataSim <- data.frame(y, x1, x2, x3) ## BIVARIATE PROBIT with Partial Observability out <- SemiParBIVProbit(list(y ~ x1 + x2, y ~ x3), data = dataSim, Model = "BPO") sem.checks(out) 34 SemiParBIVProbit.fit summary(out) outC <- SemiParBIVProbit(list(y ~ x1 + x2, y ~ x3), data = dataSim, Model = "BPO", BivD = "C0") sem.checks(outC) summary(outC) # first ten estimated probabilities for the four events from object out cbind(out$p11, out$p10, out$p00, out$p01)[1:10,] # case with smooth function # (a bit more computationally intensive) f1 <- function(x) cos(pi*2*x) + sin(pi*x) y1 <- ifelse(-1.55 + 2*x1 + f1(x2) + u[,1] > 0, 1, 0) y2 <- ifelse( 0.45 - x3 + u[,2] > 0, 1, 0) y <- y1*y2 dataSim <- data.frame(y, x1, x2, x3) out <- SemiParBIVProbit(list(y ~ x1 + s(x2), y ~ x3), data = dataSim, Model = "BPO") sem.checks(out) summary(out, cplot = TRUE) # plot estimated and true functions x2 <- sort(x2); f1.x2 <- f1(x2)[order(x2)] - mean(f1(x2)) plot(out, eq = 1, scale = 0); lines(x2, f1.x2, col = "red") # # ## End(Not run) SemiParBIVProbit.fit Description Internal fitting function. Internal Function SemiParBIVProbit.fit.post 35 Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> SemiParBIVProbit.fit.post Internal Function Description Internal post estimation function. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> SemiParBIVProbit.fit1 Internal Function Description Internal fitting function. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> SemiParBIVProbitObject Fitted SemiParBIVProbit object Description A fitted semiparametric bivariate probit object returned by function SemiParBIVProbit and of class "SemiParBIVProbit". 36 SemiParBIVProbitObject Value fit List of values and diagnostics extracted from the output of the algorithm. For instance, fit$gradient, fit$Fisher and fit$S.h return the gradient vector, Fisher information (when used) and overall penalty matrix scaled by its smoothing parameters, for the fitted bivariate probit model. See the documentation of trust for details on the diagnostics provided. gam1 Univariate fit for equation 1. See the documentation of mgcv for full details. gam2 Univariate fit for equation 2. gam2.1 Univariate fit for equation 2, estimated using an adaptation of the Heckman sample selection correction procedure when Method = "BSS". coefficients The coefficients of the fitted model. They are given in the following order: parametric and regression spline (if present) coefficients for the first equation, parametric and regression spline coefficients for the second equation, and association coefficient. weights Prior weights used during model fitting. sp Estimated smoothing parameters of the smooth components. iter.sp Number of iterations performed for the smoothing parameter estimation step. iter.if Number of iterations performed in the initial step of the algorithm. iter.inner Number of iterations performed within the smoothing parameter estimation step. rho Estimated correlation coefficient between the two equations, when a bivariate normal model is employed. theta Estimated copula parameter linking the two equations. n Sample size. n.sel Number of selected observations in the sample selection model case. X1, X2 Design matrices associated with the first and second linear predictors. X1.d2, X2.d2 Number of columns of X1 and X2. l.sp1, l.sp2 Number of smooth components in equations 1 and 2. He Penalized -hessian/Fisher. This is the same as HeSh for unpenalized models. HeSh Unpenalized -hessian/Fisher. Vb Inverse of He. This corresponds to the Bayesian variance-covariance matrix used for confidence/credible interval calculations. F This is given by Vb*HeSh. t.edf Total degrees of freedom of the estimated bivariate model. It is calculated as sum(diag(F)). edf1, edf2 Degrees of freedom for the two equations of the fitted bivariate model. It is calculated when splines are used. bs.mgfit List of values and diagnostics extracted from magic in mgcv. conv.sp If TRUE then the smoothing parameter selection algorithm stopped before reaching the maximum number of iterations allowed. wor.c Working model quantities. startSS 37 p11, p10, p01, p00 Model probabilities evaluated at (y_1 = 1, y_2 = 1), (y_1 = 1, y_2 = 0), (y_1 = 0, y_2 = 1) and (y_1 = 0, y_2 = 0). p0 Distribution function of a standardised univariate normal evaluated at y_1 = 0. This is only provided when Method = "BSS". eta1, eta2 Estimated linear predictors for the two equations. y1, y2 Responses of the two equations. xi1, xi2 Estimated shape parameters of the link functions of the two equations if an asymmetric link approach is used. logLik Value of the (unpenalized) log-likelihood evaluated at the (penalized or unpenalized) parameter estimates. pPen1, pPen2 List specifying any penalties to be applied to the parametric model terms of the model equations. good Indicator variable indicating the observations actually used in model fitting. respvec List containing response vectors. X2s Full design matrix of outcome equation in sample selection case. OR, GM Odds ratio and Gamma measure. See summary.SemiParBIVProbit for details. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> See Also SemiParBIVProbit, plot.SemiParBIVProbit, summary.SemiParBIVProbit, predict.SemiParBIVProbit startSS Internal Function Description It calculates sensible starting values for bivariate probit models with non-random sample selection. This is based on an adaptation of the Heckman two-step procedure. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References 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. 38 summary.SemiParBIVProbit summary.SemiParBIVProbit SemiParBIVProbit summary Description It takes a fitted SemiParBIVProbit object produced by SemiParBIVProbit() and produces some summaries from it. Usage ## S3 method for class 'SemiParBIVProbit' summary(object, n.sim = 100, s.meth = "svd", prob.lev = 0.05, thrs1 = 0.5, thrs2 = 0.5, cplot = FALSE, ...) Arguments object A fitted SemiParBIVProbit object as produced by SemiParBIVProbit(). n.sim The number of simulated coefficient vectors from the posterior distribution of the estimated model parameters. This is used to calculate intervals for the association parameter and other measures (e.g., gamma measure). It should be increased if more precision is required. s.meth Matrix decomposition used to determine the matrix root of the covariance matrix. See the documentation of mvtnorm for further details. prob.lev Probability of the left and right tails of the posterior distribution used for interval calculations. thrs1, thrs2 Threshold to calculate the correct prediction ratios for the two binary responses. Default value is 0.5. It only works when Method = "B". cplot If TRUE then a bivariate contour meta plot corresponding to the assumed bivariate model with estimated association parameter is produced. See documentation of CDVine for details. ... Other graphics parameters to pass on to plotting commands, when cplot=TRUE. Details As in the package mgcv, based on the results of Marra and Wood (2012), ‘Bayesian p-values’ are returned for the smooth terms. These have better frequentist performance than their frequentist counterpart. Let ˆf and Vf denote the vector of values of a smooth term evaluated at the original covariate values and the corresponding Bayesian covariance matrix, and let Vfr− denote the rank r pseudoinverse of Vf . The statistic used is T = ˆf 0 Vfr− ˆf . This is compared to a chi-squared distribution with degrees of freedom given by r, which is obtained by biased rounding of the estimated summary.SemiParBIVProbit 39 degrees of freedom. See Marra (2013) for further details. Covariate selection can also be achieved using a single penalty shrinkage approach as shown in Marra and Wood (2011). Posterior simulation is used to obtain intervals of nonlinear functions of parameters, such as the odds ratio and gamma measure discussed by Tajar et al. (2001). For more details on some of the model summaries see Radice, Marra and Wojtys (submitted). Value tableP1 Table containing parametric estimates, their standard errors, z-values and pvalues for equation 1. tableP2 As above but for equation 2. tableNP1 Table of nonparametric summaries for each smooth component including effective degrees of freedom, estimated rank, approximate Wald statistic for testing the null hypothesis that the smooth term is zero and corresponding p-value, for equation 1. tableNP2 As above but for equation 2. n Sample size. rho Estimated correlation parameter between the two equations. theta Estimated copula parameter linking the two equations. formula1, formula2 Formulas used for equations 1 and 2. l.sc1, l.sc2 Number of smooth components in equations 1 and 2. t.edf Total degrees of freedom of the estimated bivariate model. CIrs ‘Confidence’ interval for either ρ or θ, depending of the model fitted. CIl1, CIl2 ‘Confidence’ intervals for the shape parameters when an asymmetric link approach is employed. n.sel Number of selected observations in the sample selection case. good Indicator variable indicating the observations actually used in model fitting. OR, CIor Odds ratio and related CI. The odds ratio is a measure of association between binary random variables and is defined as p00p11/p10p01. In the case of independence this ratio is equal to 1. It can take values in the range (-Inf, Inf) and it does not depend on the marginal probabilities (Tajar et al., 2001). Interval is calculated using posterior simulation. GM, CIgm Gamma measure and related CI. This measure of association was proposed by Goodman and Kruskal (1954). It is defined as (OR - 1)/(OR + 1), can take values in the range (-1, 1) and does not depend on the marginal probabilities. Interval is calculated using posterior simulation. xi1, xi2 Estimated power or shape parameters of the link functions of the two equations if an asymmetric link approach is used. table.R Joint absolute frequency table of observed binary value combinations. This as well as all remaining quantities below are not provided when Method = "BSS" or Method = "BPO". table.P Joint absolute frequency table of predicted binary value combinations. 40 summary.SemiParBIVProbit table.F It is given by table.P/table.R. MR Based on the comparison between table.P and table.R, it reports the matching rate as percentage. P1, P2 Marginal probabilities for the two binary outcomes y_1 = 1 and y_2 = 1. QPS1, QPS2 Quadratic probability scores for the two binary responses as suggested by Diebold and Rudebusch (1989). Values are on the interval [0,2], with 0 indicating a perfect fit. CR1, CR2 Correct prediction ratios for the two binary responses as percentage. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References Diebold F.X. and Rudebusch G.D. (1989), Scoring the Leading Indicators. Journal of Business, 62(3), 369-391. Goodman L.A. and Kruskal W.H. (1954), Measures of Association for Cross Classification. Journal of the American Statistical Association, 49, 732-764. Marra G. (2013), On P-values for Semiparametric Bivariate Probit Models. Statistical Methodology, 10(1), 23-28. Radice R., Marra G. and M. Wojtys (submitted), Copula Regression Spline Models for Binary Outcomes. Marra G. and Wood S.N. (2011), Practical Variable Selection for Generalized Additive Models. Computational Statistics and Data Analysis, 55(7), 2372-2387. Marra G. and Wood S.N. (2012), Coverage Properties of Confidence Intervals for Generalized Additive Model Components. Scandinavian Journal of Statistics, 39(1), 53-74. Tajar M., Denuit M. and Lambert P. (2001), Copula-Type Representation for Random Couples with Bernoulli Margins. Discussion Papaer 0118, Universite Catholique De Louvain. See Also AT, est.prev, SemiParBIVProbitObject, plot.SemiParBIVProbit, predict.SemiParBIVProbit Examples ## see examples for SemiParBIVProbit VuongClarke.bcm VuongClarke.bcm 41 Vuong and Clarke tests Description The Vuong and Clarke tests are likelihood-ratio-based tests that can be used for choosing between two non-nested models. Usage VuongClarke.bcm(obj1, obj2, sig.lev = 0.05) Arguments obj1,obj2 Objects of the two fitted bivariate non-nested models. sig.lev Significance level used for testing. Details The Vuong (1989) and Clarke (2007) tests are likelihood-ratio-based tests for model selection that use the Kullback-Leibler information criterion. The implemented tests can be used for choosing between two bivariate models which are non-nested. In the Vuong test, the null hypothesis is that the two models are equally close to the actual model, whereas the alternative is that one model is closer. The test follows asymptotically a standard normal distribution under the null. Assume that the critical region is (−c, c), where c is typically set to 1.96. If the value of the test is higher than c then we reject the null hypothesis that the models are equivalent in favor of model obj1. Viceversa if the value is smaller than c. If the value falls in [−c, c] then we cannot discriminate between the two competing models given the data. In the Clarke test, if the two models are statistically equivalent then the log-likelihood ratios of the observations should be evenly distributed around zero and around half of the ratios should be larger than zero. The test follows asymptotically a binomial distribution with parameters n and 0.5. Critical values can be obtained as shown in Clarke (2007). Intuitively, model obj1 is preferred over obj2 if the value of the test is significantly larger than its expected value under the null hypothesis (n/2), and vice versa. If the value is not significantly different from n/2 then obj1 can be thought of as equivalent to obj2. Value It returns two decisions based on the tests and criteria discussed above. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> 42 working.comp References Clarke K. (2007), A Simple Distribution-Free Test for Non-Nested Model Selection. Political Analysis, 15, 347-363. Vuong Q.H. (1989), Likelihood Ratio Tests for Model Selection and Non-Nested Hypotheses. Econometrica, 57(2), 307-333. See Also SemiParBIVProbit Examples ## see examples for SemiParBIVProbit working.comp Internal Function Description It efficiently calculates the working model quantities needed to implement the automatic multiple smoothing parameter estimation procedure by exploiting the band structure of the matrices involved in the calculations. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References 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. Radice R., Marra G. and M. Wojtys (submitted), Copula Regression Spline Models for Binary Outcomes. working.comp1 working.comp1 43 Internal Function Description It efficiently calculates the working model quantities needed to implement the automatic multiple smoothing parameter estimation procedure by exploiting a result which leads to very fast and stable calculations. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> working.comp3 Internal Function Description It efficiently calculates the working model quantities needed to implement the automatic multiple smoothing parameter estimation procedure by exploiting the band structure of the matrices involved in the calculations. This function includes the association parameter of the two equations in the construction of the working linear model. Author(s) Maintainer: Giampiero Marra <giampiero.marra@ucl.ac.uk> References Radice R., Marra G. and M. Wojtys (submitted), Copula Regression Spline Models for Binary Outcomes. Index LM.bpm, 12 ∗Topic covariance matrix adjustment adjCov, 4 adjCovSD, 5 ∗Topic delta method AT, 6 est.prev, 10 ∗Topic diagnostics sem.checks, 23 ∗Topic endogeneity gt.bpm, 11 LM.bpm, 12 SemiParBIVProbit, 24 SemiParBIVProbit-package, 2 VuongClarke.bcm, 41 ∗Topic gradient test gt.bpm, 11 ∗Topic hplot plot.SemiParBIVProbit, 17 ∗Topic information criteria summary.SemiParBIVProbit, 38 ∗Topic lagrange multiplier test LM.bpm, 12 ∗Topic likelihood ratio test VuongClarke.bcm, 41 ∗Topic logLik logLik.SemiParBIVProbit, 14 ∗Topic non-random sample selection est.prev, 10 gt.bpm, 11 LM.bpm, 12 SemiParBIVProbit, 24 SemiParBIVProbit-package, 2 VuongClarke.bcm, 41 ∗Topic package SemiParBIVProbit-package, 2 ∗Topic partial observability SemiParBIVProbit, 24 SemiParBIVProbit-package, 2 ∗Topic AIC logLik.SemiParBIVProbit, 14 ∗Topic ATE AT, 6 mb, 15 ∗Topic BIC logLik.SemiParBIVProbit, 14 ∗Topic Clarke test VuongClarke.bcm, 41 ∗Topic Manski’s bounds mb, 15 ∗Topic Nonparametric bounds mb, 15 ∗Topic Vuong test VuongClarke.bcm, 41 ∗Topic Worst-case bounds mb, 15 ∗Topic asymmetric link function SemiParBIVProbit, 24 ∗Topic average treatment effect AT, 6 mb, 15 ∗Topic bayesian posterior simulation AT, 6 ∗Topic bivariate model VuongClarke.bcm, 41 ∗Topic bivariate probit model AT, 6 est.prev, 10 gt.bpm, 11 LM.bpm, 12 SemiParBIVProbit, 24 ∗Topic complex survey design adjCovSD, 5 ∗Topic confidence interval mb, 15 ∗Topic correlated model equations/errors gt.bpm, 11 44 INDEX 45 ∗Topic prediction predict.SemiParBIVProbit, 18 ∗Topic prevalence est.prev, 10 mb, 15 ∗Topic regression spline SemiParBIVProbit, 24 ∗Topic regression plot.SemiParBIVProbit, 17 SemiParBIVProbit-package, 2 summary.SemiParBIVProbit, 38 ∗Topic score test LM.bpm, 12 ∗Topic semiparametric bivariate probit modelling AT, 6 est.prev, 10 gt.bpm, 11 LM.bpm, 12 print.AT, 19 print.est.prev, 20 print.mb, 20 print.SemiParBIVProbit, 21 print.summary.SemiParBIVProbit, 22 sem.checks, 23 SemiParBIVProbit, 24 SemiParBIVProbit-package, 2 VuongClarke.bcm, 41 ∗Topic shrinkage smoother SemiParBIVProbit, 24 summary.SemiParBIVProbit, 38 ∗Topic smooth plot.SemiParBIVProbit, 17 SemiParBIVProbit, 24 SemiParBIVProbit-package, 2 summary.SemiParBIVProbit, 38 ∗Topic variable selection SemiParBIVProbit, 24 SemiParBIVProbit-package, 2 summary.SemiParBIVProbit, 38 adjCov, 4, 28 adjCovSD, 5 AIC, 14 AT, 3, 6, 15, 18, 19, 28, 40 BIC, 14 bprobgHs, 8 bprobgHsPL, 8 bprobgHsPO, 9 bprobgHsSS, 9 copgHs, 10 est.prev, 3, 10, 15, 18, 20, 28, 40 gt.bpm, 3, 11, 28 LM.bpm, 3, 12, 28 logLik, 14 logLik.SemiParBIVProbit, 14 mb, 3, 15, 21 pen, 16 penPL, 16 plot.SemiParBIVProbit, 3, 17, 19, 28, 37, 40 predict.SemiParBIVProbit, 18, 18, 28, 37, 40 print.AT, 19 print.est.prev, 20 print.mb, 20 print.SemiParBIVProbit, 21 print.summary.SemiParBIVProbit, 22 S.m, 23 S.mPL, 23 sem.checks, 23, 27, 28 SemiParBIVProbit, 3–7, 11–13, 15, 16, 18, 19, 22, 24, 24, 37, 42 SemiParBIVProbit-package, 2 SemiParBIVProbit.fit, 34 SemiParBIVProbit.fit.post, 35 SemiParBIVProbit.fit1, 35 SemiParBIVProbitObject, 28, 35, 40 startSS, 37 summary.SemiParBIVProbit, 3, 5–7, 11, 18, 19, 23, 28, 37, 38 VuongClarke.bcm, 28, 41 working.comp, 42 working.comp1, 43 working.comp3, 43
© Copyright 2024