diff --git a/.Rbuildignore b/.Rbuildignore index fd50189..5267572 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -4,3 +4,4 @@ src/deprecated ^\.Rproj\.user$ ^\.travis\.yml$ ^appveyor\.yml$ +^tests/testthat/test_krlogit_fns\.R$ diff --git a/.travis.yml b/.travis.yml index 250d97c..ba0b05d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,6 +1,37 @@ -# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r - language: R sudo: false cache: packages -warnings_are_errors: false \ No newline at end of file +warnings_are_errors: false + +matrix: + include: + - os: linux + r: release + + - os: linux + r: oldrel + after_success: + - echo skipping source packaging on linux/oldrel + + - os: linux + r: devel + after_success: + - echo skipping source packaging on linux/devel + + - os: osx + r: release + if: branch = master + + - os: osx + osx_image: xcode6.4 + r: oldrel + if: branch = master + +r_packages: +- numDeriv + +r_github_packages: +- r-lib/covr + +after_success: +- test $TRAVIS_OS_NAME == "linux" && Rscript -e 'covr::coveralls()' diff --git a/DESCRIPTION b/DESCRIPTION index 1c29651..a4f1d2e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,11 @@ Package: KRLS2 Type: Package -Title: Kernel-based Regularized Least squares (KRLS) -Version: 0.4 -Date: 2016-05-15 +Title: Kernel-based Regularized Least squares +Version: 1.1.0 +Date: 2018-02-08 Author: Jens Hainmueller (Stanford) Chad Hazlett (UCLA) Luke Sonnet (UCLA) Maintainer: Jens Hainmueller -Description: Package implements Kernel-based Regularized Least Squares (KRLS), a +Description: Implements Kernel-based Regularized Least Squares (KRLS), a machine learning method to fit multidimensional functions y=f(x) for regression and classification problems without relying on linearity or additivity assumptions. KRLS finds the best fitting function by minimizing the squared loss @@ -17,6 +17,7 @@ Imports: LinkingTo: Rcpp, RcppArmadillo License: GPL (>= 2) Suggests: - lattice, - testthat + testthat, + lattice +URL: https://www.r-project.org, https://www.stanford.edu/~jhain/ RoxygenNote: 6.0.1 diff --git a/NAMESPACE b/NAMESPACE index c192916..1fe4530 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -37,4 +37,16 @@ export(summary.krls2) export(trace_mat) import(RSpectra) importFrom(Rcpp,sourceCpp) -useDynLib(KRLS2) +importFrom(grDevices,devAskNewPage) +importFrom(graphics,arrows) +importFrom(graphics,par) +importFrom(graphics,plot) +importFrom(stats,as.formula) +importFrom(stats,optim) +importFrom(stats,optimize) +importFrom(stats,predict) +importFrom(stats,pt) +importFrom(stats,quantile) +importFrom(stats,sd) +importFrom(stats,var) +useDynLib(KRLS2, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index e5bb2dc..465b6c4 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -3,81 +3,81 @@ #' @export mult_diag <- function(x, d) { - .Call('_KRLS2_mult_diag', PACKAGE = 'KRLS2', x, d) + .Call(`_KRLS2_mult_diag`, x, d) } #' @export trace_mat <- function(x) { - .Call('_KRLS2_trace_mat', PACKAGE = 'KRLS2', x) + .Call(`_KRLS2_trace_mat`, x) } #' @export krls_gr_trunc <- function(U, D, y, w, fitted, dhat, lambda) { - .Call('_KRLS2_krls_gr_trunc', PACKAGE = 'KRLS2', U, D, y, w, fitted, dhat, lambda) + .Call(`_KRLS2_krls_gr_trunc`, U, D, y, w, fitted, dhat, lambda) } #' @export krls_hess_trunc_inv <- function(U, D, w, lambda) { - .Call('_KRLS2_krls_hess_trunc_inv', PACKAGE = 'KRLS2', U, D, w, lambda) + .Call(`_KRLS2_krls_hess_trunc_inv`, U, D, w, lambda) } #' @export krlogit_fn_trunc <- function(par, U, D, y, w, lambda) { - .Call('_KRLS2_krlogit_fn_trunc', PACKAGE = 'KRLS2', par, U, D, y, w, lambda) + .Call(`_KRLS2_krlogit_fn_trunc`, par, U, D, y, w, lambda) } #' @export krlogit_gr_trunc <- function(par, U, D, y, w, lambda) { - .Call('_KRLS2_krlogit_gr_trunc', PACKAGE = 'KRLS2', par, U, D, y, w, lambda) + .Call(`_KRLS2_krlogit_gr_trunc`, par, U, D, y, w, lambda) } #' @export partial_logit <- function(K, coef, beta0) { - .Call('_KRLS2_partial_logit', PACKAGE = 'KRLS2', K, coef, beta0) + .Call(`_KRLS2_partial_logit`, K, coef, beta0) } #' @export krlogit_hess_trunc_inv <- function(par, U, D, y, w, lambda) { - .Call('_KRLS2_krlogit_hess_trunc_inv', PACKAGE = 'KRLS2', par, U, D, y, w, lambda) + .Call(`_KRLS2_krlogit_hess_trunc_inv`, par, U, D, y, w, lambda) } #' @export euc_dist <- function(x1, x2) { - .Call('_KRLS2_euc_dist', PACKAGE = 'KRLS2', x1, x2) + .Call(`_KRLS2_euc_dist`, x1, x2) } #' @export kern_gauss_1d <- function(x1, x2, b) { - .Call('_KRLS2_kern_gauss_1d', PACKAGE = 'KRLS2', x1, x2, b) + .Call(`_KRLS2_kern_gauss_1d`, x1, x2, b) } #' @export kern_gauss <- function(x, b) { - .Call('_KRLS2_kern_gauss', PACKAGE = 'KRLS2', x, b) + .Call(`_KRLS2_kern_gauss`, x, b) } #' @export new_gauss_kern <- function(newx, oldx, b) { - .Call('_KRLS2_new_gauss_kern', PACKAGE = 'KRLS2', newx, oldx, b) + .Call(`_KRLS2_new_gauss_kern`, newx, oldx, b) } #' @export solve_for_d_ls <- function(y, U, D, lambda) { - .Call('_KRLS2_solve_for_d_ls', PACKAGE = 'KRLS2', y, U, D, lambda) + .Call(`_KRLS2_solve_for_d_ls`, y, U, D, lambda) } #' @export solve_for_d_ls_w <- function(y, U, D, w, lambda) { - .Call('_KRLS2_solve_for_d_ls_w', PACKAGE = 'KRLS2', y, U, D, w, lambda) + .Call(`_KRLS2_solve_for_d_ls_w`, y, U, D, w, lambda) } #' @export pwmfx <- function(k, x, coefhat, vcovc, p, b) { - .Call('_KRLS2_pwmfx', PACKAGE = 'KRLS2', k, x, coefhat, vcovc, p, b) + .Call(`_KRLS2_pwmfx`, k, x, coefhat, vcovc, p, b) } #' @export pwmfx_novar <- function(k, x, coefhat, p, b) { - .Call('_KRLS2_pwmfx_novar', PACKAGE = 'KRLS2', k, x, coefhat, p, b) + .Call(`_KRLS2_pwmfx_novar`, k, x, coefhat, p, b) } diff --git a/R/inference.R b/R/inference.R index 974ee95..cc48682 100644 --- a/R/inference.R +++ b/R/inference.R @@ -99,24 +99,28 @@ inference.krls2 <- function(obj, if(is.null(clusters)) { score <- matrix(nrow = n, ncol = length(obj$dhat)) for(i in 1:n) { - score[i, ] = krls_gr_trunc(obj$U[i, , drop = F], - obj$D, - y[i], - obj$w[i], - yfitted[i], - obj$dhat, - obj$lambda/n) + score[i, ] <- krls_gr_trunc( + obj$U[i, , drop = F], + obj$D, + y[i], + obj$w[i], + yfitted[i], + obj$dhat, + obj$lambda / n + ) } } else { score <- matrix(nrow = length(clusters), ncol = length(obj$dhat)) for(j in 1:length(clusters)){ - score[j, ] = krls_gr_trunc(obj$U[clusters[[j]], , drop = F], - obj$D, - y[clusters[[j]]], - obj$w[clusters[[j]]], - yfitted[clusters[[j]]], - obj$dhat, - length(clusters[[j]]) * obj$lambda/n) + score[j, ] <- krls_gr_trunc( + obj$U[clusters[[j]], , drop = F], + obj$D, + y[clusters[[j]]], + obj$w[clusters[[j]]], + yfitted[clusters[[j]]], + obj$dhat, + length(clusters[[j]]) * obj$lambda / n + ) } } diff --git a/R/kernels.R b/R/kernels.R index 9c0a9d4..c91483b 100644 --- a/R/kernels.R +++ b/R/kernels.R @@ -29,7 +29,7 @@ generateK <- function(X, U <- truncDat$Utrunc D <- truncDat$eigvals } else { - eigobj <- eigen(K) + eigobj <- eigen(K, symmetric = T) eigvaliszero <- eigobj$values == 0 if(any(eigvaliszero)) { diff --git a/R/krls2.R b/R/krls2.R index 4c562a5..5274952 100644 --- a/R/krls2.R +++ b/R/krls2.R @@ -25,9 +25,11 @@ #' via the summary() function. See summary.krls2(). #' #' @import RSpectra -#' @useDynLib KRLS2 +#' @useDynLib KRLS2, .registration = TRUE #' @importFrom Rcpp sourceCpp - +#' @importFrom grDevices devAskNewPage +#' @importFrom graphics arrows par plot +#' @importFrom stats as.formula optim optimize predict pt quantile sd var ############# # Functions # @@ -54,11 +56,13 @@ #' @param U Positive scalar that determines the upper bound of the search window for the leave-one-out optimization to find \eqn{\lambda}{lambda} with least squares loss. Default is \code{NULL} which means that the upper bound is found by using an algorithm outlined in \code{\link{lambdaline}}. Ignored with logistic loss. #' @param tol Positive scalar that determines the tolerance used in the optimization routine used to find \eqn{\lambda}{lambda} with least squares loss. Default is \code{NULL} which means that convergence is achieved when the difference in the sum of squared leave-one-out errors between the \var{i} and the \var{i+1} iteration is less than \var{N * 10^-3}. Ignored with logistic loss. #' @param truncate A boolean that defaults to \code{FALSE}. If \code{TRUE} truncates the kernel matrix, keeping as many eigenvectors as needed so that 1-\code{epsilon} of the total variance in the kernel matrix is retained. Alternatively, you can simply specify \code{epsilon} and truncation will be used. -#' @param epsilon Scalar between 0 and 1 that determines the total variance that can be lost in truncation. If not NULL, truncation is automatically set to TRUE. +#' @param epsilon Scalar between 0 and 1 that determines the total variance that can be lost in truncation. If not NULL, truncation is automatically set to TRUE. If \code{truncate == TRUE}, default is 0.001. +#' @param lastkeeper Number of columns of \code{U} to keep when \code{truncate == TRUE}. Overrides \code{epsilon}. #' @param con A list of control arguments passed to optimization for the numerical optimization of the kernel regularized logistic loss function. -#' @param returnopt A boolean that defaults to \code{FALSE}. If \code{TRUE}, returns the result of the \code{optim} method called to optimize the kernel regularized logistic loss function. +#' @param returnopt A boolean that defaults to \code{FALSE}. If \code{TRUE}, returns the result of the \code{optim} method called to optimize the kernel regularized logistic loss function. Returns \code{NULL} with leastsquares loss. #' @param printlevel A number that is either 0 (default), 1, or 2. 0 Has minimal printing, 1 prints out most diagnostics, and 2 prints out most diagnostics including \code{optim} diagnostics for each fold in the cross-validation selection of hyperparameters. #' @param warn A number that sets your \code{warn} option. We default to 1 so that warnings print as they occur. You can change this to 2 if you want all warnings to be errors, to 0 if you want all warnings to wait until the top-level call is finished, or to a negative number to ignore them. +#' @param sigma DEPRECATED. Users should now use \code{b}, included for backwards compatability. #' @details #' \code{krls} implements the Kernel-based Regularized Least Squares (KRLS) estimator as described in Hainmueller and Hazlett (2014). Please consult this reference for any details. @@ -121,7 +125,7 @@ krls <- function(# Data arguments lastkeeper = NULL, # Optimization arguments con = list(maxit=500), - returnopt = TRUE, + returnopt = FALSE, printlevel = 0, warn = 1, sigma = NULL, # to provide legacy support for old code, @@ -370,17 +374,19 @@ krls <- function(# Data arguments coefhat <- UDinv %*% out$dhat + opt <- NULL if (loss == "leastsquares") { yfitted <- Kdat$K %*% coefhat yfitted <- yfitted * y.init.sd + y.init.mean - + } else { - opt <- if(returnopt) out$opt else NULL - yfitted <- logistic(K=Kdat$K, coeff=coefhat, beta0 = out$beta0hat) + if (returnopt) { + opt <- out$opt + } } z <- list(K = Kdat$K, @@ -398,7 +404,8 @@ krls <- function(# Data arguments b = b, lambda = lambda, kernel = whichkernel, - loss = loss + loss = loss, + opt = opt ) class(z) <- "krls2" diff --git a/R/plot.krls2.R b/R/plot.krls2.R index ed19137..ca3b1aa 100644 --- a/R/plot.krls2.R +++ b/R/plot.krls2.R @@ -11,7 +11,7 @@ #' max and the min value of the predictor (instead of the range from the 1st #' to the 3rd quantile). #' -#' @param obj an object of class \code{krls2} that preferably results from a call to +#' @param x an object of class \code{krls2} that preferably results from a call to #' \code{\link{summary.krls2}}. Also accepts the output of \code{\link{krls}} but #' this is slower as it has to recompute pointwise marginal effects if the #' histograms of pointwise marginal effects are requested by \code{which}. @@ -68,7 +68,7 @@ #' @export #' plot.krls2 <- - function(obj, + function(x, which = c(1:2), main = "distributions of pointwise marginal effects", setx = "mean", @@ -78,14 +78,14 @@ plot.krls2 <- ...) { - if( class(obj)!= "krls2" ){ - warning("obj not of class 'krls2'") + if( class(x)!= "krls2" ){ + warning("`x` must be of class 'krls2'") #UseMethod("summary") return(invisible(NULL)) } - d <- ncol(obj$X) - n <- nrow(obj$X) + d <- ncol(x$X) + n <- nrow(x$X) if(length(probs)!=2){ stop("length(probs) must be 2") } @@ -98,7 +98,7 @@ plot.krls2 <- } else { if(length(setx)!=1){stop("setx must be one of mean or median")} if(sum(setx %in% c("mean","median"))<1){stop("setx must be one of mean or median")} - setx <- apply(obj$X,2,setx) + setx <- apply(x$X,2,setx) } nplots <- 0 @@ -110,25 +110,26 @@ plot.krls2 <- on.exit(devAskNewPage(oask)) } - if(is.null(colnames(obj$X))){ - colnames(obj$X) <- paste("x",1:d,sep="") + if(is.null(colnames(x$X))){ + colnames(x$X) <- paste("x",1:d,sep="") } # derivatives if(1 %in% which){ # histograms of partial derivatives - if(is.null(obj$derivatives)){ + if(is.null(x$derivatives)){ warning("running summary.krls2 because object passed does not have PWMFX. To save time, pass the result of summary.krls2 instead of krls to this plot method.") - obj <- inference.krls2(obj) + x <- inference.krls2(x) } else { - colnames(obj$derivatives) <- colnames(obj$X) + colnames(x$derivatives) <- colnames(x$X) - form <- as.formula(paste("~",paste(colnames(obj$derivatives),collapse="+"),sep="")) - require(lattice) + form <- as.formula(paste("~",paste(colnames(x$derivatives),collapse="+"),sep="")) + + requireNamespace("lattice", quietly = TRUE) - print(histogram(form, - data=data.frame(obj$derivatives), + print(lattice::histogram(form, + data=data.frame(x$derivatives), breaks=NULL, main=main ,...) @@ -140,13 +141,13 @@ plot.krls2 <- if(2 %in% which){ # conditional expectation plots lengthunique <- function(x){length(unique(x))} # vector with positions of binary variables - binaryindicator <- which(apply(obj$X,2,lengthunique)==2) - quantiles <- apply(obj$X,2,quantile,probs=probs) + binaryindicator <- which(apply(x$X,2,lengthunique)==2) + quantiles <- apply(x$X,2,quantile,probs=probs) for(i in 1:d){ if(i %in% binaryindicator){ # E[Y|X] for binary Xs - Xi <- c(min(obj$X[,i]),max(obj$X[,i])) + Xi <- c(min(x$X[,i]),max(x$X[,i])) Newdata <- matrix(rep(setx,2),ncol=d,byrow=T) Newdata[,i] <- Xi @@ -156,12 +157,12 @@ plot.krls2 <- Newdata <- matrix(rep(setx,nvalues),ncol=d,byrow=T) Newdata[,i] <- Xi } - pout <- predict(obj,newdata=Newdata,se=TRUE) + pout <- predict(x,newdata=Newdata,se=TRUE) Ylo <- pout$fit-1.96*pout$se Yhi <- pout$fit+1.96*pout$se plot( y=pout$fit,x=Xi, - xlab=colnames(obj$X)[i], + xlab=colnames(x$X)[i], ylab=c("E[Y|X]"), ylim=c(min(Ylo) -.25*c(sqrt(var(pout$fit))), max(Yhi))+.25*c(sqrt(var(pout$fit))), diff --git a/R/summary.krls2.R b/R/summary.krls2.R index 1376deb..3b3551e 100644 --- a/R/summary.krls2.R +++ b/R/summary.krls2.R @@ -10,7 +10,7 @@ summary.krls2 <- function(object, probs = c(.25, .5, .75), ...) { if( class(object)!= "krls2" ) { - warning("Object not of class 'krls2'") + warning("`object` must be of class 'krls2'") UseMethod("summary") return(invisible(NULL)) } diff --git a/appveyor.yml b/appveyor.yml index c6c1438..e20b48e 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -14,6 +14,13 @@ cache: - C:\RLibrary # Adapt as necessary starting from here +environment: + + matrix: + - R_VERSION: release + + - R_VERSION: oldrel + RTOOLS_VERSION: 33 build_script: - travis-tool.sh install_deps diff --git a/man/krls.Rd b/man/krls.Rd index 551f356..c9bfd99 100644 --- a/man/krls.Rd +++ b/man/krls.Rd @@ -9,7 +9,8 @@ krls(X, y, w = NULL, loss = "leastsquares", whichkernel = "gaussian", lambda = NULL, hyperfolds = 5, lambdastart = 10^(-6)/length(y), lambdainterval = c(10^-8, 25), L = NULL, U = NULL, tol = NULL, truncate = FALSE, epsilon = NULL, lastkeeper = NULL, con = list(maxit - = 500), returnopt = TRUE, printlevel = 0, warn = 1, sigma = NULL, ...) + = 500), returnopt = FALSE, printlevel = 0, warn = 1, sigma = NULL, + ...) } \arguments{ \item{X}{\emph{N} by \emph{D} data numeric matrix that contains the values of \emph{D} predictor variables for \eqn{i=1,\ldots,N} observations. The matrix may not contain missing values or constants. Note that no intercept is required for the least squares or logistic loss. In the case of least squares, the function operates on demeaned data and subtracting the mean of \emph{y} is equivalent to including an (unpenalized) intercept into the model. In the case of logistic loss, we automatically estimate an unpenalized intercept in the linear component of the model.} @@ -44,15 +45,19 @@ krls(X, y, w = NULL, loss = "leastsquares", whichkernel = "gaussian", \item{truncate}{A boolean that defaults to \code{FALSE}. If \code{TRUE} truncates the kernel matrix, keeping as many eigenvectors as needed so that 1-\code{epsilon} of the total variance in the kernel matrix is retained. Alternatively, you can simply specify \code{epsilon} and truncation will be used.} -\item{epsilon}{Scalar between 0 and 1 that determines the total variance that can be lost in truncation. If not NULL, truncation is automatically set to TRUE.} +\item{epsilon}{Scalar between 0 and 1 that determines the total variance that can be lost in truncation. If not NULL, truncation is automatically set to TRUE. If \code{truncate == TRUE}, default is 0.001.} + +\item{lastkeeper}{Number of columns of \code{U} to keep when \code{truncate == TRUE}. Overrides \code{epsilon}.} \item{con}{A list of control arguments passed to optimization for the numerical optimization of the kernel regularized logistic loss function.} -\item{returnopt}{A boolean that defaults to \code{FALSE}. If \code{TRUE}, returns the result of the \code{optim} method called to optimize the kernel regularized logistic loss function.} +\item{returnopt}{A boolean that defaults to \code{FALSE}. If \code{TRUE}, returns the result of the \code{optim} method called to optimize the kernel regularized logistic loss function. Returns \code{NULL} with leastsquares loss.} \item{printlevel}{A number that is either 0 (default), 1, or 2. 0 Has minimal printing, 1 prints out most diagnostics, and 2 prints out most diagnostics including \code{optim} diagnostics for each fold in the cross-validation selection of hyperparameters.} \item{warn}{A number that sets your \code{warn} option. We default to 1 so that warnings print as they occur. You can change this to 2 if you want all warnings to be errors, to 0 if you want all warnings to wait until the top-level call is finished, or to a negative number to ignore them.} + +\item{sigma}{DEPRECATED. Users should now use \code{b}, included for backwards compatability.} } \description{ This is the primary fitting function. diff --git a/man/plot.krls2.Rd b/man/plot.krls2.Rd index b1c870f..8d9f6e9 100644 --- a/man/plot.krls2.Rd +++ b/man/plot.krls2.Rd @@ -4,13 +4,13 @@ \alias{plot.krls2} \title{Plot method for KRLS and KRLogit Model Fits} \usage{ -\method{plot}{krls2}(obj, which = c(1:2), +\method{plot}{krls2}(x, which = c(1:2), main = "distributions of pointwise marginal effects", setx = "mean", ask = prod(par("mfcol")) < nplots, nvalues = 50, probs = c(0.25, 0.75), ...) } \arguments{ -\item{obj}{an object of class \code{krls2} that preferably results from a call to +\item{x}{an object of class \code{krls2} that preferably results from a call to \code{\link{summary.krls2}}. Also accepts the output of \code{\link{krls}} but this is slower as it has to recompute pointwise marginal effects if the histograms of pointwise marginal effects are requested by \code{which}.} diff --git a/src/rcpparma_fns.cpp b/src/krls_helpers.cpp similarity index 100% rename from src/rcpparma_fns.cpp rename to src/krls_helpers.cpp diff --git a/tests/testthat/test_krlogit_fns.R b/tests/testthat/test_krlogit_fns.R new file mode 100644 index 0000000..02abbf2 --- /dev/null +++ b/tests/testthat/test_krlogit_fns.R @@ -0,0 +1,90 @@ +context("KRLogit gradient/hessian correct") + +# Note that this file is ignored by buildignore/CRAN to avoid +# adding numDeriv to `suggests` in DESCRIPTION + +test_that("", { + library(numDeriv) + + set.seed(42) + + # Fit model + krlog_out <- krls( + y = mtcars$am, + X = mtcars[, c("mpg", "wt")], + loss = "logistic", + lambda = 0.2, + epsilon = 0.01, + returnopt = TRUE + ) + + par <- rnorm(ncol(krlog_out$U) + 1) # +1 for beta0 + U <- krlog_out$U + D <- krlog_out$D + y <- krlog_out$y + w <- krlog_out$w + lambda <- 0.1 + + # ----- + # Check gradient + # ----- + + # numerical deriv + num_grad <- grad( + krlogit_fn_trunc, + x = par, + U = U, + D = D, + y = y, + w = w, + lambda = lambda + ) + + # our grad fn + an_grad <- krlogit_gr_trunc( + par = par, + U = U, + D = D, + y = y, + w = w, + lambda = lambda + ) + + # looks good! + expect_equivalent( + num_grad, + an_grad + ) + + # ----- + # Check hessian + # ----- + + # numerical hess + num_hess <- hessian( + krlogit_fn_trunc, + x = par, + U = U, + D = D, + y = y, + w = w, + lambda = lambda + ) + + # our hess + an_hess <- solve(krlogit_hess_trunc_inv( + par = par, + U = U, + D = D, + y = y, + w = w, + lambda = lambda + )) + + + expect_equivalent( + num_hess, + an_hess + ) + +}) diff --git a/tests/testthat/test_krls_summary_arguments.R b/tests/testthat/test_krls_summary_arguments.R index 2ee45af..30c8dcd 100644 --- a/tests/testthat/test_krls_summary_arguments.R +++ b/tests/testthat/test_krls_summary_arguments.R @@ -1,8 +1,6 @@ -library(testthat) -library(KRLS2) -context("Test different input configurations work\n") +context("KRLS/KRLogit specifications") -n <- 50 +n <- 25 betas <- rnorm(2) X <- cbind(rbinom(n, 1, 0.4), rnorm(n)) ypure <- X %*% betas @@ -55,23 +53,23 @@ for(i in 1:nrow(test_grid)) { lambda <- lambdagrid } - cat('\n') - print(test_grid[i,]) + # cat('\n') + # print(test_grid[i,]) # don't support these, should throw error when fitting KR* if ((test_grid$weighted[i] & test_grid$loss[i] == 'leastsquares' & !test_grid$truncation[i])) { - print('expect fit error') + # print('expect fit error') expecting_fit_error <- 'weighted' #this means we will expect an error with one of these things in it (regex) } else { - print('expect no fit error') + # print('expect no fit error') expecting_fit_error <- NA #this means we will not expect an error } # don't support these, should throw error when getting pwmfx if (test_grid$whichkernel[i] != 'gaussian') { - print('expect inference error') + # print('expect inference error') expecting_inf_error <- 'Robust|kernel' #this means we will expect an error with one of these things in it (regex) } else { - print('expect no inference error') + # print('expect no inference error') expecting_inf_error <- NA #this means we will not expect an error } diff --git a/tests/testthat/test_truncation.R b/tests/testthat/test_truncation.R new file mode 100644 index 0000000..90600a4 --- /dev/null +++ b/tests/testthat/test_truncation.R @@ -0,0 +1,53 @@ +context("Truncation works") + +test_that("Truncation returns roughly same results as no truncation", { + # This largely tests if anything goes really wrong in truncation + # not that it is a good approximation + + krlso <- krls(y = mtcars$am, X = mtcars$mpg, truncate = FALSE) + + krlso_trunc <- krls(y = mtcars$am, X = mtcars$mpg, epsilon = 0.0000001) + + expect_true( + ncol(krlso_trunc$U) < ncol(krlso$U) + ) + + expect_true( + max( + inference.krls2(krlso)$derivatives - + inference.krls2(krlso_trunc)$derivatives + ) < 1e-6 + ) + + expect_true( + max( + krlso$fitted - + krlso_trunc$fitted + ) < 1e-7 + ) + + # Overfitting without truncation with logistic seems to be a serious problem + # TODO resovlve below + # krlogo <- krls(y = mtcars$am, X = mtcars$mpg, truncate = FALSE, loss = "logistic", hyperfolds = nrow(mtcars)) + # + # krlogo_trunc <- krls(y = mtcars$am, X = mtcars$mpg, epsilon = 0.0000001, loss = "logistic", hyperfolds = nrow(mtcars)) + # + # expect_true( + # ncol(krlogo_trunc$U) < ncol(krlogo$U) + # ) + # + # expect_true( + # max( + # inference.krls2(krlogo, vcov = FALSE)$derivatives - + # inference.krls2(krlogo_trunc, vcov = FALSE)$derivatives + # ) < 1e-6 + # ) + # + # expect_equal( + # max( + # krlogo$fitted - + # krlogo_trunc$fitted + # ) < 1e-7 + # ) +}) + \ No newline at end of file