Nothing
################################################################################################################
#'@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_2023_AER.pdf}{Cattaneo, Crump, Farrell and Feng (2023a)} and
#' \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2023_NonlinearBinscatter.pdf}{Cattaneo, Crump, Farrell and Feng (2023b)}.
#' 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 (two-sided) 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 integer \code{q}.
#'@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_2023_Stata.pdf}{Cattaneo, Crump, Farrell and Feng (2023c)} 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 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{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{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. 2023a: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2023_AER.pdf}{On Binscatter}. Working Paper.
#'
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2023b: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2023_NonlinearBinscatter.pdf}{Nonlinear Binscatter Methods}. Working Paper.
#'
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2023c: \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2023_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.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 (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.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.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.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")
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.