Nothing
R/binstest.R
In binsreg: Binscatter Estimation and Inference
Defines functions
summary.CCFFbinstest
print.CCFFbinstest
binstest
Documented in
binstest
print.CCFFbinstest
summary.CCFFbinstest
################################################################################################################
#'@title Data-Driven Nonparametric Shape Restriction and Parametric Model Specification Testing using Binscatter
#'@description \code{binstest} implements binscatter-based hypothesis testing procedures for parametric functional
#' forms of and nonparametric shape restrictions on the regression function of interest, following the results
#' in \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_AER.pdf}{Cattaneo, Crump, Farrell and Feng (2024a)} and
#' \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_NonlinearBinscatter.pdf}{Cattaneo, Crump, Farrell and Feng (2024b)}.
#' If the binning scheme is not set by the user,
#' the companion function \code{\link{binsregselect}} is used to implement binscatter in a
#' data-driven way and inference procedures are based on robust bias correction.
#' Binned scatter plots based on different methods can be constructed using the companion functions \code{\link{binsreg}},
#' \code{\link{binsqreg}} or \code{\link{binsglm}}.
#'@param y outcome variable. A vector.
#'@param x independent variable of interest. A vector.
#'@param w control variables. A matrix, a vector or a \code{\link{formula}}.
#'@param data an optional data frame containing variables used in the model.
#'@param estmethod estimation method. The default is \code{estmethod="reg"} for tests based on binscatter least squares regression. Other options are \code{"qreg"} for quantile regression and \code{"glm"} for generalized linear regression. If \code{estmethod="glm"}, the option \code{family} must be specified.
#'@param family a description of the error distribution and link function to be used in the generalized linear model when \code{estmethod="glm"}. (See \code{\link{family}} for details of family functions.)
#'@param quantile the quantile to be estimated. A number strictly between 0 and 1.
#'@param deriv derivative order of the regression function for estimation, testing and plotting.
#' The default is \code{deriv=0}, which corresponds to the function itself.
#'@param at value of \code{w} at which the estimated function is evaluated. The default is \code{at="mean"}, which corresponds to
#' the mean of \code{w}. Other options are: \code{at="median"} for the median of \code{w}, \code{at="zero"} for a vector of zeros.
#' \code{at} can also be a vector of the same length as the number of columns of \code{w} (if \code{w} is a matrix) or a data frame containing the same variables as specified in \code{w} (when
#' \code{data} is specified). Note that when \code{at="mean"} or \code{at="median"}, all factor variables (if specified) are excluded from the evaluation (set as zero).
#'@param nolink if true, the function within the inverse link function is reported instead of the conditional mean function for the outcome.
#'@param testmodel a vector or a logical value. It sets the degree of polynomial and the number of smoothness constraints for parametric model specification
#' testing. If \code{testmodel=c(p,s)} is specified, a piecewise polynomial of degree \code{p} with \code{s} smoothness constraints is used.
#' If \code{testmodel=T} or \code{testmodel=NULL} (default) is specified, \code{testmodel=c(1,1)} is used unless the degree \code{p} or the smoothness \code{s}
#' selection is requested via the option \code{pselect} or \code{sselect} (see more details in the explanation of \code{pselect} and \code{sselect}).
#'@param testmodelparfit a data frame or matrix which contains the evaluation grid and fitted values of the model(s) to be
#' tested against. The column contains a series of evaluation points
#' at which the binscatter model and the parametric model of interest are compared with
#' each other. Each parametric model is represented by other columns, which must
#' contain the fitted values at the corresponding evaluation points.
#'@param testmodelpoly degree of a global polynomial model to be tested against.
#'@param testshape a vector or a logical value. It sets the degree of polynomial and the number of smoothness constraints for nonparametric shape restriction
#' testing. If \code{testshape=c(p,s)} is specified, a piecewise polynomial of degree \code{p} with \code{s} smoothness constraints is used.
#' If \code{testshape=T} or \code{testshape=NULL} (default) is specified, \code{testshape=c(1,1)} is used unless the degree \code{p} or smoothness \code{s} selection
#' is requested via the option \code{pselect} or \code{sselect} (see more details in the explanation of \code{pselect} and \code{sselect}).
#'@param testshapel a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#' corresponds to one boundary of a one-sided hypothesis test to the left of the form
#' \code{H0: sup_x mu(x)<=a}.
#'@param testshaper a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#' corresponds to one boundary of a one-sided hypothesis test to the right of the form
#' \code{H0: inf_x mu(x)>=a}.
#'@param testshape2 a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#' corresponds to one boundary of a two-sided hypothesis test of the form
#' \code{H0: sup_x |mu(x)-a|=0}.
#'@param lp an Lp metric used for parametric model specification testing and/or shape restriction testing. The default is \code{lp=Inf}, which
#' corresponds to the sup-norm of the t-statistic. Other options are \code{lp=q} for a positive number \code{q>=1}. Note that
#' \code{lp=Inf} ("sup-norm") has to be used for testing one-sided shape restrictions.
#'@param bins a vector. If \code{bins=c(p,s)}, it sets the piecewise polynomial of degree \code{p} with \code{s} smoothness constraints
#' for data-driven (IMSE-optimal) selection of the partitioning/binning scheme.
#' The default is \code{bins=c(0,0)}, which corresponds to the piecewise constant.
#'@param nbins number of bins for partitioning/binning of \code{x}. If \code{nbins=T} or \code{nbins=NULL} (default) is specified, the number
#' of bins is selected via the companion command \code{\link{binsregselect}} in a data-driven, optimal way whenever possible.
#' If a vector with more than one number is specified, the number of bins is selected within this vector via the companion command \code{\link{binsregselect}}.
#'@param pselect vector of numbers within which the degree of polynomial \code{p} for point estimation is selected.
#' If the selected optimal degree is \code{p}, then piecewise polynomials of degree
#' \code{p+1} are used to conduct testing for nonparametric shape restrictions or parametric model specifications.
#' \emph{Note:} To implement the degree or smoothness selection, in addition to \code{pselect} or \code{sselect},
#' \code{nbins=#} must be specified.
#'@param sselect vector of numbers within which the number of smoothness constraints \code{s} for point estimation is selected.
#' If the selected optimal smoothness is \code{s}, then piecewise polynomials of \code{s+1} smoothness constraints
#' are used to conduct testing for nonparametric shape restrictions or parametric model specifications.
#' If not specified, for each value \code{p} supplied in the option \code{pselect}, only the piecewise polynomial
#' with the maximum smoothness is considered, i.e., \code{s=p}.
#'@param binspos position of binning knots. The default is \code{binspos="qs"}, which corresponds to quantile-spaced
#' binning (canonical binscatter). The other options are \code{"es"} for evenly-spaced binning, or
#' a vector for manual specification of the positions of inner knots (which must be within the range of
#' \code{x}).
#'@param binsmethod method for data-driven selection of the number of bins. The default is \code{binsmethod="dpi"},
#' which corresponds to the IMSE-optimal direct plug-in rule. The other option is: \code{"rot"}
#' for rule of thumb implementation.
#'@param nbinsrot initial number of bins value used to construct the DPI number of bins selector.
#' If not specified, the data-driven ROT selector is used instead.
#'@param randcut upper bound on a uniformly distributed variable used to draw a subsample for bins/degree/smoothness selection.
#' Observations for which \code{runif()<=#} are used. # must be between 0 and 1. By default, \code{max(5000, 0.01n)} observations
#' are used if the samples size \code{n>5000}.
#'@param nsims number of random draws for hypothesis testing. The default is
#' \code{nsims=500}, which corresponds to 500 draws from a standard Gaussian random vector of size
#' \code{[(p+1)*J - (J-1)*s]}. Setting at least \code{nsims=2000} is recommended to obtain the final results.
#'@param simsgrid number of evaluation points of an evenly-spaced grid within each bin used for evaluation of
#' the supremum (infimum or Lp metric) operation needed to construct hypothesis testing
#' procedures. The default is \code{simsgrid=20}, which corresponds to 20 evenly-spaced
#' evaluation points within each bin for approximating the supremum (infimum or Lp metric) operator.
#' Setting at least \code{simsgrid=50} is recommended to obtain the final results.
#'@param simsseed seed for simulation.
#'@param vce procedure to compute the variance-covariance matrix estimator. For least squares regression and generalized linear regression, the allowed options are the same as that for \code{\link{binsreg}} or \code{\link{binsqreg}}.
#' For quantile regression, the allowed options are the same as that for \code{\link{binsqreg}}.
#'@param cluster cluster ID. Used for compute cluster-robust standard errors.
#'@param asyvar if true, the standard error of the nonparametric component is computed and the uncertainty related to control
#' variables is omitted. Default is \code{asyvar=FALSE}, that is, the uncertainty related to control variables is taken into account.
#'@param dfcheck adjustments for minimum effective sample size checks, which take into account number of unique
#' values of \code{x} (i.e., number of mass points), number of clusters, and degrees of freedom of
#' the different stat models considered. The default is \code{dfcheck=c(20, 30)}.
#' See \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_Stata.pdf}{Cattaneo, Crump, Farrell and Feng (2024c)} for more details.
#'@param masspoints how mass points in \code{x} are handled. Available options:
#' \itemize{
#' \item \code{"on"} all mass point and degrees of freedom checks are implemented. Default.
#' \item \code{"noadjust"} mass point checks and the corresponding effective sample size adjustments are omitted.
#' \item \code{"nolocalcheck"} within-bin mass point and degrees of freedom checks are omitted.
#' \item \code{"off"} "noadjust" and "nolocalcheck" are set simultaneously.
#' \item \code{"veryfew"} forces the function to proceed as if \code{x} has only a few number of mass points (i.e., distinct values).
#' In other words, forces the function to proceed as if the mass point and degrees of freedom checks were failed.
#' }
#'@param weights an optional vector of weights to be used in the fitting process. Should be \code{NULL} or
#' a numeric vector. For more details, see \code{\link{lm}}.
#'@param subset optional rule specifying a subset of observations to be used.
#'@param numdist number of distinct values for selection. Used to speed up computation.
#'@param numclust number of clusters for selection. Used to speed up computation.
#'@param estmethodopt a list of optional arguments used by \code{\link[quantreg]{rq}} (for quantile regression) or \code{\link{glm}} (for fitting generalized linear models).
#'@param ... optional arguments to control bootstrapping if \code{estmethod="qreg"} and \code{vce="boot"}. See \code{\link[quantreg]{boot.rq}}.
#'@return \item{\code{testshapeL}}{Results for \code{testshapel}, including: \code{testvalL}, null boundary values;
#' \code{stat.shapeL}, test statistics; and \code{pval.shapeL}, p-value.}
#' \item{\code{testshapeR}}{Results for \code{testshaper}, including: \code{testvalR}, null boundary values;
#' \code{stat.shapeR}, test statistics; and \code{pval.shapeR}, p-value.}
#' \item{\code{testshape2}}{Results for \code{testshape2}, including: \code{testval2}, null boundary values;
#' \code{stat.shape2}, test statistics; and \code{pval.shape2}, p-value.}
#' \item{\code{testpoly}}{Results for \code{testmodelpoly}, including: \code{testpoly}, the degree of global polynomial;
#' \code{stat.poly}, test statistic; \code{pval.poly}, p-value.}
#' \item{\code{testmodel}}{Results for \code{testmodelparfit}, including: \code{stat.model}, test statistics;
#' \code{pval.model}, p-values.}
#' \item{\code{imse.var.rot}}{Variance constant in IMSE, ROT selection.}
#' \item{\code{imse.bsq.rot}}{Bias constant in IMSE, ROT selection.}
#' \item{\code{imse.var.dpi}}{Variance constant in IMSE, DPI selection.}
#' \item{\code{imse.bsq.dpi}}{Bias constant in IMSE, DPI selection.}
#' \item{\code{opt}}{ A list containing options passed to the function, as well as total sample size \code{n},
#' number of distinct values \code{Ndist} in \code{x}, number of clusters \code{Nclust}, and
#' number of bins \code{nbins}.}
#'
#'@author
#' Matias D. Cattaneo, Princeton University, Princeton, NJ. \email{cattaneo@princeton.edu}.
#'
#' Richard K. Crump, Federal Reserve Bank of New York, New York, NY. \email{richard.crump@ny.frb.org}.
#'
#' Max H. Farrell, UC Santa Barbara, Santa Barbara, CA. \email{mhfarrell@gmail.com}.
#'
#' Yingjie Feng (maintainer), Tsinghua University, Beijing, China. \email{fengyingjiepku@gmail.com}.
#'
#'@references
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2024a: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_AER.pdf}{On Binscatter}. American Economic Review 114(5): 1488-1514.
#'
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2024b: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_NonlinearBinscatter.pdf}{Nonlinear Binscatter Methods}. Working Paper.
#'
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2024c: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_Stata.pdf}{Binscatter Regressions}. Working Paper.
#'
#'@seealso \code{\link{binsreg}}, \code{\link{binsqreg}}, \code{\link{binsglm}}, \code{\link{binsregselect}}.
#'
#'@examples
#' x <- runif(500); y <- sin(x)+rnorm(500)
#' est <- binstest(y,x, testmodelpoly=1)
#' summary(est)
#'@export
binstest <- function(y, x, w=NULL, data=NULL, estmethod="reg", family=gaussian(),
quantile=NULL, deriv=0, at=NULL, nolink=F,
testmodel=NULL, testmodelparfit=NULL, testmodelpoly=NULL,
testshape=NULL, testshapel=NULL,
testshaper=NULL, testshape2=NULL, lp=Inf,
bins=NULL, nbins=NULL,
pselect=NULL, sselect=NULL,
binspos="qs", binsmethod="dpi", nbinsrot=NULL, randcut=NULL,
nsims=500, simsgrid=20, simsseed=NULL,
vce=NULL, cluster=NULL, asyvar=F,
dfcheck=c(20,30), masspoints="on", weights=NULL, subset=NULL,
numdist=NULL, numclust=NULL, estmethodopt=NULL, ...) {
# param for internal use
qrot <- 2
####################
### prepare data ###
####################
xname <- deparse(substitute(x), width.cutoff = 500L)
yname <- deparse(substitute(y), width.cutoff = 500L)
weightsname <- deparse(substitute(weights), width.cutoff = 500L)
subsetname <- deparse(substitute(subset), width.cutoff = 500L)
clustername <- deparse(substitute(cluster), width.cutoff = 500L)
# extract y, x, w, weights, subset, if needed (w as a matrix, others as vectors)
# generate design matrix for covariates W
w.factor <- NULL
if (is.null(data)) {
if (is.data.frame(y)) y <- y[,1]
if (is.data.frame(x)) x <- x[,1]
if (!is.null(w)) {
if (is.vector(w)) {
w <- as.matrix(w)
if (is.character(w)) {
w <- model.matrix(~w)[,-1,drop=F]
w.factor <- rep(T, ncol(w))
}
} else if (is.formula(w)) {
w.model <- binsreg.model.mat(w)
w <- w.model$design
w.factor <- w.model$factor.colnum
} else if (is.data.frame(w)) {
w <- as.matrix(w)
}
}
} else {
if (xname %in% names(data)) x <- data[,xname]
if (yname %in% names(data)) y <- data[,yname]
if (weightsname != "NULL") if (weightsname %in% names(data)) {
weights <- data[,weightsname]
}
if (subsetname != "NULL") if (subsetname %in% names(data)) {
subset <- data[,subsetname]
}
if (clustername != "NULL") if (clustername %in% names(data)) {
cluster <- data[,clustername]
}
if (deparse(substitute(w), width.cutoff = 500L)[1]!="NULL") {
if (is.formula(w)) {
if (is.data.frame(at)) {
if (ncol(data)!=ncol(at)) {
data[setdiff(names(at), names(data))] <- NA
at[setdiff(names(data), names(at))] <- NA
}
data <- rbind(data, at)
}
w.model <- binsreg.model.mat(w, data=data)
w <- w.model$design
w.factor <- w.model$factor.colnum
if (is.data.frame(at)) {
eval.w <- w[nrow(w),]
w <- w[-nrow(w),, drop=F]
}
} else {
if (deparse(substitute(w), width.cutoff = 500L) %in% names(data)) w <- data[,deparse(substitute(w), width.cutoff = 500L)]
w <- as.matrix(w)
if (is.character(w)) {
w <- model.matrix(~w)[,-1,drop=F]
w.factor <- rep(T, ncol(w))
}
}
}
}
if (!is.null(w)) nwvar <- ncol(w)
else nwvar <- 0
# extract subset
if (!is.null(subset)) {
y <- y[subset]
x <- x[subset]
w <- w[subset, , drop = F]
if (!is.null(cluster)) cluster <- cluster[subset]
if (!is.null(weights)) weights <- weights[subset]
}
na.ok <- complete.cases(x) & complete.cases(y)
if (!is.null(w)) na.ok <- na.ok & complete.cases(w)
if (!is.null(weights)) na.ok <- na.ok & complete.cases(weights)
if (!is.null(cluster)) na.ok <- na.ok & complete.cases(cluster)
y <- y[na.ok]
x <- x[na.ok]
w <- w[na.ok, , drop = F]
weights <- weights[na.ok]
cluster <- cluster[na.ok]
xmin <- min(x); xmax <- max(x)
# evaluation point of w
if (is.null(at)) at <- "mean"
# change method name if needed
if (is.character(family))
family <- get(family, mode = "function", envir = parent.frame())
if (is.function(family))
family <- family()
if (is.null(family$family)) {
print(family)
stop("'family' not recognized")
}
if (!is.null(family)) if (family$family!="gaussian" | family$link!="identity") {
estmethod <- "glm"
}
# extract family obj
if (estmethod=="glm") {
familyname <- family$family
linkinv <- family$linkinv
linkinv.1 <- family$mu.eta # 1st deriv of inverse link
# 2nd derivative
if (family$link=="logit") {
linkinv.2 <- function(z) linkinv.1(z)*(1-2*linkinv(z))
} else if (family$link=="probit") {
linkinv.2 <- function(z) (-z)*linkinv.1(z)
} else if (family$link=="identity") {
linkinv.2 <- function(z) 0
} else if (family$link=="log") {
linkinv.2 <- function(z) linkinv(z)
} else if (family$link=="inverse") {
linkinv.2 <- function(z) 2*z^(-3)
} else if (family$link=="1/mu^2") {
linkinv.2 <- function(z) 0.75*z^(-2.5)
} else if (family$link=="sqrt") {
linkinv.2 <- function(z) 2
} else {
print("The specified link not allowed for computing polynomial-based CIs.")
stop()
}
}
# define the estimation method
if (estmethod=="reg") {
estmethod.name <- "least squares regression"
if (is.null(vce)) vce <- "HC1"
vce.select <- vce
} else if (estmethod=="qreg") {
estmethod.name <- "quantile regression"
if (is.null(vce)) vce <- "nid"
if (vce=="iid") {
vce.select <- "const"
} else {
vce.select <- "HC1"
}
if (is.null(quantile)) quantile <- 0.5
} else if (estmethod=="glm") {
estmethod.name <- "generalized linear model"
if (is.null(vce)) vce <- "HC1"
vce.select <- vce
} else {
print("estmethod incorrectly specified.")
stop()
}
# analyze bins- and degree-related options
tshaON <- tmodON <- F
if (!is.null(testshapel)|!is.null(testshaper)|!is.null(testshape2)) tshaON <- T
if (!is.null(testmodelparfit)|!is.null(testmodelpoly)) tmodON <- T
if (is.logical(testshape)) if (!testshape) testshape <- NULL
if (is.logical(testmodel)) if (!testmodel) testmodel <- NULL
if (is.logical(testshape) & (!tshaON | !is.character(binspos))) testshape <- NULL
if (is.logical(testmodel) & (!tmodON | !is.character(binspos))) testmodel <- NULL
if (!is.character(binspos)) {
if (!is.null(nbins)|!is.null(pselect)|!is.null(sselect)) {
print("nbins, pselect or sselect incorrectly specified.")
stop()
}
}
# 4 cases: select J, select p, user specify both, or an error
if (!is.null(bins)) {
bins.p <- bins[1]
if (length(bins)==1) {
bins.s <- bins.p
} else if (length(bins)==2) {
bins.s <- bins[2]
} else {
print("bins not correctly specified.")
stop()
}
plist <- bins.p; slist <- bins.s
if (is.numeric(nbins)) if (length(nbins)==1) if (nbins!=0) {
print("bins or nbins is incorrectly specified.")
stop()
}
} else {
plist <- pselect; slist <- sselect
bins.p <- bins.s <- NULL
}
len_nbins <- length(nbins); len_p <- length(plist); len_s <- length(slist)
if (len_p==1 & len_s==0) {
slist <- plist; len_s <- 1
}
if (len_p==0 & len_s==1) {
plist <- slist; len_p <- 1
}
# 1st case: select J
selection <- ""
if (is.character(binspos)) {
if (is.logical(nbins)) if (nbins) selection <- "J"
if (len_nbins==1) if (nbins==0) selection <- "J"
if (is.null(nbins)|len_nbins>1|!is.null(bins)) selection <- "J" # turn on J selection
}
if (selection=="J") {
if (len_p>1|len_s>1) {
if (is.null(nbins)) {
print("nbins must be specified for degree/smoothness selection.")
stop()
} else {
print("Only one p and one s are allowed to select # of bins.")
stop()
}
}
if (is.null(plist)) plist <- deriv
if (is.null(slist)) slist <- plist
if (is.null(bins)) {
bins.p <- plist; bins.s <- slist
}
len_p <- len_s <- 1
if (is.null(testshape)) testshape <- c(plist+1, slist+1)
if (is.logical(testshape)) if (testshape) {
testshape <- c(plist+1, slist+1)
}
if (is.null(testmodel)) testmodel <- c(plist+1, slist+1)
if (is.logical(testmodel)) if (testmodel) {
testmodel <- c(plist+1, slist+1)
}
}
# 2nd case: select p (at least for one object)
pselectOK <- F
if (selection!="J" & (is.logical(testshape)|is.null(testshape)) & tshaON) {
pselectOK <- T
}
if (selection!="J" & (is.logical(testmodel)|is.null(testmodel)) & tmodON) {
pselectOK <- T
}
if (pselectOK & (len_p>1|len_s>1) &len_nbins==1) {
selection <- "P"
}
# 3rd case: user specified
if ((len_p<=1 & len_s<=1) & selection!="J") {
selection <- "U"
if (is.null(testshape)) {
if (len_p==1 & len_s==1) testshape <- c(plist+1, slist+1)
else testshape <- c(deriv+1, deriv+1)
}
if (is.null(testmodel)) {
if (len_p==1 & len_s==1) testmodel <- c(plist+1, slist+1)
else testmodel <- c(deriv+1, deriv+1)
}
}
if (selection=="") {
print("Degree, smoothness, or # of bins not correctly specified")
stop()
}
##################################error checking
exit <- 0
if (!is.character(binspos)) {
if (min(binspos)<=xmin|max(binspos)>=xmax) {
print("Knots out of allowed range")
exit <- 1
}
} else {
if (binspos!="es"&binspos!="qs") {
print("binspos incorrectly specified.")
exit <- 1
}
}
if (length(bins)==2) if (bins[1]<bins[2]) {
print("p<s not allowed.")
exit <- 1
}
if (lp<1) {
print("lp has to be no less than 1.")
exit <- 1
}
if (!is.null(testshapel) | !is.null(testshaper)) if (lp!=Inf) {
print("Sup-norm (lp=Inf) has to be used for testing one-sided restrictions.")
exit <- 1
}
if (length(testshape)==2) if (testshape[1]<testshape[2]) {
print("p<s not allowed.")
exit <- 1
}
if (length(testmodel)==2) if (testmodel[1]<testmodel[2]) {
print("p<s not allowed.")
exit <- 1
}
if (length(testshape)>=1) if(!is.logical(testshape)) if (testshape[1]<deriv) {
print("p<deriv not allowed.")
exit <- 1
}
if (length(testmodel)>=1) if(!is.logical(testmodel)) if (testmodel[1]<deriv) {
print("p<deriv not allowed.")
exit <- 1
}
if (binsmethod!="dpi" & binsmethod!="rot") {
print("bin selection method incorrectly specified.")
exit <- 1
}
if (masspoints == "veryfew") {
print("veryfew not allowed for testing.")
exit <- 1
}
if ((length(testshape)>=1)&(length(bins)>=1)) if (!is.logical(testshape)) if (testshape[1]<=bins[1]) {
warning("p for testing <= p for bin selection not suggested.")
}
if ((length(testmodel)>=1)&(length(bins)>=1)) if (!is.logical(testmodel)) if (testmodel[1]<=bins[1]) {
warning("p for testing <= p for bin selection not suggested.")
}
if (exit>0) stop()
if (nsims<2000|simsgrid<50) {
print("Note: Setting at least nsims=2000 and simsgrid=50 is recommended to obtain the final results.")
}
#################################################
localcheck <- massadj <- T
if (masspoints=="on") {
localcheck <- T; massadj <- T
} else if (masspoints=="off") {
localcheck <- F; massadj <- F
} else if (masspoints=="noadjust") {
localcheck <- T; massadj <- F
} else if (masspoints=="nolocalcheck") {
localcheck <- F; massadj <- T
}
# effective size
eN <- N <- length(x)
Ndist <- NA
if (massadj) {
if (!is.null(numdist)) {
Ndist <- numdist
} else {
Ndist <- length(unique(x))
}
eN <- min(eN, Ndist)
}
Nclust <- NA
if (!is.null(cluster)) {
if (!is.null(numclust)) {
Nclust <- numclust
} else {
Nclust <- length(unique(cluster))
}
eN <- min(eN, Nclust)
}
# prepare params
tsha.p <- testshape[1]; tsha.s <- testshape[2]
if (is.logical(tsha.p)) tsha.p <- NULL
if (!is.null(tsha.s)) if (is.na(tsha.s)) tsha.s <- tsha.p
tmod.p <- testmodel[1]; tmod.s <- testmodel[2]
if (is.logical(tmod.p)) tmod.p <- NULL
if (!is.null(tmod.s)) if (is.na(tmod.s)) tmod.s <- tmod.p
if (selection=="J") {
if (!is.null(tsha.p)) if (tsha.p<=plist) {
tsha.p <- bins.p+1; tsha.s <- tsha.p
warning("Degree for testshape has been changed. It must be greater than the degree for bin selection.")
}
if (!is.null(tmod.p)) if (tmod.p<=plist) {
tmod.p <- bins.p+1; tmod.s <- tmod.p
warning("Degree for testmodel has been changed. It must be greater than the degree for bin selection.")
}
}
if (selection=="U") {
warning("Testing procedures are valid when nbins is much larger than the IMSE-optimal choice. Compare your choice with the IMSE-optimal one obtained by binsregselect().")
}
nL <- length(testshapel); nR <- length(testshaper); nT <- length(testshape2)
ntestshape <- nL + nR + nT
if (selection=="U") {
selectmethod <- "User-specified"
} else {
if (binsmethod=="dpi") {
selectmethod <- "IMSE direct plug-in"
} else {
selectmethod <- "IMSE rule-of-thumb"
}
if (selection=="J") selectmethod <- paste(selectmethod, "(select # of bins)")
if (selection=="P") selectmethod <- paste(selectmethod, "(select degree and smoothness)")
}
knot <- NULL
if (!is.character(binspos)) {
nbins <- length(binspos)+1
knot <- c(xmin, binspos, xmax)
position <- "User-specified"
es <- F
} else {
if (binspos == "es") {
es <- T
position <- "Evenly-spaced"
} else {
es <- F
position <- "Quantile-spaced"
}
}
#### bin selection if needed ########
imse.v.rot <- imse.v.dpi <- imse.b.rot <- imse.b.dpi <- NA
if (selection!="U") {
# check if rot can be implemented
if (is.null(bins.p)) binspcheck <- 6
else binspcheck <- bins.p
if (is.null(nbinsrot)) {
if (eN <= dfcheck[1]+binspcheck+1+qrot) {
print("Too small effective sample size for bin selection.")
stop()
}
}
if (is.na(Ndist)) {
Ndist.sel <- NULL
} else {
Ndist.sel <- Ndist
}
if (is.na(Nclust)) {
Nclust.sel <- NULL
} else {
Nclust.sel <- Nclust
}
randcut1k <- randcut
if (is.null(randcut) & N>5000) {
randcut1k <- max(5000/N, 0.01)
warning("To speed up computation, bin/degree selection uses a subsample of roughly max(5000, 0.01n) observations if the sample size n>5000. To use the full sample, set randcut=1.")
}
if (selection=="J") {
binselect <- binsregselect(y, x, w, deriv=deriv,
bins=c(bins.p,bins.s), binspos=binspos, nbins=nbins,
binsmethod=binsmethod, nbinsrot=nbinsrot,
vce=vce.select, cluster=cluster, randcut=randcut1k,
dfcheck=dfcheck, masspoints=masspoints, weights=weights,
numdist=Ndist.sel, numclust=Nclust.sel)
if (is.na(binselect$nbinsrot.regul)) {
print("Bin selection fails.")
stop()
}
if (binsmethod == "rot") {
nbins <- binselect$nbinsrot.regul
imse.v.rot <- binselect$imse.var.rot
imse.b.rot <- binselect$imse.bsq.rot
} else if (binsmethod == "dpi") {
nbins <- binselect$nbinsdpi
imse.v.dpi <- binselect$imse.var.dpi
imse.b.dpi <- binselect$imse.bsq.dpi
if (is.na(nbins)) {
warning("DPI selection fails. ROT choice used.")
nbins <- binselect$nbinsrot.regul
imse.v.rot <- binselect$imse.var.rot
imse.b.rot <- binselect$imse.bsq.rot
}
}
} else if (selection=="P") {
binselect <- binsregselect(y, x, w, deriv=deriv,
binspos=binspos, nbins=nbins,
pselect=plist, sselect=slist,
binsmethod=binsmethod, nbinsrot=nbinsrot,
vce=vce.select, cluster=cluster, randcut=randcut1k,
dfcheck=dfcheck, masspoints=masspoints, weights=weights,
numdist=Ndist.sel, numclust=Nclust.sel)
if (is.na(binselect$prot.regul)) {
print("Bin selection fails.")
stop()
}
if (binsmethod == "rot") {
bins.p <- binselect$prot.regul
bins.s <- binselect$srot.regul
imse.v.rot <- binselect$imse.var.rot
imse.b.rot <- binselect$imse.bsq.rot
} else if (binsmethod == "dpi") {
bins.p <- binselect$pdpi
bins.s <- binselect$sdpi
imse.v.dpi <- binselect$imse.var.dpi
imse.b.dpi <- binselect$imse.bsq.dpi
if (is.na(bins.p)) {
warning("DPI selection fails. ROT choice used.")
bins.p <- binselect$prot.regul
bins.s <- binselect$srot.regul
imse.v.rot <- binselect$imse.var.rot
imse.b.rot <- binselect$imse.bsq.rot
}
}
if (is.logical(testshape)|is.null(testshape)) {
tsha.p <- bins.p+1; tsha.s <- bins.s+1
} else {
if (tsha.p<=bins.p) {
tsha.p <- bins.p+1; tsha.s <- tsha.p
warning("Degree for testshape has been changed. It must be greater than the IMSE-optimal degree.")
}
}
if (is.logical(testmodel)|is.null(testmodel)) {
tmod.p <- bins.p+1; tmod.s <- bins.s+1
} else {
if (tmod.p<=bins.p) {
tmod.p <- bins.p+1; tmod.s <- tmod.p
warning("Degree for testmodel has been changed. It must be greater than the IMSE-optimal degree.")
}
}
}
}
# check eff size for testing
tsha.fewobs <- F
if (ntestshape!=0) if ((nbins-1)*(tsha.p-tsha.s+1)+tsha.p+1+dfcheck[2]>=eN) {
tsha.fewobs <- T
warning("Too small eff. sample size for testing shape.")
}
tmod.fewobs <- F
if (!is.null(testmodelparfit)|!is.null(testmodelpoly)) if ((nbins-1)*(tmod.p-tmod.s+1)+tmod.p+1+dfcheck[2]>=eN) {
tmod.fewobs <- T
warning("Too small eff. sample size for testing models.")
}
# Generate knot
if (is.null(knot)) {
if (es) {
knot <- genKnot.es(xmin, xmax, nbins)
} else {
knot <- genKnot.qs(x, nbins)
}
}
# check local mass points
knot <- c(knot[1], unique(knot[-1]))
if (nbins!=length(knot)-1) {
warning("Repeated knots. Some bins dropped.")
nbins <- length(knot)-1
}
if (localcheck) {
uniqmin <- binsreg.checklocalmass(x, nbins, es, knot=knot) # mimic STATA
if (ntestshape != 0) {
if (uniqmin < tsha.p+1) {
tsha.fewobs <- T
warning("Some bins have too few distinct values of x for testing shape.")
}
}
if (!is.null(testmodelparfit)|!is.null(testmodelpoly)) {
if (uniqmin < tmod.p+1) {
tmod.fewobs <- T
warning("Some bins have too few distinct values of x for testing models.")
}
}
}
# adjust w variables
if (!is.null(w)) {
if (is.character(at)) {
if (at=="mean") {
eval.w <- colWeightedMeans(x=w, w=weights)
if (!is.null(w.factor)) eval.w[w.factor] <- 0
} else if (at=="median") {
eval.w <- colWeightedMedians(x=w, w=weights)
if (!is.null(w.factor)) eval.w[w.factor] <- 0
} else if (at=="zero") {
eval.w <- rep(0, nwvar)
}
} else if (is.vector(at)) {
eval.w <- at
} else if (is.data.frame(at)) {
eval.w <- eval.w
}
} else {
eval.w <- NULL
}
# seed
if (!is.null(simsseed)) set.seed(simsseed)
# prepare grid for uniform inference
x.grid <- binsreg.grid(knot, simsgrid)$eval
#####################################
####### Shape restriction test ######
#####################################
stat.shapeL <- pval.shapeL <- NA
stat.shapeR <- pval.shapeR <- NA
stat.shape2 <- pval.shape2 <- NA
if (nL>0) {
stat.shapeL <- matrix(NA, nL, 2)
pval.shapeL <- c()
}
if (nR>0) {
stat.shapeR <- matrix(NA, nR, 2)
pval.shapeR <- c()
}
if (nT>0) {
stat.shape2 <- matrix(NA, nT, 2)
pval.shape2 <- c()
}
if (ntestshape != 0 & !tsha.fewobs) {
B <- binsreg.spdes(eval=x, p=tsha.p, s=tsha.s, knot=knot, deriv=0)
k <- ncol(B)
P <- cbind(B, w)
if (estmethod=="reg") {
model <- lm(y ~ P-1, weights=weights)
} else if (estmethod=="qreg") {
model <- do.call(rq, c(list(formula=y ~ P-1, tau=quantile, weights=weights), estmethodopt))
} else if (estmethod=="glm") {
model <- do.call(glm, c(list(formula=y ~ P-1, family=family, weights=weights), estmethodopt))
}
beta <- model$coeff[1:k]
basis.sha <- binsreg.spdes(eval=x.grid, p=tsha.p, s=tsha.s, knot=knot, deriv=deriv)
if (estmethod=="glm" & (!nolink)) {
pred.sha <- binsreg.pred(X=basis.sha, model=model, type="all",
vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
basis.0 <- binsreg.spdes(eval=x.grid, p=tsha.p, s=tsha.s, knot=knot, deriv=0)
fit.0 <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
pred.sha.0 <- linkinv.1(fit.0)
if (asyvar | deriv==0) {
pred.sha$se <- pred.sha.0 * pred.sha$se
if (deriv == 0) pred.sha$fit <- linkinv(pred.sha$fit)
if (deriv == 1) pred.sha$fit <- pred.sha.0 * pred.sha$fit
} else {
basis.sha.0 <- basis.0
basis.sha.1 <- basis.sha
if (!is.null(eval.w)) {
basis.sha.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
basis.sha.1 <- cbind(basis.sha.1, outer(rep(1, nrow(basis.sha.1)), rep(0, nwvar)))
}
basis.all <- linkinv.2(fit.0)*pred.sha$fit*basis.sha.0 + pred.sha.0*basis.sha.1
pred.sha$fit <- pred.sha.0 * pred.sha$fit
pred.sha$se <- binsreg.pred(basis.all, model=model, type="se",
vce=vce, cluster=cluster, avar=T)$se
}
} else {
if (estmethod=="qreg") {
pred.sha <- binsreg.pred(basis.sha, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w,
is.qreg=TRUE, avar=asyvar, ...)
} else {
pred.sha <- binsreg.pred(basis.sha, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
}
}
fit.sha <- pred.sha$fit
se.sha <- pred.sha$se
pos <- !is.na(beta)
k.new <- sum(pos)
if (estmethod=="qreg") vcv.sha <- binsreg.vcov(model, type=vce, cluster=cluster, is.qreg=TRUE, ...)[1:k.new, 1:k.new]
else vcv.sha <- binsreg.vcov(model, type=vce, cluster=cluster)[1:k.new, 1:k.new]
for (j in 1:ntestshape) {
if (j <= nL) {
stat.shapeL[j,2] <- 1
stat.shapeL[j,1] <- max((fit.sha - testshapel[j]) / se.sha)
} else if (j <= nL+nR) {
stat.shapeR[j-nL,2] <- 2
stat.shapeR[j-nL,1] <- min((fit.sha - testshaper[j-nL]) / se.sha)
} else {
stat.shape2[j-nL-nR,2] <- 3
if (is.infinite(lp)) {
stat.shape2[j-nL-nR,1] <- max(abs((fit.sha - testshape2[j-nL-nR]) / se.sha))
} else {
stat.shape2[j-nL-nR,1] <- mean(abs((fit.sha - testshape2[j-nL-nR]) / se.sha)^lp)^(1/lp)
}
}
}
stat.shape <- NULL
if (nL > 0) {
stat.shape <- rbind(stat.shape, stat.shapeL)
}
if (nR > 0) {
stat.shape <- rbind(stat.shape, stat.shapeR)
}
if (nT > 0) {
stat.shape <- rbind(stat.shape, stat.shape2)
}
# Compute p-val
Sigma.root <- lssqrtm(vcv.sha)
num <- basis.sha[,pos,drop=F] %*% Sigma.root
denom.sha <- sqrt(rowSums((basis.sha[, pos, drop=F] %*% vcv.sha) * basis.sha[, pos, drop=F]))
pval.shape <- binsreg.pval(num, denom.sha, nsims, tstat=stat.shape, side=NULL, lp=lp)$pval
if (nL!=0) {
stat.shapeL <- stat.shapeL[,1]
pval.shapeL <- pval.shape[1:nL]
}
if (nR!=0) {
stat.shapeR <- stat.shapeR[,1]
pval.shapeR <- pval.shape[(nL+1):(nL+nR)]
}
if (nT!=0) {
stat.shape2 <- stat.shape2[,1]
pval.shape2 <- pval.shape[(nL+nR+1):ntestshape]
}
}
############################
#### Specification Test ####
############################
stat.poly <- pval.poly <- NA
stat.mod <- pval.mod <- NA
if (!is.null(testmodelpoly)) {
stat.poly <- matrix(NA, 1,2)
pval.poly <- c()
}
if (!is.null(testmodelparfit)) {
stat.mod <- matrix(NA,ncol(testmodelparfit)-1,2)
pval.mod <- c()
}
if ((!is.null(testmodelparfit)|!is.null(testmodelpoly))&!tmod.fewobs) {
if (tmod.p==tsha.p & tmod.s==tsha.s & exists("vcv.sha")) {
exist.mod <- T
#beta.mod <- beta.sha
vcv.mod <- vcv.sha
fit.mod <- fit.sha
se.mod <- se.sha
basis.mod <- basis.sha
denom.mod <- denom.sha
pos <- pos
k.new <- k.new
model <- model
} else {
exist.mod <- F
B <- binsreg.spdes(eval=x, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
k <- ncol(B)
P <- cbind(B, w)
if (estmethod=="reg") {
model <- lm(y ~ P-1, weights=weights)
} else if (estmethod=="qreg") {
model <- do.call(rq, c(list(formula=y ~ P-1, tau=quantile, weights=weights), estmethodopt))
} else {
model <- do.call(glm, c(list(formula=y ~ P-1, family=family, weights=weights), estmethodopt))
}
beta <- model$coeff[1:k]
pos <- !is.na(beta)
# beta.mod <- beta[pos]
k.new <- sum(pos)
if (estmethod=="qreg") vcv.mod <- binsreg.vcov(model, type=vce, cluster=cluster, is.qreg = TRUE, ...)[1:k.new, 1:k.new]
else vcv.mod <- binsreg.vcov(model, type=vce, cluster=cluster)[1:k.new, 1:k.new]
}
######################
# Test poly reg
if (!is.null(testmodelpoly)) {
if (!exist.mod) {
basis.mod <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=deriv)
if (estmethod=="glm" & (!nolink)) {
pred.mod <- binsreg.pred(X=basis.mod, model=model, type="all",
vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
basis.0 <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
fit.0 <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
pred.mod.0 <- linkinv.1(fit.0)
if (asyvar | deriv==0) {
pred.mod$se <- pred.mod.0 * pred.mod$se
if (deriv == 0) pred.mod$fit <- linkinv(pred.mod$fit)
if (deriv == 1) pred.mod$fit <- pred.mod.0 * pred.mod$fit
} else {
basis.mod.0 <- basis.0
basis.mod.1 <- basis.mod
if (!is.null(eval.w)) {
basis.mod.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
basis.mod.1 <- cbind(basis.mod.1, outer(rep(1, nrow(basis.mod.1)), rep(0, nwvar)))
}
basis.all <- linkinv.2(fit.0)*pred.mod$fit*basis.mod.0 + pred.mod.0*basis.mod.1
pred.mod$fit <- pred.mod.0 * pred.mod$fit
pred.mod$se <- binsreg.pred(basis.all, model=model, type="se",
vce=vce, cluster=cluster, avar=T)$se
}
} else {
if (estmethod=="qreg") {
pred.mod <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w,
is.qreg=TRUE, avar=asyvar, ...)
} else {
pred.mod <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
}
}
fit.mod <- pred.mod$fit
se.mod <- pred.mod$se
denom.mod <- sqrt(rowSums((basis.mod[, pos, drop=F] %*% vcv.mod) * basis.mod[, pos, drop=F]))
}
# Run a poly reg
x.p <- matrix(NA, N, testmodelpoly+1)
for (j in 1:(testmodelpoly+1)) x.p[,j] <- x^(j-1)
P.poly <- cbind(x.p, w)
if (estmethod=="reg") {
model.poly <- lm(y~P.poly-1, weights=weights)
} else if (estmethod=="qreg") {
model.poly <- do.call(rq, c(list(formula=y~P.poly-1, tau=quantile, weights=weights), estmethodopt))
} else if (estmethod=="glm") {
model.poly <- do.call(glm, c(list(formula=y~P.poly-1, family=family, weights=weights), estmethodopt))
}
beta.poly <- model.poly$coefficients
poly.fit <- 0
for (j in deriv:testmodelpoly) {
poly.fit <- poly.fit + x.grid^(j-deriv)*beta.poly[j+1]*factorial(j)/factorial(j-deriv)
}
if (!is.null(eval.w) & deriv==0) poly.fit <- poly.fit + sum(eval.w * beta.poly[-(1:(testmodelpoly+1))])
stat.poly[1,2] <- 3
if (is.infinite(lp)) stat.poly[1,1] <- max(abs((fit.mod - poly.fit) / se.mod))
else stat.poly[1,1] <- mean(abs((fit.mod - poly.fit) / se.mod)^lp)^(1/lp)
Sigma.root <- lssqrtm(vcv.mod)
num <- basis.mod[,pos,drop=F] %*% Sigma.root
pval.poly <- binsreg.pval(num, denom.mod, nsims, tstat=stat.poly, side=NULL, lp=lp)$pval
stat.poly <- stat.poly[,1]
}
##################################
# Test model in another data frame
if (!is.null(testmodelparfit)) {
x.grid <- testmodelparfit[,1]
basis.mod <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=deriv)
if (estmethod=="glm" & (!nolink)) {
pred.mod <- binsreg.pred(X=basis.mod, model=model, type="all",
vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
basis.0 <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
fit.0 <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
pred.mod.0 <- linkinv.1(fit.0)
if (asyvar | deriv==0) {
pred.mod$se <- pred.mod.0 * pred.mod$se
if (deriv == 0) pred.mod$fit <- linkinv(pred.mod$fit)
if (deriv == 1) pred.mod$fit <- pred.mod.0 * pred.mod$fit
} else {
basis.mod.0 <- basis.0
basis.mod.1 <- basis.mod
if (!is.null(eval.w)) {
basis.mod.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
basis.mod.1 <- cbind(basis.mod.1, outer(rep(1, nrow(basis.mod.1)), rep(0, nwvar)))
}
basis.all <- linkinv.2(fit.0)*pred.mod$fit*basis.mod.0 + pred.mod.0*basis.mod.1
pred.mod$fit <- pred.mod.0 * pred.mod$fit
pred.mod$se <- binsreg.pred(basis.all, model=model, type="se",
vce=vce, cluster=cluster, avar=T)$se
}
} else {
if (estmethod=="qreg") {
pred.mod <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w,
is.qreg=TRUE, avar=asyvar, ...)
} else {
pred.mod <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
}
}
fit.mod <- pred.mod$fit
se.mod <- pred.mod$se
denom.mod <- sqrt(rowSums((basis.mod[, pos, drop=F] %*% vcv.mod) * basis.mod[, pos, drop=F]))
for (j in 2:ncol(testmodelparfit)) {
stat.mod[j-1,2] <- 3
if (is.infinite(lp)) stat.mod[j-1,1] <- max(abs((fit.mod - testmodelparfit[,j]) / se.mod))
else stat.mod[j-1,1] <- mean(abs((fit.mod - testmodelparfit[,j]) / se.mod)^lp)^(1/lp)
}
Sigma.root <- lssqrtm(vcv.mod)
num <- basis.mod[,pos,drop=F] %*% Sigma.root
pval.mod <- binsreg.pval(num, denom.mod, nsims, tstat=stat.mod, side=NULL, lp=lp)$pval
stat.mod <- stat.mod[,1]
}
}
######################
#######output#########
######################
out <- list(testshapeL=list(testvalL=testshapel, stat.shapeL=stat.shapeL, pval.shapeL=pval.shapeL),
testshapeR=list(testvalR=testshaper, stat.shapeR=stat.shapeR, pval.shapeR=pval.shapeR),
testshape2=list(testval2=testshape2, stat.shape2=stat.shape2, pval.shape2=pval.shape2),
testpoly=list(testpoly=testmodelpoly, stat.poly=stat.poly, pval.poly=pval.poly),
testmodel=list(stat.model=stat.mod, pval.model=pval.mod),
imse.var.dpi=imse.v.dpi, imse.bsq.dpi=imse.b.dpi,
imse.var.rot=imse.v.rot, imse.bsq.rot=imse.b.rot,
opt = list(bins.p=bins.p, bins.s=bins.s, deriv=deriv,
testshape=c(tsha.p, tsha.s), testmodel=c(tmod.p, tmod.s),
binspos=position, binsmethod=selectmethod,
n=N, Ndist=Ndist, Nclust=Nclust,
nbins=nbins, knot=knot, lp=lp,
estmethod=estmethod.name, quantile=quantile, family=family))
out$call <- match.call()
class(out) <- "CCFFbinstest"
return(out)
}
##########################################################################
#' Internal function.
#'
#' @param x Class \code{CCFFbinstest} objects.
#'
#' @keywords internal
#' @export
#'
print.CCFFbinstest <- function(x, ...) {
cat("Call: binstest\n\n")
cat(paste("Sample size (n) = ", x$opt$n, "\n", sep=""))
cat(paste("# of distinct values (Ndist) = ", x$opt$Ndist, "\n", sep=""))
cat(paste("# of clusters (Nclust) = ", x$opt$Nclust, "\n", sep=""))
cat(paste("Estimation Method (estmethod) = ", x$opt$estmethod, "\n", sep=""))
cat(paste("Derivative (deriv) = ", x$opt$deriv, "\n", sep=""))
cat(paste("Bin/Degree selection:", "\n"))
cat(paste(" Method (binsmethod) = ", x$opt$binsmethod,"\n", sep=""))
cat(paste(" Placement (binspos) = ", x$opt$binspos, "\n", sep=""))
cat(paste(" degree (p) = ", x$opt$bins.p, "\n", sep=""))
cat(paste(" smooth (s) = ", x$opt$bins.s, "\n", sep=""))
cat(paste(" # of bins (nbins) = ", x$opt$nbins, "\n", sep=""))
cat("\n")
}
################################################################################
#' Internal function.
#'
#' @param object Class \code{CCFFbinstest} objects.
#'
#' @keywords internal
#' @export
summary.CCFFbinstest <- function(object, ...) {
x <- object
args <- list(...)
cat("Call: binstest\n\n")
cat(paste("Sample size (n) = ", x$opt$n, "\n", sep=""))
cat(paste("# of distinct values (Ndist) = ", x$opt$Ndist, "\n", sep=""))
cat(paste("# of clusters (Nclust) = ", x$opt$Nclust, "\n", sep=""))
cat(paste("Estimation Method (estmethod) = ", x$opt$estmethod, "\n", sep=""))
cat(paste("Derivative (deriv) = ", x$opt$deriv, "\n", sep=""))
cat(paste("Bin/Degree selection:", "\n"))
cat(paste(" Method (binsmethod) = ", x$opt$binsmethod,"\n", sep=""))
cat(paste(" Placement (binspos) = ", x$opt$binspos, "\n", sep=""))
cat(paste(" degree (p) = ", x$opt$bins.p, "\n", sep=""))
cat(paste(" smooth (s) = ", x$opt$bins.s, "\n", sep=""))
cat(paste(" # of bins (nbins) = ", x$opt$nbins, "\n", sep=""))
cat("\n")
if (!is.null(x$testshapeL$testvalL)|!is.null(x$testshapeR$testvalR)|!is.null(x$testshape2$testval2)) {
cat(paste("Shape Restriction Tests:")); cat("\n")
cat(paste("degree (p) = ", x$opt$testshape[1], sep="")); cat(paste(" "))
cat(paste("smooth (s) = ", x$opt$testshape[2], sep="")); cat("\n")
if (!is.null(x$testshapeL$testvalL)) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: sup mu <=", width = 15, justify = "left"))
cat(format("sup T", width = 12, justify = "right"))
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
for (i in 1:length(x$testshapeL$testvalL)) {
cat(format(paste(x$testshapeL$testvalL[i]), width = 15, justify = "centre"))
cat(format(sprintf("%3.3f", x$testshapeL$stat.shapeL[i]), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testshapeL$pval.shapeL[i]), width = 12, justify = "right"))
cat("\n")
}
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
cat("\n")
}
if (!is.null(x$testshapeR$testvalR)) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: inf mu >=", width = 15, justify = "left"))
cat(format("inf T", width = 12, justify = "right"))
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
for (i in 1:length(x$testshapeR$testvalR)) {
cat(format(paste(x$testshapeR$testvalR[i]), width = 15, justify = "centre"))
cat(format(sprintf("%3.3f", x$testshapeR$stat.shapeR[i]), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testshapeR$pval.shapeR[i]), width = 12, justify = "right"))
cat("\n")
}
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
cat("\n")
}
if (!is.null(x$testshape2$testval2)) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: mu =", width = 15, justify = "left"))
if (is.infinite(x$opt$lp)) {
cat(format("sup |T|", width = 12, justify = "right"))
} else {
cat(format(paste("L", x$opt$lp, " of |T|", sep=""), width = 12, justify = "right"))
}
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
for (i in 1:length(x$testshape2$testval2)) {
cat(format(paste(x$testshape2$testval2[i]), width = 15, justify = "centre"))
cat(format(sprintf("%3.3f", x$testshape2$stat.shape2[i]), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testshape2$pval.shape2[i]), width = 12, justify = "right"))
cat("\n")
}
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
cat("\n")
}
}
if (!is.null(x$testpoly$testpoly)|!all(is.na(x$testmodel$stat.model))) {
cat(paste("Model Specification Tests:")); cat("\n")
cat(paste("degree (p) = ", x$opt$testmodel[1], sep="")); cat(paste(" "))
cat(paste("smooth (s) = ", x$opt$testmodel[2], sep="")); cat("\n")
if (!is.null(x$testpoly$testpoly)) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: mu =", width = 15, justify = "left"))
if (is.infinite(x$opt$lp)) {
cat(format("sup |T|", width = 12, justify = "right"))
} else {
cat(format(paste("L", x$opt$lp, " of |T|", sep=""), width = 12, justify = "right"))
}
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
cat(format(paste("poly. p=", x$testpoly$testpoly, sep=""), width = 15, justify = "left"))
cat(format(sprintf("%3.3f", x$testpoly$stat.poly), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testpoly$pval.poly), width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
}
if (!all(is.na(x$testmodel$stat.model))) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: mu =", width = 15, justify = "left"))
if (is.infinite(x$opt$lp)) {
cat(format("sup |T|", width = 12, justify = "right"))
} else {
cat(format(paste("L", x$opt$lp, " of |T|", sep=""), width = 12, justify = "right"))
}
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
for (i in 1:length(x$testmodel$stat.model)) {
cat(format(paste("model ", i, sep=""), width = 15, justify = "centre"))
cat(format(sprintf("%3.3f", x$testmodel$stat.model[i]), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testmodel$pval.model[i]), width = 12, justify = "right"))
cat("\n")
}
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
}
}
}
Try the binsreg package in your browser
Any scripts or data that you put into this service are public.
binsreg documentation built on Sept. 11, 2024, 7:54 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions summary.CCFFbinstest print.CCFFbinstest binstest
Documented in binstest print.CCFFbinstest summary.CCFFbinstest
################################################################################################################
#'@title Data-Driven Nonparametric Shape Restriction and Parametric Model Specification Testing using Binscatter
#'@description \code{binstest} implements binscatter-based hypothesis testing procedures for parametric functional
#' forms of and nonparametric shape restrictions on the regression function of interest, following the results
#' in \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_AER.pdf}{Cattaneo, Crump, Farrell and Feng (2024a)} and
#' \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_NonlinearBinscatter.pdf}{Cattaneo, Crump, Farrell and Feng (2024b)}.
#' If the binning scheme is not set by the user,
#' the companion function \code{\link{binsregselect}} is used to implement binscatter in a
#' data-driven way and inference procedures are based on robust bias correction.
#' Binned scatter plots based on different methods can be constructed using the companion functions \code{\link{binsreg}},
#' \code{\link{binsqreg}} or \code{\link{binsglm}}.
#'@param y outcome variable. A vector.
#'@param x independent variable of interest. A vector.
#'@param w control variables. A matrix, a vector or a \code{\link{formula}}.
#'@param data an optional data frame containing variables used in the model.
#'@param estmethod estimation method. The default is \code{estmethod="reg"} for tests based on binscatter least squares regression. Other options are \code{"qreg"} for quantile regression and \code{"glm"} for generalized linear regression. If \code{estmethod="glm"}, the option \code{family} must be specified.
#'@param family a description of the error distribution and link function to be used in the generalized linear model when \code{estmethod="glm"}. (See \code{\link{family}} for details of family functions.)
#'@param quantile the quantile to be estimated. A number strictly between 0 and 1.
#'@param deriv derivative order of the regression function for estimation, testing and plotting.
#' The default is \code{deriv=0}, which corresponds to the function itself.
#'@param at value of \code{w} at which the estimated function is evaluated. The default is \code{at="mean"}, which corresponds to
#' the mean of \code{w}. Other options are: \code{at="median"} for the median of \code{w}, \code{at="zero"} for a vector of zeros.
#' \code{at} can also be a vector of the same length as the number of columns of \code{w} (if \code{w} is a matrix) or a data frame containing the same variables as specified in \code{w} (when
#' \code{data} is specified). Note that when \code{at="mean"} or \code{at="median"}, all factor variables (if specified) are excluded from the evaluation (set as zero).
#'@param nolink if true, the function within the inverse link function is reported instead of the conditional mean function for the outcome.
#'@param testmodel a vector or a logical value. It sets the degree of polynomial and the number of smoothness constraints for parametric model specification
#' testing. If \code{testmodel=c(p,s)} is specified, a piecewise polynomial of degree \code{p} with \code{s} smoothness constraints is used.
#' If \code{testmodel=T} or \code{testmodel=NULL} (default) is specified, \code{testmodel=c(1,1)} is used unless the degree \code{p} or the smoothness \code{s}
#' selection is requested via the option \code{pselect} or \code{sselect} (see more details in the explanation of \code{pselect} and \code{sselect}).
#'@param testmodelparfit a data frame or matrix which contains the evaluation grid and fitted values of the model(s) to be
#' tested against. The column contains a series of evaluation points
#' at which the binscatter model and the parametric model of interest are compared with
#' each other. Each parametric model is represented by other columns, which must
#' contain the fitted values at the corresponding evaluation points.
#'@param testmodelpoly degree of a global polynomial model to be tested against.
#'@param testshape a vector or a logical value. It sets the degree of polynomial and the number of smoothness constraints for nonparametric shape restriction
#' testing. If \code{testshape=c(p,s)} is specified, a piecewise polynomial of degree \code{p} with \code{s} smoothness constraints is used.
#' If \code{testshape=T} or \code{testshape=NULL} (default) is specified, \code{testshape=c(1,1)} is used unless the degree \code{p} or smoothness \code{s} selection
#' is requested via the option \code{pselect} or \code{sselect} (see more details in the explanation of \code{pselect} and \code{sselect}).
#'@param testshapel a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#' corresponds to one boundary of a one-sided hypothesis test to the left of the form
#' \code{H0: sup_x mu(x)<=a}.
#'@param testshaper a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#' corresponds to one boundary of a one-sided hypothesis test to the right of the form
#' \code{H0: inf_x mu(x)>=a}.
#'@param testshape2 a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#' corresponds to one boundary of a two-sided hypothesis test of the form
#' \code{H0: sup_x |mu(x)-a|=0}.
#'@param lp an Lp metric used for parametric model specification testing and/or shape restriction testing. The default is \code{lp=Inf}, which
#' corresponds to the sup-norm of the t-statistic. Other options are \code{lp=q} for a positive number \code{q>=1}. Note that
#' \code{lp=Inf} ("sup-norm") has to be used for testing one-sided shape restrictions.
#'@param bins a vector. If \code{bins=c(p,s)}, it sets the piecewise polynomial of degree \code{p} with \code{s} smoothness constraints
#' for data-driven (IMSE-optimal) selection of the partitioning/binning scheme.
#' The default is \code{bins=c(0,0)}, which corresponds to the piecewise constant.
#'@param nbins number of bins for partitioning/binning of \code{x}. If \code{nbins=T} or \code{nbins=NULL} (default) is specified, the number
#' of bins is selected via the companion command \code{\link{binsregselect}} in a data-driven, optimal way whenever possible.
#' If a vector with more than one number is specified, the number of bins is selected within this vector via the companion command \code{\link{binsregselect}}.
#'@param pselect vector of numbers within which the degree of polynomial \code{p} for point estimation is selected.
#' If the selected optimal degree is \code{p}, then piecewise polynomials of degree
#' \code{p+1} are used to conduct testing for nonparametric shape restrictions or parametric model specifications.
#' \emph{Note:} To implement the degree or smoothness selection, in addition to \code{pselect} or \code{sselect},
#' \code{nbins=#} must be specified.
#'@param sselect vector of numbers within which the number of smoothness constraints \code{s} for point estimation is selected.
#' If the selected optimal smoothness is \code{s}, then piecewise polynomials of \code{s+1} smoothness constraints
#' are used to conduct testing for nonparametric shape restrictions or parametric model specifications.
#' If not specified, for each value \code{p} supplied in the option \code{pselect}, only the piecewise polynomial
#' with the maximum smoothness is considered, i.e., \code{s=p}.
#'@param binspos position of binning knots. The default is \code{binspos="qs"}, which corresponds to quantile-spaced
#' binning (canonical binscatter). The other options are \code{"es"} for evenly-spaced binning, or
#' a vector for manual specification of the positions of inner knots (which must be within the range of
#' \code{x}).
#'@param binsmethod method for data-driven selection of the number of bins. The default is \code{binsmethod="dpi"},
#' which corresponds to the IMSE-optimal direct plug-in rule. The other option is: \code{"rot"}
#' for rule of thumb implementation.
#'@param nbinsrot initial number of bins value used to construct the DPI number of bins selector.
#' If not specified, the data-driven ROT selector is used instead.
#'@param randcut upper bound on a uniformly distributed variable used to draw a subsample for bins/degree/smoothness selection.
#' Observations for which \code{runif()<=#} are used. # must be between 0 and 1. By default, \code{max(5000, 0.01n)} observations
#' are used if the samples size \code{n>5000}.
#'@param nsims number of random draws for hypothesis testing. The default is
#' \code{nsims=500}, which corresponds to 500 draws from a standard Gaussian random vector of size
#' \code{[(p+1)*J - (J-1)*s]}. Setting at least \code{nsims=2000} is recommended to obtain the final results.
#'@param simsgrid number of evaluation points of an evenly-spaced grid within each bin used for evaluation of
#' the supremum (infimum or Lp metric) operation needed to construct hypothesis testing
#' procedures. The default is \code{simsgrid=20}, which corresponds to 20 evenly-spaced
#' evaluation points within each bin for approximating the supremum (infimum or Lp metric) operator.
#' Setting at least \code{simsgrid=50} is recommended to obtain the final results.
#'@param simsseed seed for simulation.
#'@param vce procedure to compute the variance-covariance matrix estimator. For least squares regression and generalized linear regression, the allowed options are the same as that for \code{\link{binsreg}} or \code{\link{binsqreg}}.
#' For quantile regression, the allowed options are the same as that for \code{\link{binsqreg}}.
#'@param cluster cluster ID. Used for compute cluster-robust standard errors.
#'@param asyvar if true, the standard error of the nonparametric component is computed and the uncertainty related to control
#' variables is omitted. Default is \code{asyvar=FALSE}, that is, the uncertainty related to control variables is taken into account.
#'@param dfcheck adjustments for minimum effective sample size checks, which take into account number of unique
#' values of \code{x} (i.e., number of mass points), number of clusters, and degrees of freedom of
#' the different stat models considered. The default is \code{dfcheck=c(20, 30)}.
#' See \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_Stata.pdf}{Cattaneo, Crump, Farrell and Feng (2024c)} for more details.
#'@param masspoints how mass points in \code{x} are handled. Available options:
#' \itemize{
#' \item \code{"on"} all mass point and degrees of freedom checks are implemented. Default.
#' \item \code{"noadjust"} mass point checks and the corresponding effective sample size adjustments are omitted.
#' \item \code{"nolocalcheck"} within-bin mass point and degrees of freedom checks are omitted.
#' \item \code{"off"} "noadjust" and "nolocalcheck" are set simultaneously.
#' \item \code{"veryfew"} forces the function to proceed as if \code{x} has only a few number of mass points (i.e., distinct values).
#' In other words, forces the function to proceed as if the mass point and degrees of freedom checks were failed.
#' }
#'@param weights an optional vector of weights to be used in the fitting process. Should be \code{NULL} or
#' a numeric vector. For more details, see \code{\link{lm}}.
#'@param subset optional rule specifying a subset of observations to be used.
#'@param numdist number of distinct values for selection. Used to speed up computation.
#'@param numclust number of clusters for selection. Used to speed up computation.
#'@param estmethodopt a list of optional arguments used by \code{\link[quantreg]{rq}} (for quantile regression) or \code{\link{glm}} (for fitting generalized linear models).
#'@param ... optional arguments to control bootstrapping if \code{estmethod="qreg"} and \code{vce="boot"}. See \code{\link[quantreg]{boot.rq}}.
#'@return \item{\code{testshapeL}}{Results for \code{testshapel}, including: \code{testvalL}, null boundary values;
#' \code{stat.shapeL}, test statistics; and \code{pval.shapeL}, p-value.}
#' \item{\code{testshapeR}}{Results for \code{testshaper}, including: \code{testvalR}, null boundary values;
#' \code{stat.shapeR}, test statistics; and \code{pval.shapeR}, p-value.}
#' \item{\code{testshape2}}{Results for \code{testshape2}, including: \code{testval2}, null boundary values;
#' \code{stat.shape2}, test statistics; and \code{pval.shape2}, p-value.}
#' \item{\code{testpoly}}{Results for \code{testmodelpoly}, including: \code{testpoly}, the degree of global polynomial;
#' \code{stat.poly}, test statistic; \code{pval.poly}, p-value.}
#' \item{\code{testmodel}}{Results for \code{testmodelparfit}, including: \code{stat.model}, test statistics;
#' \code{pval.model}, p-values.}
#' \item{\code{imse.var.rot}}{Variance constant in IMSE, ROT selection.}
#' \item{\code{imse.bsq.rot}}{Bias constant in IMSE, ROT selection.}
#' \item{\code{imse.var.dpi}}{Variance constant in IMSE, DPI selection.}
#' \item{\code{imse.bsq.dpi}}{Bias constant in IMSE, DPI selection.}
#' \item{\code{opt}}{ A list containing options passed to the function, as well as total sample size \code{n},
#' number of distinct values \code{Ndist} in \code{x}, number of clusters \code{Nclust}, and
#' number of bins \code{nbins}.}
#'
#'@author
#' Matias D. Cattaneo, Princeton University, Princeton, NJ. \email{cattaneo@princeton.edu}.
#'
#' Richard K. Crump, Federal Reserve Bank of New York, New York, NY. \email{richard.crump@ny.frb.org}.
#'
#' Max H. Farrell, UC Santa Barbara, Santa Barbara, CA. \email{mhfarrell@gmail.com}.
#'
#' Yingjie Feng (maintainer), Tsinghua University, Beijing, China. \email{fengyingjiepku@gmail.com}.
#'
#'@references
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2024a: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_AER.pdf}{On Binscatter}. American Economic Review 114(5): 1488-1514.
#'
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2024b: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_NonlinearBinscatter.pdf}{Nonlinear Binscatter Methods}. Working Paper.
#'
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2024c: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2024_Stata.pdf}{Binscatter Regressions}. Working Paper.
#'
#'@seealso \code{\link{binsreg}}, \code{\link{binsqreg}}, \code{\link{binsglm}}, \code{\link{binsregselect}}.
#'
#'@examples
#' x <- runif(500); y <- sin(x)+rnorm(500)
#' est <- binstest(y,x, testmodelpoly=1)
#' summary(est)
#'@export
binstest <- function(y, x, w=NULL, data=NULL, estmethod="reg", family=gaussian(),
quantile=NULL, deriv=0, at=NULL, nolink=F,
testmodel=NULL, testmodelparfit=NULL, testmodelpoly=NULL,
testshape=NULL, testshapel=NULL,
testshaper=NULL, testshape2=NULL, lp=Inf,
bins=NULL, nbins=NULL,
pselect=NULL, sselect=NULL,
binspos="qs", binsmethod="dpi", nbinsrot=NULL, randcut=NULL,
nsims=500, simsgrid=20, simsseed=NULL,
vce=NULL, cluster=NULL, asyvar=F,
dfcheck=c(20,30), masspoints="on", weights=NULL, subset=NULL,
numdist=NULL, numclust=NULL, estmethodopt=NULL, ...) {
# param for internal use
qrot <- 2
####################
### prepare data ###
####################
xname <- deparse(substitute(x), width.cutoff = 500L)
yname <- deparse(substitute(y), width.cutoff = 500L)
weightsname <- deparse(substitute(weights), width.cutoff = 500L)
subsetname <- deparse(substitute(subset), width.cutoff = 500L)
clustername <- deparse(substitute(cluster), width.cutoff = 500L)
# extract y, x, w, weights, subset, if needed (w as a matrix, others as vectors)
# generate design matrix for covariates W
w.factor <- NULL
if (is.null(data)) {
if (is.data.frame(y)) y <- y[,1]
if (is.data.frame(x)) x <- x[,1]
if (!is.null(w)) {
if (is.vector(w)) {
w <- as.matrix(w)
if (is.character(w)) {
w <- model.matrix(~w)[,-1,drop=F]
w.factor <- rep(T, ncol(w))
}
} else if (is.formula(w)) {
w.model <- binsreg.model.mat(w)
w <- w.model$design
w.factor <- w.model$factor.colnum
} else if (is.data.frame(w)) {
w <- as.matrix(w)
}
}
} else {
if (xname %in% names(data)) x <- data[,xname]
if (yname %in% names(data)) y <- data[,yname]
if (weightsname != "NULL") if (weightsname %in% names(data)) {
weights <- data[,weightsname]
}
if (subsetname != "NULL") if (subsetname %in% names(data)) {
subset <- data[,subsetname]
}
if (clustername != "NULL") if (clustername %in% names(data)) {
cluster <- data[,clustername]
}
if (deparse(substitute(w), width.cutoff = 500L)[1]!="NULL") {
if (is.formula(w)) {
if (is.data.frame(at)) {
if (ncol(data)!=ncol(at)) {
data[setdiff(names(at), names(data))] <- NA
at[setdiff(names(data), names(at))] <- NA
}
data <- rbind(data, at)
}
w.model <- binsreg.model.mat(w, data=data)
w <- w.model$design
w.factor <- w.model$factor.colnum
if (is.data.frame(at)) {
eval.w <- w[nrow(w),]
w <- w[-nrow(w),, drop=F]
}
} else {
if (deparse(substitute(w), width.cutoff = 500L) %in% names(data)) w <- data[,deparse(substitute(w), width.cutoff = 500L)]
w <- as.matrix(w)
if (is.character(w)) {
w <- model.matrix(~w)[,-1,drop=F]
w.factor <- rep(T, ncol(w))
}
}
}
}
if (!is.null(w)) nwvar <- ncol(w)
else nwvar <- 0
# extract subset
if (!is.null(subset)) {
y <- y[subset]
x <- x[subset]
w <- w[subset, , drop = F]
if (!is.null(cluster)) cluster <- cluster[subset]
if (!is.null(weights)) weights <- weights[subset]
}
na.ok <- complete.cases(x) & complete.cases(y)
if (!is.null(w)) na.ok <- na.ok & complete.cases(w)
if (!is.null(weights)) na.ok <- na.ok & complete.cases(weights)
if (!is.null(cluster)) na.ok <- na.ok & complete.cases(cluster)
y <- y[na.ok]
x <- x[na.ok]
w <- w[na.ok, , drop = F]
weights <- weights[na.ok]
cluster <- cluster[na.ok]
xmin <- min(x); xmax <- max(x)
# evaluation point of w
if (is.null(at)) at <- "mean"
# change method name if needed
if (is.character(family))
family <- get(family, mode = "function", envir = parent.frame())
if (is.function(family))
family <- family()
if (is.null(family$family)) {
print(family)
stop("'family' not recognized")
}
if (!is.null(family)) if (family$family!="gaussian" | family$link!="identity") {
estmethod <- "glm"
}
# extract family obj
if (estmethod=="glm") {
familyname <- family$family
linkinv <- family$linkinv
linkinv.1 <- family$mu.eta # 1st deriv of inverse link
# 2nd derivative
if (family$link=="logit") {
linkinv.2 <- function(z) linkinv.1(z)*(1-2*linkinv(z))
} else if (family$link=="probit") {
linkinv.2 <- function(z) (-z)*linkinv.1(z)
} else if (family$link=="identity") {
linkinv.2 <- function(z) 0
} else if (family$link=="log") {
linkinv.2 <- function(z) linkinv(z)
} else if (family$link=="inverse") {
linkinv.2 <- function(z) 2*z^(-3)
} else if (family$link=="1/mu^2") {
linkinv.2 <- function(z) 0.75*z^(-2.5)
} else if (family$link=="sqrt") {
linkinv.2 <- function(z) 2
} else {
print("The specified link not allowed for computing polynomial-based CIs.")
stop()
}
}
# define the estimation method
if (estmethod=="reg") {
estmethod.name <- "least squares regression"
if (is.null(vce)) vce <- "HC1"
vce.select <- vce
} else if (estmethod=="qreg") {
estmethod.name <- "quantile regression"
if (is.null(vce)) vce <- "nid"
if (vce=="iid") {
vce.select <- "const"
} else {
vce.select <- "HC1"
}
if (is.null(quantile)) quantile <- 0.5
} else if (estmethod=="glm") {
estmethod.name <- "generalized linear model"
if (is.null(vce)) vce <- "HC1"
vce.select <- vce
} else {
print("estmethod incorrectly specified.")
stop()
}
# analyze bins- and degree-related options
tshaON <- tmodON <- F
if (!is.null(testshapel)|!is.null(testshaper)|!is.null(testshape2)) tshaON <- T
if (!is.null(testmodelparfit)|!is.null(testmodelpoly)) tmodON <- T
if (is.logical(testshape)) if (!testshape) testshape <- NULL
if (is.logical(testmodel)) if (!testmodel) testmodel <- NULL
if (is.logical(testshape) & (!tshaON | !is.character(binspos))) testshape <- NULL
if (is.logical(testmodel) & (!tmodON | !is.character(binspos))) testmodel <- NULL
if (!is.character(binspos)) {
if (!is.null(nbins)|!is.null(pselect)|!is.null(sselect)) {
print("nbins, pselect or sselect incorrectly specified.")
stop()
}
}
# 4 cases: select J, select p, user specify both, or an error
if (!is.null(bins)) {
bins.p <- bins[1]
if (length(bins)==1) {
bins.s <- bins.p
} else if (length(bins)==2) {
bins.s <- bins[2]
} else {
print("bins not correctly specified.")
stop()
}
plist <- bins.p; slist <- bins.s
if (is.numeric(nbins)) if (length(nbins)==1) if (nbins!=0) {
print("bins or nbins is incorrectly specified.")
stop()
}
} else {
plist <- pselect; slist <- sselect
bins.p <- bins.s <- NULL
}
len_nbins <- length(nbins); len_p <- length(plist); len_s <- length(slist)
if (len_p==1 & len_s==0) {
slist <- plist; len_s <- 1
}
if (len_p==0 & len_s==1) {
plist <- slist; len_p <- 1
}
# 1st case: select J
selection <- ""
if (is.character(binspos)) {
if (is.logical(nbins)) if (nbins) selection <- "J"
if (len_nbins==1) if (nbins==0) selection <- "J"
if (is.null(nbins)|len_nbins>1|!is.null(bins)) selection <- "J" # turn on J selection
}
if (selection=="J") {
if (len_p>1|len_s>1) {
if (is.null(nbins)) {
print("nbins must be specified for degree/smoothness selection.")
stop()
} else {
print("Only one p and one s are allowed to select # of bins.")
stop()
}
}
if (is.null(plist)) plist <- deriv
if (is.null(slist)) slist <- plist
if (is.null(bins)) {
bins.p <- plist; bins.s <- slist
}
len_p <- len_s <- 1
if (is.null(testshape)) testshape <- c(plist+1, slist+1)
if (is.logical(testshape)) if (testshape) {
testshape <- c(plist+1, slist+1)
}
if (is.null(testmodel)) testmodel <- c(plist+1, slist+1)
if (is.logical(testmodel)) if (testmodel) {
testmodel <- c(plist+1, slist+1)
}
}
# 2nd case: select p (at least for one object)
pselectOK <- F
if (selection!="J" & (is.logical(testshape)|is.null(testshape)) & tshaON) {
pselectOK <- T
}
if (selection!="J" & (is.logical(testmodel)|is.null(testmodel)) & tmodON) {
pselectOK <- T
}
if (pselectOK & (len_p>1|len_s>1) &len_nbins==1) {
selection <- "P"
}
# 3rd case: user specified
if ((len_p<=1 & len_s<=1) & selection!="J") {
selection <- "U"
if (is.null(testshape)) {
if (len_p==1 & len_s==1) testshape <- c(plist+1, slist+1)
else testshape <- c(deriv+1, deriv+1)
}
if (is.null(testmodel)) {
if (len_p==1 & len_s==1) testmodel <- c(plist+1, slist+1)
else testmodel <- c(deriv+1, deriv+1)
}
}
if (selection=="") {
print("Degree, smoothness, or # of bins not correctly specified")
stop()
}
##################################error checking
exit <- 0
if (!is.character(binspos)) {
if (min(binspos)<=xmin|max(binspos)>=xmax) {
print("Knots out of allowed range")
exit <- 1
}
} else {
if (binspos!="es"&binspos!="qs") {
print("binspos incorrectly specified.")
exit <- 1
}
}
if (length(bins)==2) if (bins[1]<bins[2]) {
print("p<s not allowed.")
exit <- 1
}
if (lp<1) {
print("lp has to be no less than 1.")
exit <- 1
}
if (!is.null(testshapel) | !is.null(testshaper)) if (lp!=Inf) {
print("Sup-norm (lp=Inf) has to be used for testing one-sided restrictions.")
exit <- 1
}
if (length(testshape)==2) if (testshape[1]<testshape[2]) {
print("p<s not allowed.")
exit <- 1
}
if (length(testmodel)==2) if (testmodel[1]<testmodel[2]) {
print("p<s not allowed.")
exit <- 1
}
if (length(testshape)>=1) if(!is.logical(testshape)) if (testshape[1]<deriv) {
print("p<deriv not allowed.")
exit <- 1
}
if (length(testmodel)>=1) if(!is.logical(testmodel)) if (testmodel[1]<deriv) {
print("p<deriv not allowed.")
exit <- 1
}
if (binsmethod!="dpi" & binsmethod!="rot") {
print("bin selection method incorrectly specified.")
exit <- 1
}
if (masspoints == "veryfew") {
print("veryfew not allowed for testing.")
exit <- 1
}
if ((length(testshape)>=1)&(length(bins)>=1)) if (!is.logical(testshape)) if (testshape[1]<=bins[1]) {
warning("p for testing <= p for bin selection not suggested.")
}
if ((length(testmodel)>=1)&(length(bins)>=1)) if (!is.logical(testmodel)) if (testmodel[1]<=bins[1]) {
warning("p for testing <= p for bin selection not suggested.")
}
if (exit>0) stop()
if (nsims<2000|simsgrid<50) {
print("Note: Setting at least nsims=2000 and simsgrid=50 is recommended to obtain the final results.")
}
#################################################
localcheck <- massadj <- T
if (masspoints=="on") {
localcheck <- T; massadj <- T
} else if (masspoints=="off") {
localcheck <- F; massadj <- F
} else if (masspoints=="noadjust") {
localcheck <- T; massadj <- F
} else if (masspoints=="nolocalcheck") {
localcheck <- F; massadj <- T
}
# effective size
eN <- N <- length(x)
Ndist <- NA
if (massadj) {
if (!is.null(numdist)) {
Ndist <- numdist
} else {
Ndist <- length(unique(x))
}
eN <- min(eN, Ndist)
}
Nclust <- NA
if (!is.null(cluster)) {
if (!is.null(numclust)) {
Nclust <- numclust
} else {
Nclust <- length(unique(cluster))
}
eN <- min(eN, Nclust)
}
# prepare params
tsha.p <- testshape[1]; tsha.s <- testshape[2]
if (is.logical(tsha.p)) tsha.p <- NULL
if (!is.null(tsha.s)) if (is.na(tsha.s)) tsha.s <- tsha.p
tmod.p <- testmodel[1]; tmod.s <- testmodel[2]
if (is.logical(tmod.p)) tmod.p <- NULL
if (!is.null(tmod.s)) if (is.na(tmod.s)) tmod.s <- tmod.p
if (selection=="J") {
if (!is.null(tsha.p)) if (tsha.p<=plist) {
tsha.p <- bins.p+1; tsha.s <- tsha.p
warning("Degree for testshape has been changed. It must be greater than the degree for bin selection.")
}
if (!is.null(tmod.p)) if (tmod.p<=plist) {
tmod.p <- bins.p+1; tmod.s <- tmod.p
warning("Degree for testmodel has been changed. It must be greater than the degree for bin selection.")
}
}
if (selection=="U") {
warning("Testing procedures are valid when nbins is much larger than the IMSE-optimal choice. Compare your choice with the IMSE-optimal one obtained by binsregselect().")
}
nL <- length(testshapel); nR <- length(testshaper); nT <- length(testshape2)
ntestshape <- nL + nR + nT
if (selection=="U") {
selectmethod <- "User-specified"
} else {
if (binsmethod=="dpi") {
selectmethod <- "IMSE direct plug-in"
} else {
selectmethod <- "IMSE rule-of-thumb"
}
if (selection=="J") selectmethod <- paste(selectmethod, "(select # of bins)")
if (selection=="P") selectmethod <- paste(selectmethod, "(select degree and smoothness)")
}
knot <- NULL
if (!is.character(binspos)) {
nbins <- length(binspos)+1
knot <- c(xmin, binspos, xmax)
position <- "User-specified"
es <- F
} else {
if (binspos == "es") {
es <- T
position <- "Evenly-spaced"
} else {
es <- F
position <- "Quantile-spaced"
}
}
#### bin selection if needed ########
imse.v.rot <- imse.v.dpi <- imse.b.rot <- imse.b.dpi <- NA
if (selection!="U") {
# check if rot can be implemented
if (is.null(bins.p)) binspcheck <- 6
else binspcheck <- bins.p
if (is.null(nbinsrot)) {
if (eN <= dfcheck[1]+binspcheck+1+qrot) {
print("Too small effective sample size for bin selection.")
stop()
}
}
if (is.na(Ndist)) {
Ndist.sel <- NULL
} else {
Ndist.sel <- Ndist
}
if (is.na(Nclust)) {
Nclust.sel <- NULL
} else {
Nclust.sel <- Nclust
}
randcut1k <- randcut
if (is.null(randcut) & N>5000) {
randcut1k <- max(5000/N, 0.01)
warning("To speed up computation, bin/degree selection uses a subsample of roughly max(5000, 0.01n) observations if the sample size n>5000. To use the full sample, set randcut=1.")
}
if (selection=="J") {
binselect <- binsregselect(y, x, w, deriv=deriv,
bins=c(bins.p,bins.s), binspos=binspos, nbins=nbins,
binsmethod=binsmethod, nbinsrot=nbinsrot,
vce=vce.select, cluster=cluster, randcut=randcut1k,
dfcheck=dfcheck, masspoints=masspoints, weights=weights,
numdist=Ndist.sel, numclust=Nclust.sel)
if (is.na(binselect$nbinsrot.regul)) {
print("Bin selection fails.")
stop()
}
if (binsmethod == "rot") {
nbins <- binselect$nbinsrot.regul
imse.v.rot <- binselect$imse.var.rot
imse.b.rot <- binselect$imse.bsq.rot
} else if (binsmethod == "dpi") {
nbins <- binselect$nbinsdpi
imse.v.dpi <- binselect$imse.var.dpi
imse.b.dpi <- binselect$imse.bsq.dpi
if (is.na(nbins)) {
warning("DPI selection fails. ROT choice used.")
nbins <- binselect$nbinsrot.regul
imse.v.rot <- binselect$imse.var.rot
imse.b.rot <- binselect$imse.bsq.rot
}
}
} else if (selection=="P") {
binselect <- binsregselect(y, x, w, deriv=deriv,
binspos=binspos, nbins=nbins,
pselect=plist, sselect=slist,
binsmethod=binsmethod, nbinsrot=nbinsrot,
vce=vce.select, cluster=cluster, randcut=randcut1k,
dfcheck=dfcheck, masspoints=masspoints, weights=weights,
numdist=Ndist.sel, numclust=Nclust.sel)
if (is.na(binselect$prot.regul)) {
print("Bin selection fails.")
stop()
}
if (binsmethod == "rot") {
bins.p <- binselect$prot.regul
bins.s <- binselect$srot.regul
imse.v.rot <- binselect$imse.var.rot
imse.b.rot <- binselect$imse.bsq.rot
} else if (binsmethod == "dpi") {
bins.p <- binselect$pdpi
bins.s <- binselect$sdpi
imse.v.dpi <- binselect$imse.var.dpi
imse.b.dpi <- binselect$imse.bsq.dpi
if (is.na(bins.p)) {
warning("DPI selection fails. ROT choice used.")
bins.p <- binselect$prot.regul
bins.s <- binselect$srot.regul
imse.v.rot <- binselect$imse.var.rot
imse.b.rot <- binselect$imse.bsq.rot
}
}
if (is.logical(testshape)|is.null(testshape)) {
tsha.p <- bins.p+1; tsha.s <- bins.s+1
} else {
if (tsha.p<=bins.p) {
tsha.p <- bins.p+1; tsha.s <- tsha.p
warning("Degree for testshape has been changed. It must be greater than the IMSE-optimal degree.")
}
}
if (is.logical(testmodel)|is.null(testmodel)) {
tmod.p <- bins.p+1; tmod.s <- bins.s+1
} else {
if (tmod.p<=bins.p) {
tmod.p <- bins.p+1; tmod.s <- tmod.p
warning("Degree for testmodel has been changed. It must be greater than the IMSE-optimal degree.")
}
}
}
}
# check eff size for testing
tsha.fewobs <- F
if (ntestshape!=0) if ((nbins-1)*(tsha.p-tsha.s+1)+tsha.p+1+dfcheck[2]>=eN) {
tsha.fewobs <- T
warning("Too small eff. sample size for testing shape.")
}
tmod.fewobs <- F
if (!is.null(testmodelparfit)|!is.null(testmodelpoly)) if ((nbins-1)*(tmod.p-tmod.s+1)+tmod.p+1+dfcheck[2]>=eN) {
tmod.fewobs <- T
warning("Too small eff. sample size for testing models.")
}
# Generate knot
if (is.null(knot)) {
if (es) {
knot <- genKnot.es(xmin, xmax, nbins)
} else {
knot <- genKnot.qs(x, nbins)
}
}
# check local mass points
knot <- c(knot[1], unique(knot[-1]))
if (nbins!=length(knot)-1) {
warning("Repeated knots. Some bins dropped.")
nbins <- length(knot)-1
}
if (localcheck) {
uniqmin <- binsreg.checklocalmass(x, nbins, es, knot=knot) # mimic STATA
if (ntestshape != 0) {
if (uniqmin < tsha.p+1) {
tsha.fewobs <- T
warning("Some bins have too few distinct values of x for testing shape.")
}
}
if (!is.null(testmodelparfit)|!is.null(testmodelpoly)) {
if (uniqmin < tmod.p+1) {
tmod.fewobs <- T
warning("Some bins have too few distinct values of x for testing models.")
}
}
}
# adjust w variables
if (!is.null(w)) {
if (is.character(at)) {
if (at=="mean") {
eval.w <- colWeightedMeans(x=w, w=weights)
if (!is.null(w.factor)) eval.w[w.factor] <- 0
} else if (at=="median") {
eval.w <- colWeightedMedians(x=w, w=weights)
if (!is.null(w.factor)) eval.w[w.factor] <- 0
} else if (at=="zero") {
eval.w <- rep(0, nwvar)
}
} else if (is.vector(at)) {
eval.w <- at
} else if (is.data.frame(at)) {
eval.w <- eval.w
}
} else {
eval.w <- NULL
}
# seed
if (!is.null(simsseed)) set.seed(simsseed)
# prepare grid for uniform inference
x.grid <- binsreg.grid(knot, simsgrid)$eval
#####################################
####### Shape restriction test ######
#####################################
stat.shapeL <- pval.shapeL <- NA
stat.shapeR <- pval.shapeR <- NA
stat.shape2 <- pval.shape2 <- NA
if (nL>0) {
stat.shapeL <- matrix(NA, nL, 2)
pval.shapeL <- c()
}
if (nR>0) {
stat.shapeR <- matrix(NA, nR, 2)
pval.shapeR <- c()
}
if (nT>0) {
stat.shape2 <- matrix(NA, nT, 2)
pval.shape2 <- c()
}
if (ntestshape != 0 & !tsha.fewobs) {
B <- binsreg.spdes(eval=x, p=tsha.p, s=tsha.s, knot=knot, deriv=0)
k <- ncol(B)
P <- cbind(B, w)
if (estmethod=="reg") {
model <- lm(y ~ P-1, weights=weights)
} else if (estmethod=="qreg") {
model <- do.call(rq, c(list(formula=y ~ P-1, tau=quantile, weights=weights), estmethodopt))
} else if (estmethod=="glm") {
model <- do.call(glm, c(list(formula=y ~ P-1, family=family, weights=weights), estmethodopt))
}
beta <- model$coeff[1:k]
basis.sha <- binsreg.spdes(eval=x.grid, p=tsha.p, s=tsha.s, knot=knot, deriv=deriv)
if (estmethod=="glm" & (!nolink)) {
pred.sha <- binsreg.pred(X=basis.sha, model=model, type="all",
vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
basis.0 <- binsreg.spdes(eval=x.grid, p=tsha.p, s=tsha.s, knot=knot, deriv=0)
fit.0 <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
pred.sha.0 <- linkinv.1(fit.0)
if (asyvar | deriv==0) {
pred.sha$se <- pred.sha.0 * pred.sha$se
if (deriv == 0) pred.sha$fit <- linkinv(pred.sha$fit)
if (deriv == 1) pred.sha$fit <- pred.sha.0 * pred.sha$fit
} else {
basis.sha.0 <- basis.0
basis.sha.1 <- basis.sha
if (!is.null(eval.w)) {
basis.sha.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
basis.sha.1 <- cbind(basis.sha.1, outer(rep(1, nrow(basis.sha.1)), rep(0, nwvar)))
}
basis.all <- linkinv.2(fit.0)*pred.sha$fit*basis.sha.0 + pred.sha.0*basis.sha.1
pred.sha$fit <- pred.sha.0 * pred.sha$fit
pred.sha$se <- binsreg.pred(basis.all, model=model, type="se",
vce=vce, cluster=cluster, avar=T)$se
}
} else {
if (estmethod=="qreg") {
pred.sha <- binsreg.pred(basis.sha, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w,
is.qreg=TRUE, avar=asyvar, ...)
} else {
pred.sha <- binsreg.pred(basis.sha, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
}
}
fit.sha <- pred.sha$fit
se.sha <- pred.sha$se
pos <- !is.na(beta)
k.new <- sum(pos)
if (estmethod=="qreg") vcv.sha <- binsreg.vcov(model, type=vce, cluster=cluster, is.qreg=TRUE, ...)[1:k.new, 1:k.new]
else vcv.sha <- binsreg.vcov(model, type=vce, cluster=cluster)[1:k.new, 1:k.new]
for (j in 1:ntestshape) {
if (j <= nL) {
stat.shapeL[j,2] <- 1
stat.shapeL[j,1] <- max((fit.sha - testshapel[j]) / se.sha)
} else if (j <= nL+nR) {
stat.shapeR[j-nL,2] <- 2
stat.shapeR[j-nL,1] <- min((fit.sha - testshaper[j-nL]) / se.sha)
} else {
stat.shape2[j-nL-nR,2] <- 3
if (is.infinite(lp)) {
stat.shape2[j-nL-nR,1] <- max(abs((fit.sha - testshape2[j-nL-nR]) / se.sha))
} else {
stat.shape2[j-nL-nR,1] <- mean(abs((fit.sha - testshape2[j-nL-nR]) / se.sha)^lp)^(1/lp)
}
}
}
stat.shape <- NULL
if (nL > 0) {
stat.shape <- rbind(stat.shape, stat.shapeL)
}
if (nR > 0) {
stat.shape <- rbind(stat.shape, stat.shapeR)
}
if (nT > 0) {
stat.shape <- rbind(stat.shape, stat.shape2)
}
# Compute p-val
Sigma.root <- lssqrtm(vcv.sha)
num <- basis.sha[,pos,drop=F] %*% Sigma.root
denom.sha <- sqrt(rowSums((basis.sha[, pos, drop=F] %*% vcv.sha) * basis.sha[, pos, drop=F]))
pval.shape <- binsreg.pval(num, denom.sha, nsims, tstat=stat.shape, side=NULL, lp=lp)$pval
if (nL!=0) {
stat.shapeL <- stat.shapeL[,1]
pval.shapeL <- pval.shape[1:nL]
}
if (nR!=0) {
stat.shapeR <- stat.shapeR[,1]
pval.shapeR <- pval.shape[(nL+1):(nL+nR)]
}
if (nT!=0) {
stat.shape2 <- stat.shape2[,1]
pval.shape2 <- pval.shape[(nL+nR+1):ntestshape]
}
}
############################
#### Specification Test ####
############################
stat.poly <- pval.poly <- NA
stat.mod <- pval.mod <- NA
if (!is.null(testmodelpoly)) {
stat.poly <- matrix(NA, 1,2)
pval.poly <- c()
}
if (!is.null(testmodelparfit)) {
stat.mod <- matrix(NA,ncol(testmodelparfit)-1,2)
pval.mod <- c()
}
if ((!is.null(testmodelparfit)|!is.null(testmodelpoly))&!tmod.fewobs) {
if (tmod.p==tsha.p & tmod.s==tsha.s & exists("vcv.sha")) {
exist.mod <- T
#beta.mod <- beta.sha
vcv.mod <- vcv.sha
fit.mod <- fit.sha
se.mod <- se.sha
basis.mod <- basis.sha
denom.mod <- denom.sha
pos <- pos
k.new <- k.new
model <- model
} else {
exist.mod <- F
B <- binsreg.spdes(eval=x, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
k <- ncol(B)
P <- cbind(B, w)
if (estmethod=="reg") {
model <- lm(y ~ P-1, weights=weights)
} else if (estmethod=="qreg") {
model <- do.call(rq, c(list(formula=y ~ P-1, tau=quantile, weights=weights), estmethodopt))
} else {
model <- do.call(glm, c(list(formula=y ~ P-1, family=family, weights=weights), estmethodopt))
}
beta <- model$coeff[1:k]
pos <- !is.na(beta)
# beta.mod <- beta[pos]
k.new <- sum(pos)
if (estmethod=="qreg") vcv.mod <- binsreg.vcov(model, type=vce, cluster=cluster, is.qreg = TRUE, ...)[1:k.new, 1:k.new]
else vcv.mod <- binsreg.vcov(model, type=vce, cluster=cluster)[1:k.new, 1:k.new]
}
######################
# Test poly reg
if (!is.null(testmodelpoly)) {
if (!exist.mod) {
basis.mod <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=deriv)
if (estmethod=="glm" & (!nolink)) {
pred.mod <- binsreg.pred(X=basis.mod, model=model, type="all",
vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
basis.0 <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
fit.0 <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
pred.mod.0 <- linkinv.1(fit.0)
if (asyvar | deriv==0) {
pred.mod$se <- pred.mod.0 * pred.mod$se
if (deriv == 0) pred.mod$fit <- linkinv(pred.mod$fit)
if (deriv == 1) pred.mod$fit <- pred.mod.0 * pred.mod$fit
} else {
basis.mod.0 <- basis.0
basis.mod.1 <- basis.mod
if (!is.null(eval.w)) {
basis.mod.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
basis.mod.1 <- cbind(basis.mod.1, outer(rep(1, nrow(basis.mod.1)), rep(0, nwvar)))
}
basis.all <- linkinv.2(fit.0)*pred.mod$fit*basis.mod.0 + pred.mod.0*basis.mod.1
pred.mod$fit <- pred.mod.0 * pred.mod$fit
pred.mod$se <- binsreg.pred(basis.all, model=model, type="se",
vce=vce, cluster=cluster, avar=T)$se
}
} else {
if (estmethod=="qreg") {
pred.mod <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w,
is.qreg=TRUE, avar=asyvar, ...)
} else {
pred.mod <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
}
}
fit.mod <- pred.mod$fit
se.mod <- pred.mod$se
denom.mod <- sqrt(rowSums((basis.mod[, pos, drop=F] %*% vcv.mod) * basis.mod[, pos, drop=F]))
}
# Run a poly reg
x.p <- matrix(NA, N, testmodelpoly+1)
for (j in 1:(testmodelpoly+1)) x.p[,j] <- x^(j-1)
P.poly <- cbind(x.p, w)
if (estmethod=="reg") {
model.poly <- lm(y~P.poly-1, weights=weights)
} else if (estmethod=="qreg") {
model.poly <- do.call(rq, c(list(formula=y~P.poly-1, tau=quantile, weights=weights), estmethodopt))
} else if (estmethod=="glm") {
model.poly <- do.call(glm, c(list(formula=y~P.poly-1, family=family, weights=weights), estmethodopt))
}
beta.poly <- model.poly$coefficients
poly.fit <- 0
for (j in deriv:testmodelpoly) {
poly.fit <- poly.fit + x.grid^(j-deriv)*beta.poly[j+1]*factorial(j)/factorial(j-deriv)
}
if (!is.null(eval.w) & deriv==0) poly.fit <- poly.fit + sum(eval.w * beta.poly[-(1:(testmodelpoly+1))])
stat.poly[1,2] <- 3
if (is.infinite(lp)) stat.poly[1,1] <- max(abs((fit.mod - poly.fit) / se.mod))
else stat.poly[1,1] <- mean(abs((fit.mod - poly.fit) / se.mod)^lp)^(1/lp)
Sigma.root <- lssqrtm(vcv.mod)
num <- basis.mod[,pos,drop=F] %*% Sigma.root
pval.poly <- binsreg.pval(num, denom.mod, nsims, tstat=stat.poly, side=NULL, lp=lp)$pval
stat.poly <- stat.poly[,1]
}
##################################
# Test model in another data frame
if (!is.null(testmodelparfit)) {
x.grid <- testmodelparfit[,1]
basis.mod <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=deriv)
if (estmethod=="glm" & (!nolink)) {
pred.mod <- binsreg.pred(X=basis.mod, model=model, type="all",
vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
basis.0 <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
fit.0 <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
pred.mod.0 <- linkinv.1(fit.0)
if (asyvar | deriv==0) {
pred.mod$se <- pred.mod.0 * pred.mod$se
if (deriv == 0) pred.mod$fit <- linkinv(pred.mod$fit)
if (deriv == 1) pred.mod$fit <- pred.mod.0 * pred.mod$fit
} else {
basis.mod.0 <- basis.0
basis.mod.1 <- basis.mod
if (!is.null(eval.w)) {
basis.mod.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
basis.mod.1 <- cbind(basis.mod.1, outer(rep(1, nrow(basis.mod.1)), rep(0, nwvar)))
}
basis.all <- linkinv.2(fit.0)*pred.mod$fit*basis.mod.0 + pred.mod.0*basis.mod.1
pred.mod$fit <- pred.mod.0 * pred.mod$fit
pred.mod$se <- binsreg.pred(basis.all, model=model, type="se",
vce=vce, cluster=cluster, avar=T)$se
}
} else {
if (estmethod=="qreg") {
pred.mod <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w,
is.qreg=TRUE, avar=asyvar, ...)
} else {
pred.mod <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
}
}
fit.mod <- pred.mod$fit
se.mod <- pred.mod$se
denom.mod <- sqrt(rowSums((basis.mod[, pos, drop=F] %*% vcv.mod) * basis.mod[, pos, drop=F]))
for (j in 2:ncol(testmodelparfit)) {
stat.mod[j-1,2] <- 3
if (is.infinite(lp)) stat.mod[j-1,1] <- max(abs((fit.mod - testmodelparfit[,j]) / se.mod))
else stat.mod[j-1,1] <- mean(abs((fit.mod - testmodelparfit[,j]) / se.mod)^lp)^(1/lp)
}
Sigma.root <- lssqrtm(vcv.mod)
num <- basis.mod[,pos,drop=F] %*% Sigma.root
pval.mod <- binsreg.pval(num, denom.mod, nsims, tstat=stat.mod, side=NULL, lp=lp)$pval
stat.mod <- stat.mod[,1]
}
}
######################
#######output#########
######################
out <- list(testshapeL=list(testvalL=testshapel, stat.shapeL=stat.shapeL, pval.shapeL=pval.shapeL),
testshapeR=list(testvalR=testshaper, stat.shapeR=stat.shapeR, pval.shapeR=pval.shapeR),
testshape2=list(testval2=testshape2, stat.shape2=stat.shape2, pval.shape2=pval.shape2),
testpoly=list(testpoly=testmodelpoly, stat.poly=stat.poly, pval.poly=pval.poly),
testmodel=list(stat.model=stat.mod, pval.model=pval.mod),
imse.var.dpi=imse.v.dpi, imse.bsq.dpi=imse.b.dpi,
imse.var.rot=imse.v.rot, imse.bsq.rot=imse.b.rot,
opt = list(bins.p=bins.p, bins.s=bins.s, deriv=deriv,
testshape=c(tsha.p, tsha.s), testmodel=c(tmod.p, tmod.s),
binspos=position, binsmethod=selectmethod,
n=N, Ndist=Ndist, Nclust=Nclust,
nbins=nbins, knot=knot, lp=lp,
estmethod=estmethod.name, quantile=quantile, family=family))
out$call <- match.call()
class(out) <- "CCFFbinstest"
return(out)
}
##########################################################################
#' Internal function.
#'
#' @param x Class \code{CCFFbinstest} objects.
#'
#' @keywords internal
#' @export
#'
print.CCFFbinstest <- function(x, ...) {
cat("Call: binstest\n\n")
cat(paste("Sample size (n) = ", x$opt$n, "\n", sep=""))
cat(paste("# of distinct values (Ndist) = ", x$opt$Ndist, "\n", sep=""))
cat(paste("# of clusters (Nclust) = ", x$opt$Nclust, "\n", sep=""))
cat(paste("Estimation Method (estmethod) = ", x$opt$estmethod, "\n", sep=""))
cat(paste("Derivative (deriv) = ", x$opt$deriv, "\n", sep=""))
cat(paste("Bin/Degree selection:", "\n"))
cat(paste(" Method (binsmethod) = ", x$opt$binsmethod,"\n", sep=""))
cat(paste(" Placement (binspos) = ", x$opt$binspos, "\n", sep=""))
cat(paste(" degree (p) = ", x$opt$bins.p, "\n", sep=""))
cat(paste(" smooth (s) = ", x$opt$bins.s, "\n", sep=""))
cat(paste(" # of bins (nbins) = ", x$opt$nbins, "\n", sep=""))
cat("\n")
}
################################################################################
#' Internal function.
#'
#' @param object Class \code{CCFFbinstest} objects.
#'
#' @keywords internal
#' @export
summary.CCFFbinstest <- function(object, ...) {
x <- object
args <- list(...)
cat("Call: binstest\n\n")
cat(paste("Sample size (n) = ", x$opt$n, "\n", sep=""))
cat(paste("# of distinct values (Ndist) = ", x$opt$Ndist, "\n", sep=""))
cat(paste("# of clusters (Nclust) = ", x$opt$Nclust, "\n", sep=""))
cat(paste("Estimation Method (estmethod) = ", x$opt$estmethod, "\n", sep=""))
cat(paste("Derivative (deriv) = ", x$opt$deriv, "\n", sep=""))
cat(paste("Bin/Degree selection:", "\n"))
cat(paste(" Method (binsmethod) = ", x$opt$binsmethod,"\n", sep=""))
cat(paste(" Placement (binspos) = ", x$opt$binspos, "\n", sep=""))
cat(paste(" degree (p) = ", x$opt$bins.p, "\n", sep=""))
cat(paste(" smooth (s) = ", x$opt$bins.s, "\n", sep=""))
cat(paste(" # of bins (nbins) = ", x$opt$nbins, "\n", sep=""))
cat("\n")
if (!is.null(x$testshapeL$testvalL)|!is.null(x$testshapeR$testvalR)|!is.null(x$testshape2$testval2)) {
cat(paste("Shape Restriction Tests:")); cat("\n")
cat(paste("degree (p) = ", x$opt$testshape[1], sep="")); cat(paste(" "))
cat(paste("smooth (s) = ", x$opt$testshape[2], sep="")); cat("\n")
if (!is.null(x$testshapeL$testvalL)) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: sup mu <=", width = 15, justify = "left"))
cat(format("sup T", width = 12, justify = "right"))
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
for (i in 1:length(x$testshapeL$testvalL)) {
cat(format(paste(x$testshapeL$testvalL[i]), width = 15, justify = "centre"))
cat(format(sprintf("%3.3f", x$testshapeL$stat.shapeL[i]), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testshapeL$pval.shapeL[i]), width = 12, justify = "right"))
cat("\n")
}
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
cat("\n")
}
if (!is.null(x$testshapeR$testvalR)) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: inf mu >=", width = 15, justify = "left"))
cat(format("inf T", width = 12, justify = "right"))
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
for (i in 1:length(x$testshapeR$testvalR)) {
cat(format(paste(x$testshapeR$testvalR[i]), width = 15, justify = "centre"))
cat(format(sprintf("%3.3f", x$testshapeR$stat.shapeR[i]), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testshapeR$pval.shapeR[i]), width = 12, justify = "right"))
cat("\n")
}
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
cat("\n")
}
if (!is.null(x$testshape2$testval2)) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: mu =", width = 15, justify = "left"))
if (is.infinite(x$opt$lp)) {
cat(format("sup |T|", width = 12, justify = "right"))
} else {
cat(format(paste("L", x$opt$lp, " of |T|", sep=""), width = 12, justify = "right"))
}
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
for (i in 1:length(x$testshape2$testval2)) {
cat(format(paste(x$testshape2$testval2[i]), width = 15, justify = "centre"))
cat(format(sprintf("%3.3f", x$testshape2$stat.shape2[i]), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testshape2$pval.shape2[i]), width = 12, justify = "right"))
cat("\n")
}
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
cat("\n")
}
}
if (!is.null(x$testpoly$testpoly)|!all(is.na(x$testmodel$stat.model))) {
cat(paste("Model Specification Tests:")); cat("\n")
cat(paste("degree (p) = ", x$opt$testmodel[1], sep="")); cat(paste(" "))
cat(paste("smooth (s) = ", x$opt$testmodel[2], sep="")); cat("\n")
if (!is.null(x$testpoly$testpoly)) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: mu =", width = 15, justify = "left"))
if (is.infinite(x$opt$lp)) {
cat(format("sup |T|", width = 12, justify = "right"))
} else {
cat(format(paste("L", x$opt$lp, " of |T|", sep=""), width = 12, justify = "right"))
}
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
cat(format(paste("poly. p=", x$testpoly$testpoly, sep=""), width = 15, justify = "left"))
cat(format(sprintf("%3.3f", x$testpoly$stat.poly), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testpoly$pval.poly), width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
}
if (!all(is.na(x$testmodel$stat.model))) {
cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
cat(format("H0: mu =", width = 15, justify = "left"))
if (is.infinite(x$opt$lp)) {
cat(format("sup |T|", width = 12, justify = "right"))
} else {
cat(format(paste("L", x$opt$lp, " of |T|", sep=""), width = 12, justify = "right"))
}
cat(format("p value", width = 12, justify = "right"))
cat("\n")
cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
cat("\n")
for (i in 1:length(x$testmodel$stat.model)) {
cat(format(paste("model ", i, sep=""), width = 15, justify = "centre"))
cat(format(sprintf("%3.3f", x$testmodel$stat.model[i]), width = 12, justify = "right"))
cat(format(sprintf("%3.3f", x$testmodel$pval.model[i]), width = 12, justify = "right"))
cat("\n")
}
cat(paste(rep("-", 15 + 12 + 12), collapse=""))
cat("\n")
}
}
}
Try the binsreg package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.