Running 'test-bsplineD.R' [0s]
  Running 'test-bsplines.R' [1s]
  Running 'test-btensor.R' [167s]
  Running 'test-build_tensor.R' [168s]
  Running 'test-cn.R' [8s]
  Running 'test-cnr.R' [5s]
  Running 'test-coef_vcov.R' [2s]
  Running 'test-cp.R' [200s]
  Running 'test-cp_diff.R' [156s]
  Running 'test-cp_value.R' [178s]
  Running 'test-cpr.R' [2s]
  Running 'test-datasets.R' [0s]
  Running 'test-generate_cp_formula_data.R' [159s]
  Running 'test-get_spline.R' [172s]
  Running 'test-get_surface.R' [0s]
  Running 'test-iknots_or_df.R' [0s]
  Running 'test-influence_of_iknots.R' [172s]
  Running 'test-insert_a_knot.R' [0s]
  Running 'test-knot_expr.R' [0s]
  Running 'test-loglikelihood.R' [170s]
  Running 'test-matrix_rank.R' [1s]
  Running 'test-modify-formula-and-data.R' [163s]
  Running 'test-newknots.R' [171s]
  Running 'test-order_statistics.R' [156s]
  Running 'test-plot.cpr_bs.R' [165s]
  Running 'test-plot.cpr_cn.R' [161s]
  Running 'test-plot.cpr_cnr.R' [0s]
  Running 'test-plot.cpr_cp.R' [158s]
  Running 'test-plot.cpr_cpr.R' [0s]
  Running 'test-predict.R' [12s]
  Running 'test-recover-known-spline.R' [7s]
  Running 'test-summary.cpr_cn.R' [159s]
  Running 'test-summary.cpr_cnr.R' [5s]
  Running 'test-summary.cpr_cp.R' [6s]
  Running 'test-summary.cpr_cpr.R' [158s]
  Running 'test-trimmed_quantile.R' [0s]
  Running 'test-update_bsplines.R' [8s]
  Running 'test-wiggle.R' [159s]
Running the tests in 'tests/test-btensor.R' failed.
Complete output:
  > library(cpr)
  > 
  > test_warning <- function(x = NULL, msg = NULL) {
  +   stopifnot(!is.null(x))
  +   stopifnot(!is.null(msg))
  +   stopifnot(inherits(x, "warning"))
  +   stopifnot(identical(x$message,  msg))
  +   invisible(TRUE)
  + }
  > 
  > test_error <- function(x, msg = NULL) {
  +   stopifnot(!is.null(x))
  +   stopifnot(!is.null(msg))
  +   stopifnot(inherits(x, "error"))
  +   stopifnot(identical(x$message, msg))
  +   invisible(TRUE)
  + }
  > 
  > ################################################################################
  > # 1D btensor matrix wraps x as a list with warning
  > e <- new.env()
  > with(e, {
  +   x <- tryCatch(btensor(x = mtcars$hp, bknots = list(c(50, 400))), warning = function(w){w})
  +   test_warning(x, "wrapping x into a list.")
  + 
  +   x <- suppressWarnings(
  +                         btensor(x = mtcars$hp, bknots = list(c(50, 400)))
  +   )
  +   x1 <- btensor(x = list(mtcars$hp), bknots = list(c(50, 400)))
  + 
  +   stopifnot(isTRUE(all.equal(x, x1, check.attributes = FALSE)))
  + 
  + })
  > 
  > ################################################################################
  > # 1D btensor matrix is constructed as expected
  > 
  > e <- new.env()
  > with(e, {
  + 
  +   bm <-
  +     btensor(x = list(mtcars$hp), iknots = list(c(100, 150)), bknots = list(c(50,350)))
  + 
  +   mm <-
  +     model.matrix( ~ 0 +
  +                  splines::bs(mtcars$hp, knots = c(100, 150), Boundary.knots = c(50, 350), intercept = TRUE)
  +     )
  + 
  +   stopifnot(isTRUE(all.equal(mm, unclass(bm), check.attributes = FALSE)))
  + 
  + })
  > 
  > ################################################################################
  > # 2D btensor matrix is constructed as expected
  > 
  > e <- new.env()
  > with(e, {
  +   bm <-
  +     btensor(x = list(mtcars$disp, mtcars$hp),
  +             iknots = list(numeric(0), c(100, 150)),
  +             bknots = list(c(70, 475), c(50, 350))
  +     )
  + 
  +   mm <-
  +     model.matrix( ~ 0 +
  +                  splines::bs(mtcars$disp, intercept = TRUE, Boundary.knots = c(70, 475)) :
  +                  splines::bs(mtcars$hp, knots = c(100, 150), Boundary.knots = c(50, 350), intercept = TRUE)
  +     )
  + 
  +   stopifnot(isTRUE(all.equal(mm, unclass(bm), check.attributes = FALSE)))
  + 
  + })
  > 
  > ################################################################################
  > # 3D btensor matrix is constructed as expected
  > 
  > e <- new.env()
  > with(e, {
  +   bm <-
  +     btensor(x = list(mtcars$disp, mtcars$hp, mtcars$mpg),
  +             iknots = list(numeric(0), c(100, 150), c(12.2, 16.3, 21.9)),
  +             bknots = list(c(70, 475), c(50, 350), c(10, 35))
  +     )
  + 
  +   mm <-
  +     model.matrix( ~ 0 +
  +                  splines::bs(mtcars$disp, intercept = TRUE, Boundary.knots = c(70, 475)) :
  +                  splines::bs(mtcars$hp, knots = c(100, 150), Boundary.knots = c(50, 350), intercept = TRUE) :
  +                  splines::bs(mtcars$mpg, knots = c(12.2, 16.3, 21.9), Boundary.knots = c(10, 35), intercept = TRUE))
  + 
  +   stopifnot(isTRUE(all.equal(mm, unclass(bm), check.attributes = FALSE)))
  + })
  > 
  > 
  > ################################################################################
  > # bknots are constructed as exptected
  > e <- new.env()
  > with(e, {
  +   bm <-
  +     suppressWarnings(
  +        btensor(x = list(mtcars$disp, mtcars$hp, mtcars$mpg),
  +                iknots = list(numeric(0), c(100, 150), c(12.2, 16.3, 21.9)))
  +     )
  + 
  +   stopifnot(isTRUE(
  +               all.equal(
  +                         lapply(attr(bm, "bspline_list"), attr, which = "bknots"),
  +                         lapply(list(mtcars$disp, mtcars$hp, mtcars$mpg), range)
  +               )
  +     )
  +   )
  + })
  > 
  > e <- new.env()
  > with(e, {
  +   x <-
  +     tryCatch(
  +              btensor(x = list(mtcars$disp, mtcars$hp, mtcars$mpg), bknots = c(12, 12, 3))
  +              , error = function(e) {e})
  +   test_error(x, "is.list(bknots) is not TRUE")
  + })
  > 
  > e <- new.env()
  > with(e, {
  +   x <- tryCatch(
  +     btensor(x = list(mtcars$disp, mtcars$hp, mtcars$mpg), bknots = list(12, 12, 3))
  +     , error = function(e) {e})
  +   test_error(x, "all(sapply(bknots, length) == 2) is not TRUE")
  + })
  > 
  > e <- new.env()
  > with(e, {
  +   x <- tryCatch(btensor(x = list(mtcars$disp, mtcars$hp, mtcars$mpg), bknots = list(c(12, 12), c(3, NA))), error = function(e) {e})
  +   test_error(x, "length(bknots) == length(x) is not TRUE")
  + })
  > 
  > ################################################################################
  > # tests for behavior of the order argument
  > e <- new.env()
  > with(e, {
  +        test <- tryCatch(btensor(x = list(spdg$age, spdg$ttm), order = c(3, 2)), error = function(e) e)
  +        test_error(test, "is.list(order) is not TRUE")
  + 
  +        test <- tryCatch(btensor(x = list(spdg$age, spdg$ttm), order = list(c(3, 2))), error = function(e) e)
  +        test_error(test, "length(order) == length(x) is not TRUE")
  + 
  +        test <- tryCatch(btensor(x = list(spdg$age, spdg$ttm), order = list(a = c(3, 2), 2)), error = function(e) e)
  +        test_error(test, "all(sapply(order, length) == 1) is not TRUE")
  + 
  +        x0 <- btensor(x = list(spdg$age, spdg$ttm),
  +                      bknots = list(c(44, 53), c(-9, -1)),
  +                      order = list(2.9, 4.1))
  + 
  +        x1 <- btensor(x = list(spdg$age, spdg$ttm),
  +                      bknots = list(c(44, 53), c(-9, -1)),
  +                      order = list(2, 4))
  + 
  +        stopifnot(isTRUE(all.equal(x0, x1, check.attributes = FALSE)))
  + })
  > 
  > ################################################################################
  > # warning thrown when both iknots and df are provided
  > e <- new.env()
  > with(e, {
  +        test <- tryCatch(btensor(x = list(spdg$age, spdg$ttm),
  +                                 iknots = list(c(46, 48), c(-8)),
  +                                 df = list(5, 5),
  +                                 bknots = list(c(44, 53), c(-9, -1)),
  +                                 order = list(3, 2)), warning = function(w) w)
  +        test_warning(test, "Both iknots and df defined, using iknots")
  + 
  +        x0 <- suppressWarnings(btensor(x = list(spdg$age, spdg$ttm),
  +                      iknots = list(c(46, 48), c(-8)),
  +                      df = list(5, 5),
  +                      bknots = list(c(44, 53), c(-9, -1)),
  +                      order = list(3, 2)))
  + 
  +        x1 <- btensor(x = list(spdg$age, spdg$ttm),
  +                      iknots = list(c(46, 48), c(-8)),
  +                      bknots = list(c(44, 53), c(-9, -1)),
  +                      order = list(3, 2))
  + 
  +        stopifnot(isTRUE(all.equal(x0, x1, check.attributes = FALSE)))
  + 
  + })
  > 
  > ################################################################################
  > # error when length of x, iknots, bknots, order are not all the same
  > e <- new.env()
  > with(e, {
  + 
  +   test <- tryCatch(btensor(x = list(spdg$age, spdg$ttm),
  +                           iknots = list(c(46, 48)),
  +                           bknots = list(c(44, 53), c(-9, -1)),
  +                           order = list(3, 2)), error = function(e) e)
  +   test_error(test, "Length of x, iknots, bknots, and order must be the same.")
  + 
  +   test <- tryCatch(btensor(x = list(spdg$age, spdg$ttm),
  +                           df = list(5),
  +                           bknots = list(c(44, 53), c(-9, -1)),
  +                           order = list(3, 2)), error = function(e) e)
  +   test_error(test, "length(df) == length(x) is not TRUE")
  + 
  + })
  > 
  > ################################################################################
  > # expected construction when given df instead of iknots
  > e <- new.env()
  > with(e, {
  + 
  +        x0 <- btensor(list(spdg$age, spdg$ttm),
  +                           df = list(5, 5),
  +                           bknots = list(c(44, 53), c(-9, -1)),
  +                           order = list(4, 3))
  +        x1 <- btensor(list(spdg$age, spdg$ttm),
  +                      iknots = list(trimmed_quantile(spdg$age, probs = 0.5),
  +                                    trimmed_quantile(spdg$ttm, probs = 1:2/3)),
  +                           #df = list(5, 5),
  +                           bknots = list(c(44, 53), c(-9, -1)),
  +                           order = list(4, 3))
  + 
  +        stopifnot(isTRUE(all.equal(x0, x1, check.attributes = FALSE)))
  + 
  + })
  > 
  > e <- new.env()
  > with(e, {
  + 
  +        x0 <- tryCatch(btensor(list(spdg$age, spdg$ttm),
  +                           df = list(5, 5),
  +                           bknots = list(c(44, 53), c(-9, -1)),
  +                           order = list(4, 6)),
  +                       warning = function(w) w)
  +        test_warning(x0, "df being set to order")
  + 
  +        x0 <- btensor(list(spdg$age, spdg$ttm),
  +                           df = list(5, 5),
  +                           bknots = list(c(44, 53), c(-9, -1)),
  +                           order = list(4, 6))
  +        x1 <- btensor(list(spdg$age, spdg$ttm),
  +                      iknots = list(trimmed_quantile(spdg$age, probs = 0.5),
  +                                    numeric(0)),
  +                           bknots = list(c(44, 53), c(-9, -1)),
  +                           order = list(4, 6))
  + 
  +        stopifnot(isTRUE(all.equal(x0, x1, check.attributes = FALSE)))
  + 
  + })
  Warning message:
  In (function (xx, dd, oo)  : df being set to order
  > 
  > e <- new.env()
  > with(e, {
  + 
  +        x0 <- btensor(list(spdg$age, spdg$ttm),
  +                      df = list(4, 5),
  +                      bknots = list(c(44, 53), c(-9, -1)),
  +                      order = list(4, 5))
  + 
  +        x1 <- btensor(list(spdg$age, spdg$ttm),
  +                      iknots = list(numeric(0), numeric(0)),
  +                      bknots = list(c(44, 53), c(-9, -1)),
  +                      order = list(4, 5))
  + 
  +        stopifnot(isTRUE(all.equal(x0, x1, check.attributes = FALSE)))
  + 
  + })
  > 
  > ################################################################################
  > ##                              printing method                               ##
  > e <- new.env()
  > with(e, {
  +   bm <-
  +     btensor(x = list(mtcars$disp, mtcars$hp),
  +             iknots = list(numeric(0), c(100, 150)),
  +             bknots = list(c(70, 475), c(50, 350))
  +     )
  + 
  +   stopifnot(identical(bm, print(bm)))
  + 
  +   bmcap <- capture.output(print(bm))
  +   stopifnot(identical(bmcap[1], "Tensor Product Matrix dims: [32 x 24]"))
  + 
  + })
  Tensor Product Matrix dims: [32 x 24]
  
  > 
  > ################################################################################
  > ##                                End of File                                 ##
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.87    0.25    1.11 
Running the tests in 'tests/test-build_tensor.R' failed.
Complete output:
  > library(cpr)
  > 
  > A <- matrix(1:4, nrow = 10, ncol = 20)
  > B <- matrix(1:6, nrow = 10, ncol = 6)
  > 
  > # verify that passing in one matrix returns that matrix
  > x <- tryCatch(build_tensor(A, 1:10), error = function(e) e)
  > stopifnot(inherits(x, "error"))
  > stopifnot(identical(x$message, "All arguments passed to build_tensor need to be matrices."))
  > 
  > # verify that passing in one matrix returns that matrix
  > tA <- build_tensor(A)
  > stopifnot(identical(tA, A))
  > 
  > # Two ways of building the same tensor product
  > tAB  <- build_tensor(A, B)
  > tAB2 <- do.call(build_tensor,list(A, B))
  > stopifnot(all.equal(tAB, tAB2))
  > 
  > # direction matters
  > tBA <- build_tensor(B, A)
  > stopifnot(!isTRUE(all.equal(tBA, tAB)))
  > 
  > # three and four matrix tensor product
  > tABA  <- build_tensor(A, B, A)
  > tABAA <- build_tensor(A, B, A, A)
  > 
  > # verify to established matrices
  > rdsAB   <- readRDS("tensorAB.rds")
  > rdsBA   <- readRDS("tensorBA.rds")
  > rdsABA  <- readRDS("tensorABA.rds")
  > rdsABAA <- readRDS("tensorABAA.rds")
  > 
  > stopifnot(identical(tAB, rdsAB))
  > stopifnot(identical(tBA, rdsBA))
  > stopifnot(identical(tABA, rdsABA))
  > stopifnot(identical(tABAA, rdsABAA))
  > 
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.28    0.10    0.36 
Running the tests in 'tests/test-cp.R' failed.
Complete output:
  > library(cpr)
  > require(lme4)
  Loading required package: lme4
  Loading required package: Matrix
  > require(geepack)
  Loading required package: geepack
  > ################################################################################
  > # Verify that cp.cpr_bs and cp.formula both return similar objects
  > e <- new.env()
  > with(e, {
  + 
  +   xvec <- runif(500, 0, 6)
  +   bknots <- c(0, 6)
  +   dat <- data.frame(x = xvec, y = sin((xvec - 2)/pi) + 1.4 * cos(xvec/pi))
  +   acp <- cp(y ~ bsplines(x, df = 8, bknots = bknots), data = dat)
  +   theta <- coef(lm(y ~ bsplines(x, df = 8, bknots = bknots) - 1, data = dat))
  +   bcp <- cp(bsplines(xvec, df =8, bknots = bknots), theta)
  + 
  +   stopifnot(isTRUE(all.equal(names(acp), names(bcp))))
  + 
  +   stopifnot(identical(
  +     names(acp)
  +     ,
  +     c("cp", "xi", "iknots", "bknots", "order", "call", "keep_fit", "fit", "theta", "coefficients", "vcov", "vcov_theta", "loglik", "rss", "rse")
  +   ))
  + 
  + })
  > 
  > ################################################################################
  > # Verify that an error is thrown if bsplines is not used as expected in the
  > # formula
  > e <- new.env()
  > with(e, {
  +   test <- tryCatch(cp(log10(pdg) ~ age + ttm, data = spdg), error = function(e) e)
  +   stopifnot(inherits(test, "error"))
  +   stopifnot(identical(test$message, "bsplines() must appear first, once, and with no effect modifiers, on the right hand side of the formula."))
  + })
  > 
  > e <- new.env()
  > with(e, {
  +   test <- tryCatch(cp(log10(pdg) ~ age + bsplines(ttm), data = spdg), error = function(e) e)
  +   stopifnot(inherits(test, "error"))
  +   stopifnot(identical(test$message, "bsplines() must appear first, once, and with no effect modifiers, on the right hand side of the formula."))
  + })
  > 
  > e <- new.env()
  > with(e, {
  +   test <- tryCatch(cp(log10(pdg) ~ bsplines(ttm)*age, data = spdg), error = function(e) e)
  +   stopifnot(inherits(test, "error"))
  +   stopifnot(identical(test$message, "bsplines() must appear first, once, and with no effect modifiers, on the right hand side of the formula."))
  + })
  > 
  > ################################################################################
  > # Verify that a control polygon can be build from a lm
  > e <- new.env()
  > with(e, {
  +   xvec <- seq(0, 5.9999, length = 500)
  +   bknots <- c(0, 6)
  +   dat <- data.frame(x = xvec, y = sin((xvec - 2)/pi) + 1.4 * cos(xvec/pi))
  +   cp3 <- cp(y ~ bsplines(x, bknots = bknots), data = dat)
  + 
  +   stopifnot(
  +     isTRUE(
  +       all.equal(
  +         cp3$cp
  +         ,
  +         structure(list(xi_star = c(0, 2, 4, 6),
  +                        theta = c(0.797026006387093, 1.36601191348564, 1.19010324873104, 0.482219646221653)),
  +                   class = "data.frame", row.names = c(NA, -4L))
  +       )
  +     )
  +   )
  + 
  +   stopifnot(
  +     isTRUE(
  +       all.equal(
  +         summary(cp3)
  +         ,
  +         structure(list(dfs = 4L, n_iknots = 0L, iknots = structure(list( numeric(0)), class = "AsIs"), loglik = 2218.47902453217, rss = 0.00409829163655395, rse = 0.0028744886068859, wiggle = structure(0.0685628056150533, abs.error = 7.61200054267312e-16, subdivisions = 1L, message = "OK"), fdsc = 1), row.names = c(NA, -1L), class = c("cpr_summary_cpr_cp", "data.frame"))
  +       )
  +     )
  +   )
  + })
  > 
  > ################################################################################
  > # verify that a control ploygon can be build via lmer
  > e <- new.env()
  > with(e, {
  + 
  +   lmer_cp <- cp(log10(pdg) ~ bsplines(day, bknots = c(-1, 1)) + (1 | id)
  +             , data = spdg
  +             , method = lmer)
  + 
  +   stopifnot(
  +     isTRUE(
  +       all.equal(
  +         lmer_cp$cp
  +         ,
  +         structure(list(xi_star = c(-1, -0.333333333333333, 0.333333333333333, 1), theta = c(0.192421846281352, -2.18938687153398, 1.98993907207642, 0.0777684276192437)), class = "data.frame", row.names = c(NA, -4L))
  +       )
  +     )
  +   )
  + 
  +   stopifnot(
  +     isTRUE(
  +       all.equal(
  +         summary(lmer_cp)
  +         ,
  +         structure(list(dfs = 4L, n_iknots = 0L, iknots = structure(list( numeric(0)), class = "AsIs"), loglik = 7909.94604102081, rss = 622.553425369456, rse = 0.159004352365617, wiggle = structure(60.2815339858114, abs.error = 6.69259469907717e-13, subdivisions = 1L, message = "OK"), fdsc = 2), row.names = c(NA, -1L), class = c("cpr_summary_cpr_cp", "data.frame"))
  +       )
  +     )
  +   )
  + 
  + })
  > 
  > ################################################################################
  > # Verify that a control polygon can be build from a gee
  > e <- new.env()
  > with(e, {
  + 
  +   gee_cp <- cp(log10(pdg) ~ bsplines(day, bknots = c(-1, 1))
  +                , data = spdg
  +                , method = geeglm
  +                , method.args = list(id = as.name("id"), corstr = "ar1")
  +   )
  + 
  +   stopifnot(
  +     isTRUE(
  +       all.equal(
  +         gee_cp$cp
  +         ,
  +         structure(list(xi_star = c(-1, -0.333333333333333, 0.333333333333333, 1), theta = c(0.0066448178731621, -1.91923731782891, 1.93439356253956, 0.0233810025421451)), class = "data.frame", row.names = c(NA, -4L))
  +       )
  +     )
  +   )
  + 
  +   stopifnot(
  +     isTRUE(
  +       all.equal(
  +         summary(gee_cp)
  +         ,
  +         structure(list(dfs = 4L, n_iknots = 0L, iknots = structure(list( numeric(0)), class = "AsIs"), loglik = -1464.93520368629, rss = 2929.87040737258, rse = 0.344941068561338, wiggle = structure(49.9755793513228, abs.error = 5.54840388648201e-13, subdivisions = 1L, message = "OK"), fdsc = 2), row.names = c(NA, -1L), class = c("cpr_summary_cpr_cp", "data.frame"))
  +       )
  +     )
  +   )
  + 
  +   gee_cp <- cp(log10(pdg) ~ bsplines(day, df = 10, bknots = c(-1, 1))
  +                , data = spdg
  +                , method = geepack::geeglm
  +                , method.args = list(id = as.name("id"), corstr = "ar1")
  +   )
  + 
  +   gee_cpr <- cpr(gee_cp)
  +   expected_summary <- structure(list(dfs = 4:10, n_iknots = 0:6, iknots = structure(list( numeric(0), -0.188465250965251, c(-0.188465250965251, 0.522392938868911), c(-0.387864823348694, -0.188465250965251, 0.522392938868911), c(-0.791871921182266, -0.387864823348694, -0.188465250965251, 0.522392938868911), c(-0.791871921182266, -0.591390091390091, -0.387864823348694, -0.188465250965251, 0.522392938868911), c(-0.791871921182266, -0.591390091390091, -0.387864823348694, -0.188465250965251, 0.0576230492196879, 0.522392938868911)), class = "AsIs"), loglik = c(-1464.93520368629, -1376.00702167216, -1366.38011307588, -1362.38424749121, -1362.23759222167, -1362.19963612715, -1362.29607387444), rss = c(2929.87040737258, 2752.01404334433, 2732.76022615176, 2724.76849498243, 2724.47518444335, 2724.39927225431, 2724.59214774888), rse = c(0.344941068561338, 0.334314212839968, 0.333149449966508, 0.332668714884735, 0.332657564802185, 0.332659686294428, 0.33267821811044), wiggle = structure(c(49.9755793513228, 57.2632257630784, 36.42646977805, 39.72463162264, 42.8679143932625, 40.3142809497938, 40.19072967787), abs.error = 5.54840388648201e-13, subdivisions = 1L, message = "OK"), fdsc = c(2, 2, 2, 2, 2, 2, 2), `Pr(>w_(1))` = c(NA, 0, 0, 0, 9.01352735449557e-07, 0.00711887714037818, 0.0103830482941717)), row.names = c(NA, -7L), class = c("cpr_summary_cpr_cpr", "cpr_summary_cpr_cp", "data.frame"), elbow = structure(c(3, 2, 3, 2, 3, 2), dim = 2:3, dimnames = list(c("quadratic", "linear"), c("loglik", "rss", "rse"))))
  +   stopifnot(isTRUE(all.equal(expected_summary, summary(gee_cpr))))
  + 
  + })
  
    |                                                                            
    |                                                                      |   0%
    |                                                                            
    |==========                                                            |  14%
    |                                                                            
    |====================                                                  |  29%
    |                                                                            
    |==============================                                        |  43%
    |                                                                            
    |========================================                              |  57%
    |                                                                            
    |==================================================                    |  71%
    |                                                                            
    |============================================================          |  86%
    |                                                                            
    |======================================================================| 100%
  > 
  > ################################################################################
  > #                                rank deficient?                               #
  > e <- new.env()
  > with(e, {
  + 
  +   # First, a good fit
  +   cp0 <- cp(pdg ~ bsplines(day, bknots = c(-1, 1)) + ttm , data = spdg)
  + 
  +   stopifnot(inherits(cp0, "cpr_cp"))
  + 
  +   # Now update to something that is rank deficient
  +   cp1 <- tryCatch(update_bsplines(cp0, iknots = c(0, 0, 0, 0, 0)), warning = function(w) w)
  + 
  +   stopifnot(inherits(cp1, 'warning'))
  +   stopifnot(identical(cp1$message, 'Design Matrix is rank deficient. keep_fit being set to TRUE.'))
  + 
  +   cp2 <- tryCatch(cp(pdg ~ bsplines(day, iknots = c(0, 0, 0, 0, 0), bknots = c(-1, 1)) + ttm , data = spdg), warning = function(w) w)
  + 
  +   stopifnot(inherits(cp2, 'warning'))
  +   stopifnot(identical(cp2$message, 'Design Matrix is rank deficient. keep_fit being set to TRUE.'))
  + 
  + })
  > 
  > ################################################################################
  > ##                              printing method                               ##
  > e <- new.env()
  > with(e, {
  +   xvec <- runif(500, 0, 6)
  +   bknots <- c(0, 6)
  +   dat <- data.frame(x = xvec, y = sin((xvec - 2)/pi) + 1.4 * cos(xvec/pi))
  +   acp <- cp(y ~ bsplines(x, df = 8, bknots = bknots), data = dat)
  +   theta <- coef(lm(y ~ bsplines(x, df = 8, bknots = bknots) - 1, data = dat))
  +   bcp <- cp(bsplines(xvec, df =8, bknots = bknots), theta)
  + 
  +   # verify the value is returned from the print call
  +   stopifnot(identical(acp, print(acp)))
  +   stopifnot(identical(bcp, print(bcp)))
  + 
  +   acpcap <- capture.output(print(acp))
  +   expected <- capture.output(print(acp$cp))
  +   stopifnot(identical(acpcap, expected))
  + 
  +   bcpcap <- capture.output(print(bcp))
  +   expected <- capture.output(print(bcp$cp))
  +   stopifnot(identical(bcpcap, expected))
  + 
  + })
      xi_star     theta
  1 0.0000000 0.8055015
  2 0.4466616 0.9199736
  3 1.2760577 1.0859647
  4 2.4694993 1.1625714
  5 3.6220090 1.0861496
  6 4.7926129 0.8607775
  7 5.5991713 0.6214524
  8 6.0000000 0.4904246
      xi_star     theta
  1 0.0000000 0.8055015
  2 0.4466616 0.9199736
  3 1.2760577 1.0859647
  4 2.4694993 1.1625714
  5 3.6220090 1.0861496
  6 4.7926129 0.8607775
  7 5.5991713 0.6214524
  8 6.0000000 0.4904246
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
    32.07    0.50   32.59 
Running the tests in 'tests/test-cp_diff.R' failed.
Complete output:
  > library(cpr)
  > 
  > xvec <- runif(n = 500, min = 0, max = 6)
  > 
  > # Define the basis matrix
  > bmat1 <- bsplines(x = xvec, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
  > bmat2 <- bsplines(x = xvec, bknots = c(0, 6))
  > 
  > # Define the control vertices ordinates
  > theta1 <- c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5)
  > theta2 <- c(1, 3.4, -2, 1.7)
  > 
  > # build the two control polygons
  > cp1 <- cp(bmat1, theta1)
  > cp2 <- cp(bmat2, theta2)
  > 
  > stopifnot(
  +   isTRUE(
  +     all.equal(
  +       cp_diff(cp1, cp2)
  +       ,
  +       c(0, 1.4, -1.5, -1.28, -1.92, -0.42, 0.241666666666666, -1.225,  0.2)
  +     )
  +   )
  + )
  > 
  > 
  > ################################################################################
  > #                                 End of file                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.31    0.01    0.32 
Running the tests in 'tests/test-cp_value.R' failed.
Complete output:
  > library(cpr)
  > 
  > xvec <- runif(n = 500, min = 0, max = 6)
  > bmat1 <- bsplines(x = xvec, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
  > theta1 <- c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5)
  > cp1 <- cp(bmat1, theta1)
  > 
  > x <- c(0.2, 0.8, 1.3, 1.73, 2.15, 3.14, 4.22, 4.88, 5.3, 5.9)
  > 
  > stopifnot(
  +   isTRUE(
  +     all.equal(
  +       cp_value(cp1, x = x)
  +       ,
  +       c(0.4, 3.26666666666667, 3.92608695652174, 4.135, 3.925, 1.432, -0.600540540540541, -0.510999999999999, 1.19, 1.6)
  +     )
  +   )
  + )
  > 
  > 
  > bmat2 <- bsplines(x = xvec, bknots = c(0, 6))
  > theta2 <- c(1, 3.4, -2, 1.7)
  > cp2 <- cp(bmat2, theta2)
  > 
  > stopifnot(
  +   isTRUE(
  +     all.equal(
  +       cp_diff(cp1, cp2)
  +       ,
  +       c(0, 1.4, -1.5, -1.28, -1.92, -0.42, 0.241666666666666, -1.225, 0.2)
  +     )
  +   )
  + )
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.26    0.06    0.31 
Running the tests in 'tests/test-generate_cp_formula_data.R' failed.
Complete output:
  > library(cpr)
  > 
  > ################################################################################
  > # vertify generate_cp_formula_data is in the namespace but not exported
  > stopifnot(!grepl("generate_cp_formula_data", ls("package:cpr")))
  > 
  > cpr_namespace <- ls( getNamespace("cpr"), all.names = TRUE)
  > stopifnot(sum(cpr_namespace == "generate_cp_formula_data") == 1L)
  > stopifnot(sum(grepl("^generate_cp_formula_data", cpr_namespace)) == 1L)
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   data <-
  +     data.frame(
  +                  x1 = runif(20)
  +                , x2 = runif(20)
  +                , x3 = runif(20)
  +                , xf = factor(rep(c("l1","l2","l3","l4"), each = 5))
  +                , xc = rep(c("c1","c2","c3","c4", "c5"), each = 4)
  +                , pid = gl(n = 2, k = 10)
  +                , pid2 = rep(1:2, each = 10)
  +     )
  + 
  +   f <- ~ bsplines(x1, bknots = c(0,1)) + x2 + xf + xc + (x3 | pid2)
  + 
  +   cpr:::generate_cp_formula_data(f, data)
  + 
  +   stopifnot(isTRUE(
  +     all.equal(
  +               f_for_use
  +               ,
  +               . ~ bsplines(x1, bknots = c(0, 1)) + x2 + (x3 | pid2) + xfl2 + xfl3 + xfl4 + xcc2 + xcc3 + xcc4 + xcc5 - 1
  +               )
  +   ))
  + 
  +   stopifnot(isTRUE(identical(
  +     names(data_for_use)
  +     ,
  +     c("x1", "x2", "x3", "pid", "pid2", "xfl2", "xfl3", "xfl4", "xcc2", "xcc3", "xcc4", "xcc5")
  +   )))
  + 
  + })
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.32    0.07    0.39 
Running the tests in 'tests/test-get_spline.R' failed.
Complete output:
  > library(cpr)
  > 
  > x1 <- seq(0, 5.9999, length = 500)
  > x2 <- seq(0, 6 - sqrt(.Machine$double.eps), length = 123)
  > bmat1 <- bsplines(x = x1, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
  > bmat2 <- bsplines(x = x2, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
  > theta <- c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5)
  > 
  > cp1 <- cp(bmat1, theta)
  > cp2 <- cp(bmat2, theta)
  > 
  > spline <- get_spline(cp1, n = 123)
  > 
  > stopifnot(isTRUE(
  +   all.equal(
  +             spline
  +             ,
  +             data.frame(x = x2, y = as.numeric(bmat2 %*% theta))
  +   )
  + ))
  > 
  > # control net
  > initial_cn <- cn(log10(pdg) ~ btensor(list(day, age, ttm)
  +                         , df = list(8, 6, 3)
  +                         , bknots = list(c(-1, 1), c(44, 53), c(-9, -1))
  +                         , order = list(3, 4, 2)
  +                         )
  +           , data = spdg)
  > 
  > if (interactive()) {
  +   old_par <- par()
  +   par(mfrow = c(3, 1))
  +   get_spline(initial_cn, n = 100) |> plot(type = "l")
  +   get_spline(initial_cn, margin = 2, n = 100) |> plot(type = "l")
  +   get_spline(initial_cn, margin = 3, n = 100) |> plot(type = "l")
  +   par(old_par)
  + 
  +   #get_spline(initial_cn, n = 10) |> dput()
  +   #get_spline(initial_cn, margin = 2, n = 10) |> dput()
  +   #get_spline(initial_cn, margin = 3, n = 10) |> dput()
  + }
  > 
  > expected <- structure(list(x = c(-1, -0.777777779433462, -0.555555558866925, -0.333333338300387, -0.111111117733849, 0.111111102832688, 0.333333323399226, 0.555555543965764, 0.777777764532301, 0.999999985098839), y = c(-0.133762887106663, -0.320552726835096, -0.453930169776096, -0.461840210183064, -0.332093945041971, 0.139354295999598, 0.496456762997947, 0.686489715816332, 0.534159573348981, 0.0146225055203891)), class = "data.frame", row.names = c(NA, -10L))
  > stopifnot(isTRUE(all.equal(expected, get_spline(initial_cn, n = 10))))
  > 
  > expected <- structure(list(x = c(44, 44.9999999983443, 45.9999999966886, 46.9999999950329, 47.9999999933773, 48.9999999917216, 49.9999999900659, 50.9999999884102, 51.9999999867545, 52.9999999850988), y = c(0.686285778255235, 0.00445867954416262, -0.196819616993946, -0.14570083255892, -0.0703366883505854, -0.139325216264028, -0.22418621303892, -0.188094354192221, -0.0451158549211667, 0.175094646687039)), class = "data.frame", row.names = c(NA, -10L))
  > stopifnot(isTRUE(all.equal(expected, get_spline(initial_cn, margin = 2, n = 10))))
  > 
  > expected <- structure(list(x = c(-9, -8.1111111127668, -7.22222222553359, -6.33333333830039, -5.44444445106718, -4.55555556383398, -3.66666667660077, -2.77777778936757, -1.88888890213437, -1.00000001490116), y = c(-0.121538778218393, -0.114302628952598, -0.107066479686803, -0.0998303304210085, -0.0925941811552136, -0.0864796293184603, -0.0812970496621639, -0.0761144700058674, -0.070931890349571, -0.0657493106932745)), class = "data.frame", row.names = c(NA, -10L))
  > stopifnot(isTRUE(all.equal(expected, get_spline(initial_cn, margin = 3, n = 10) )))
  > 
  > ################################################################################
  > # expected an error if length(margin) > 1
  > 
  > test <- tryCatch(get_spline(initial_cn, margin = 2:3, n = 100), error = function(e) e)
  > stopifnot(inherits(test, "error"))
  > stopifnot(identical(test$message, "use get_surface when length(margin) > 1."))
  > 
  > ################################################################################
  > # expect an error if derivative is not 0, 1, or 2
  > test <- tryCatch(get_spline(cp1, derivative = 0:1), error = function(e) e)
  > stopifnot(inherits(test, "error"))
  > stopifnot(identical(test$message, "length(derivative) == 1L is not TRUE"))
  > 
  > test <- tryCatch(get_spline(cp1, derivative = 3), error = function(e) e)
  > stopifnot(inherits(test, "error"))
  > stopifnot(identical(test$message, "derivative needs to be 0, 1, or 2"))
  > 
  > ################################################################################
  > # derivatives are as expected
  > f0 <- function(x) {
  +   #(x + 2) * (x - 1) * (x - 3)
  +   x^3 - 2 * x^2 - 5 * x + 6
  + }
  > f1 <- function(x) {
  +   3 * x^2 - 4 * x - 5
  + }
  > f2 <- function(x) {
  +   6 * x - 4
  + }
  > 
  > x <- sort(runif(n = 100, min = -3, max = 5))
  > bknots = c(-3, 5)
  > bmat <- bsplines(x, bknots = bknots)
  > theta <- coef(lm(f0(x) ~ bsplines(x, bknots = bknots) + 0) )
  > 
  > cp0 <- cp(bmat, theta)
  > spline0 <- get_spline(cp0, derivative = 0)
  > spline1 <- get_spline(cp0, derivative = 1)
  > spline2 <- get_spline(cp0, derivative = 2)
  > 
  > stopifnot(isTRUE(all.equal(f0(spline0$x), spline0$y)))
  > stopifnot(isTRUE(all.equal(f1(spline1$x), spline1$y)))
  > stopifnot(isTRUE(all.equal(f2(spline2$x), spline2$y)))
  > 
  > ################################################################################
  > # Standard errors
  > test <- tryCatch(get_spline(cp0, se = TRUE), warning = function(w) w)
  > stopifnot(inherits(test, "warning"))
  > stopifnot(identical(test$message, "vcov_theta of cp0 is NULL"))
  > 
  > set.seed(742)
  > x <- sort(runif(n = 2222, min = -3, max = 5))
  > sdata <- data.frame(x = x, y = f0(x) + rnorm(2222, sd = 25))
  > cp0 <- cp(y ~ bsplines(x, bknots = bknots), data = sdata)
  > spline0 <- get_spline(cp0, se = TRUE, derivative = 0, n = 20)
  > 
  > stopifnot(
  +   isTRUE(
  +     all.equal(
  +       spline0$se
  +       ,
  +       predict(cp0, newdata = data.frame(x = spline0$x), se.fit = TRUE)$se.fit
  +       ,
  +       check.attributes = FALSE
  +     )
  +   )
  + )
  > 
  > spline1 <- get_spline(cp0, se = TRUE, derivative = 1, n = 20)
  > spline2 <- get_spline(cp0, se = TRUE, derivative = 2, n = 20)
  > 
  > if (interactive()) {
  +   old_par <- par()
  + 
  +   par(mfrow = c(1, 3))
  + 
  +   plot(x, f0(x), type = "l", main = "spline")
  +   lines(spline0$x, spline0$y, pch = 2, col = 'blue')
  +   lines(spline0$x, spline0$y + 1.96 * spline0$se, pch = 1, col = 'red')
  +   lines(spline0$x, spline0$y - 1.96 * spline0$se, pch = 1, col = 'green')
  + 
  +   plot(x, f1(x), type = "l", main = "first derivative")
  +   lines(spline1$x, spline1$y, pch = 2, col = 'blue')
  +   lines(spline1$x, spline1$y + 1.96 * spline1$se, pch = 1, col = 'red')
  +   lines(spline1$x, spline1$y - 1.96 * spline1$se, pch = 1, col = 'green')
  + 
  +   plot(x, f2(x), type = "l", main = "second derivative")
  +   lines(spline2$x, spline2$y, pch = 2, col = 'blue')
  +   lines(spline2$x, spline2$y + 1.96 * spline2$se, pch = 1, col = 'red')
  +   lines(spline2$x, spline2$y - 1.96 * spline2$se, pch = 1, col = 'green')
  + 
  +   par(old_par)
  + }
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     3.31    0.79    4.04 
Running the tests in 'tests/test-influence_of_iknots.R' failed.
Complete output:
  > library(cpr)
  > set.seed(42)
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   x <- runif(n = 100, 0, 6)#seq(0 + 1/5000, 6 - 1/5000, length.out = 100)
  +   bmat <- bsplines(x, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
  +   theta <- matrix(c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5), ncol = 1)
  +   DF <- data.frame(x = x, truth = as.numeric(bmat %*% theta))
  +   DF$y <- as.numeric(bmat %*% theta + rnorm(nrow(bmat), sd = 0.3))
  + 
  +   initial_cp <-
  +     cp(y ~ bsplines(x, iknots = c(1, 1.5, 2.3, 3.0, 4, 4.5), bknots = c(0, 6)), data = DF)
  + 
  +   iok <- influence_of_iknots(initial_cp)
  + 
  +   stopifnot(inherits(iok, "cpr_influence_of_iknots"))
  + 
  +   # if the following test for the names of the object changes make sure to
  +   # update the @return section of the source file.
  +   stopifnot(identical(names(iok), c("original_cp", "coarsened_cps", "restored_cps", "d", "influence", "chisq")))
  + 
  +   s <- summary(iok)
  + 
  +   stopifnot(identical(dim(s), c(6L, 8L)))
  +   stopifnot(identical(names(s), c("j", "iknot", "influence", "influence_rank", "chisq", "chisq_rank", "p_value", "os_p_value")))
  + })
  > 
  > ################################################################################
  > # verbose vs non-versbose
  > e <- new.env()
  > with(e, {
  +   x <- runif(n = 100, 0, 6)
  +   bmat <- bsplines(x, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
  +   theta <- matrix(c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5), ncol = 1)
  +   DF <- data.frame(x = x, truth = as.numeric(bmat %*% theta))
  +   DF$y <- as.numeric(bmat %*% theta + rnorm(nrow(bmat), sd = 0.3))
  + 
  +   initial_cp <-
  +     cp(y ~ bsplines(x, iknots = c(1, 1.5, 2.3, 3.0, 4, 4.5), bknots = c(0, 6)), data = DF)
  + 
  +   iok0 <- influence_of_iknots(initial_cp, verbose = FALSE, cl = 0)
  +   iok1 <- influence_of_iknots(initial_cp, verbose = FALSE, cl = 1)
  +   iok2 <- influence_of_iknots(initial_cp, verbose = FALSE, cl = 2)
  +   iok0v <- influence_of_iknots(initial_cp, verbose = TRUE, cl = 0)
  +   iok1v <- influence_of_iknots(initial_cp, verbose = TRUE, cl = 1)
  +   iok2v <- influence_of_iknots(initial_cp, verbose = TRUE, cl = 2)
  + 
  +   stopifnot(isTRUE(all.equal(iok0, iok1)))
  +   stopifnot(isTRUE(all.equal(iok0, iok2)))
  +   stopifnot(isTRUE(all.equal(iok0, iok0v)))
  +   stopifnot(isTRUE(all.equal(iok0, iok1v)))
  +   stopifnot(isTRUE(all.equal(iok0, iok2v)))
  + 
  + })
  
  There are 6 internal knots to evaluate
  
    Coarsening theta (step 1 of 6)
    getting meta data for coarsend basis matrices (step 2 of 6)
    generating coarsened control polygons (step 3 of 6)
    building theta hat (step 4 of 6)
    building restored control polygons (step 5 of 6)
    building test statistics (step 6 of 6)
  
  There are 6 internal knots to evaluate
  
    Coarsening theta (step 1 of 6)
    getting meta data for coarsend basis matrices (step 2 of 6)
    generating coarsened control polygons (step 3 of 6)
    building theta hat (step 4 of 6)
    building restored control polygons (step 5 of 6)
    building test statistics (step 6 of 6)
  
  There are 6 internal knots to evaluate
  
    Coarsening theta (step 1 of 6)
    getting meta data for coarsend basis matrices (step 2 of 6)
    generating coarsened control polygons (step 3 of 6)
    building theta hat (step 4 of 6)
    building restored control polygons (step 5 of 6)
    building test statistics (step 6 of 6)
  > 
  > ################################################################################
  > # CPR
  > # verbose vs non-versbose
  > e <- new.env()
  > with(e, {
  +   x <- runif(n = 100, 0, 6)
  +   bmat <- bsplines(x, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
  +   theta <- matrix(c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5), ncol = 1)
  +   DF <- data.frame(x = x, truth = as.numeric(bmat %*% theta))
  +   DF$y <- as.numeric(bmat %*% theta + rnorm(nrow(bmat), sd = 0.3))
  + 
  +   initial_cp <-
  +     cp(y ~ bsplines(x, iknots = c(1, 1.5, 2.3, 3.0, 4, 4.5), bknots = c(0, 6)), data = DF)
  +   cpr_run <- cpr(initial_cp)
  + 
  +   ioik_cpr <- influence_of_iknots(cpr_run)
  + 
  +   ioik_list <- lapply(cpr_run, influence_of_iknots)
  + 
  +   stopifnot(identical(ioik_cpr, ioik_cpr))
  + 
  +   stopifnot(
  +     identical(
  +       attr(cpr_run, "ioik")
  +       ,
  +       lapply(ioik_list, summary)
  +     )
  +   )
  + 
  +   stopifnot(
  +     all.equal(
  +       do.call(rbind, attr(cpr_run, "ioik"))
  +       ,
  +       summary(ioik_cpr)[, 1:8]
  +       ,
  +       check.attributes = FALSE
  +     )
  +   )
  + 
  +   stopifnot(
  +     identical(
  +       summary(ioik_cpr)[, "index"]
  +       ,
  +       c(1, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7)
  +     )
  +   )
  + 
  + })
  
    |                                                                            
    |                                                                      |   0%
    |                                                                            
    |==========                                                            |  14%
    |                                                                            
    |====================                                                  |  29%
    |                                                                            
    |==============================                                        |  43%
    |                                                                            
    |========================================                              |  57%
    |                                                                            
    |==================================================                    |  71%
    |                                                                            
    |============================================================          |  86%
    |                                                                            
    |======================================================================| 100%
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  + 
  + #  acn <- cn(pdg ~ btensor(list(day, age)
  + #                          , df = list(10, 8)
  + #                          , bknots = list(c(-1, 1), c(44, 53))
  + #                          ) + ttm
  + #            , data = spdg)
  + #
  + #  str(acn, max.level = 1)
  + #
  + #  influence_of_iknots(acn)
  + #
  + #  acnr <- cnr(acn)
  + #  stop("TEST NEEDS TO BE WRITTEN")
  + 
  + 
  + 
  + })
  NULL
  > 
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     1.92    0.17    2.09 
Running the tests in 'tests/test-loglikelihood.R' failed.
Complete output:
  > # Tests for extracting coefficients and vcov matrix from regression fits
  > library(cpr)
  > require(lme4)
  Loading required package: lme4
  Loading required package: Matrix
  > 
  > ################################################################################
  > # There is one method of interest, and it is non exported.  There are several S3
  > # methods to check.
  > stopifnot(!grepl("loglikelihood", ls("package:cpr")))
  > cpr_namespace <- ls( getNamespace("cpr"), all.names = TRUE)
  > stopifnot(sum(cpr_namespace == "loglikelihood") == 1L)
  > stopifnot(sum(cpr_namespace == "loglikelihood.default") == 1L)
  > stopifnot(sum(cpr_namespace == "loglikelihood.geeglm") == 1L)
  > stopifnot(sum(grepl("^loglikelihood", cpr_namespace)) == 3L)
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   fit <- lm(mpg ~ wt, data = mtcars)
  +   stopifnot("recover loglikelihood from lm" =
  +             isTRUE(all.equal(cpr:::loglikelihood(fit), as.numeric(stats::logLik(fit)))))
  + })
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   fit <- glm(I(mpg < 20) ~ wt, data = mtcars, family = quasibinomial())
  +   stopifnot("recover quasiloglikelihood from glm" =
  +             isTRUE(all.equal(cpr:::loglikelihood(fit), as.numeric(stats::logLik(fit)))))
  + })
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   require(lme4)
  +   fit <- lmer(mpg ~ wt + (1|cyl), data = mtcars)
  +   stopifnot("recover loglikelihood from lme4" =
  +             isTRUE(all.equal(cpr:::loglikelihood(fit), as.numeric(stats::logLik(fit)))))
  + })
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   require(geepack)
  +   fit <- geeglm(mpg ~ wt, id = cyl, data = mtcars, family = gaussian())
  +   y <- fit$y
  +   mu <- fit$fitted.values
  +   stopifnot("recover loglikelihood from gaussian geeglm" =
  +             isTRUE(all.equal(cpr:::loglikelihood(fit), -sum((y - mu)^2 / 2))))
  + })
  Loading required package: geepack
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   require(geepack)
  +   fit <- geeglm(I(mpg < 20) ~ wt, id = cyl, data = mtcars, family = binomial())
  +   y <- fit$y
  +   mu <- fit$fitted.values
  +   stopifnot("recover loglikelihood from binomial geeglm" =
  +             isTRUE(all.equal(cpr:::loglikelihood(fit),
  +                              sum(y * log(mu/(1 - mu)) + log(1 - mu)))))
  + })
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   require(geepack)
  +   fit <- geeglm(as.integer(mpg) ~ wt, id = cyl, data = mtcars, family = poisson())
  +   y <- fit$y
  +   mu <- fit$fitted.values
  +   stopifnot("recover loglikelihood from poisson geeglm" =
  +             isTRUE(all.equal(cpr:::loglikelihood(fit),
  +                              sum(y * log(mu) - mu))))
  + })
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   require(geepack)
  +   fit <- geeglm(mpg ~ wt, id = cyl, data = mtcars, family = Gamma())
  +   y <- fit$y
  +   mu <- fit$fitted.values
  +   stopifnot("recover loglikelihood from Gamma geeglm" =
  +             isTRUE(all.equal(cpr:::loglikelihood(fit), -sum(y/mu + log(mu)))))
  + })
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     2.18    0.34    2.53 
Running the tests in 'tests/test-modify-formula-and-data.R' failed.
Complete output:
  > # what do I need?
  > #
  > # I need a way to set no intercept without a factor variable causing a rank
  > # deficient matrix
  > 
  > library(cpr)
  > 
  > ################################################################################
  > # There is one method of interest, and it is non exported.  There are several S3
  > # methods to check.
  > stopifnot(!grepl("generate_cp_formula_data", ls("package:cpr")))
  > cpr_namespace <- ls( getNamespace("cpr"), all.names = TRUE)
  > stopifnot(sum(cpr_namespace == "generate_cp_formula_data") == 1L)
  > stopifnot(sum(grepl("generate_cp_formula_data", cpr_namespace)) == 1L)
  > 
  > ################################################################################
  > # test formula and data construction
  > 
  > e <- new.env()
  > with(e, {
  +   data <-
  +     data.frame(
  +                  x1 = runif(20)
  +                , x2 = runif(20)
  +                , x3 = runif(20)
  +                , xf = factor(rep(c("l1","l2","l3","l4"), each = 5))
  +                , xc = rep(c("c1","c2","c3","c4", "c5"), each = 4)
  +                , pid = gl(n = 2, k = 10)
  +                , pid2 = rep(1:2, each = 10)
  +     )
  + 
  +   f <- ~ bsplines(x1, bknots = c(0,1)) + x2 + xf + xc + (x3 | pid2)
  + 
  +   cpr:::generate_cp_formula_data(f, data)
  + 
  +   stopifnot(isTRUE(
  +     all.equal(
  +               f_for_use
  +               ,
  +               . ~ bsplines(x1, bknots = c(0, 1)) + x2 + (x3 | pid2) + xfl2 + xfl3 + xfl4 + xcc2 + xcc3 + xcc4 + xcc5 - 1
  +               )
  +   ))
  + 
  +   stopifnot(isTRUE(identical(
  +     names(data_for_use)
  +     ,
  +     c("x1", "x2", "x3", "pid", "pid2", "xfl2", "xfl3", "xfl4", "xcc2", "xcc3", "xcc4", "xcc5")
  +   )))
  + 
  + })
  > 
  > e <- new.env()
  > with(e, {
  +   data <-
  +     data.frame(
  +                  x1 = runif(20)
  +                , x2 = runif(20)
  +                , x3 = runif(20)
  +                , xf = factor(rep(c("l1","l2","l3","l4"), each = 5))
  +                , xc = rep(c("c1","c2","c3","c4", "c5"), each = 4)
  +                , pid = gl(n = 2, k = 10)
  +                , pid2 = rep(1:2, each = 10)
  +     )
  + 
  +   f <- ~ bsplines(x1, bknots = c(0,1)) + x2 + xf + xc
  + 
  +   cpr:::generate_cp_formula_data(f, data)
  + 
  +   stopifnot(isTRUE(
  +     all.equal(
  +               f_for_use
  +               ,
  +               . ~ bsplines(x1, bknots = c(0, 1)) + x2 + xfl2 + xfl3 + xfl4 + xcc2 + xcc3 + xcc4 + xcc5 - 1
  +               )
  +   ))
  + 
  +   stopifnot(isTRUE(identical(
  +     names(data_for_use)
  +     ,
  +     c("x1", "x2", "x3", "pid", "pid2", "xfl2", "xfl3", "xfl4", "xcc2", "xcc3", "xcc4", "xcc5")
  +   )))
  + 
  + })
  > 
  > 
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.21    0.12    0.32 
Running the tests in 'tests/test-newknots.R' failed.
Complete output:
  > library(cpr)
  > 
  > # testing of the newknots method
  > 
  > ################################################################################
  > # check that the method is not exported
  > stopifnot(!grepl("newknots", ls("package:cpr")))
  > cpr_namespace <- ls( getNamespace("cpr"), all.names = TRUE)
  > stopifnot(sum(cpr_namespace == "newknots") == 1L)
  > stopifnot(sum(grepl("^newknots$", cpr_namespace)) == 1L)
  > 
  > ################################################################################
  > 
  > cp0 <- cp(log(pdg) ~ bsplines(day, iknots = c(-.25, 0, 0.25), bknots = c(-1, 1)), data = spdg)
  > 
  > expected <- log(pdg) ~ bsplines(day, iknots = c(-0.85, 0, 0.25, 0.3), bknots = c(-1, 1))
  > 
  > test <- cpr:::newknots(cp0$call$formula, c(-0.85, 0, 0.25, 0.3))
  > 
  > stopifnot(isTRUE( all.equal(expected, test)))
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     1.57    0.25    1.81 
Running the tests in 'tests/test-order_statistics.R' failed.
Complete output:
  > library(cpr)
  > 
  > 
  > ################################################################################
  > # Simple tests - similar to Example 1 for order_statistics
  > 
  > e <- new.env()
  > with(e, {
  +   simulated_data <- matrix(rnorm(n = 54 * 15000), ncol = 54)
  + 
  +   # find all the minimums for each of the simulated samples of size 54
  +   mins <- apply(simulated_data, 1, min)
  + 
  +   # get the density values
  +   x <- seq(-5, 0, length.out = ncol(simulated_data))
  +   d <- d_order_statistic(x, n = 54, j = 1, distribution = "norm")
  +   p <- p_order_statistic(q = x,  n = 54, j = 1, distribution = "norm")
  + 
  +   stopifnot(inherits(d, "numeric"))
  +   stopifnot(length(d) == length(x))
  +   stopifnot(inherits(p, "numeric"))
  +   stopifnot(length(p) == length(x))
  + 
  +   # plot the histogram and density -- This part of the example and should be
  +   # visually a good match
  +   if (interactive()) {
  +     par(mfrow = c(1, 2))
  + 
  +     hist(mins, freq = FALSE, breaks = 100)
  +     points(x, d, type = "l", col = "red")
  + 
  +     plot(ecdf(mins))
  +     points(x, p, col = "red")
  +   }
  + })
  > rm(e)
  > 
  > ################################################################################
  > # A NA is returned when in the j or x,q arguments
  > e <- new.env()
  > with(e, {
  +   d <- d_order_statistic(-3, n = 5, j = c(1, NA, 2), distribution = "norm")
  +   p <- p_order_statistic(-3, n = 5, j = c(1, NA, 2), distribution = "norm")
  +   stopifnot( is.na(d) == c(FALSE, TRUE, FALSE) )
  +   stopifnot( is.na(p) == c(FALSE, TRUE, FALSE) )
  + 
  +   d <- d_order_statistic(c(-3, 2, NA), n = 5, j = c(1, NA, 2), distribution = "norm")
  +   p <- p_order_statistic(c(-3, 2, NA), n = 5, j = c(1, NA, 2), distribution = "norm")
  +   stopifnot( is.na(d) == c(FALSE, TRUE, TRUE) )
  +   stopifnot( is.na(p) == c(FALSE, TRUE, TRUE) )
  + })
  > rm(e)
  > 
  > 
  > ################################################################################
  > # j-order statistics
  > # if x or q, and j have the same length, then a return
  > 
  > e <- new.env()
  > with(e, {
  +   x <- j <- c(1, 2)
  +   d <- d_order_statistic(x = x, n = 2, j = j, distribution = "norm")
  +   p <- p_order_statistic(q = x, n = 2, j = j, distribution = "norm")
  +   stopifnot(inherits(d, "numeric"))
  +   stopifnot(length(d) == length(x))
  +   stopifnot(inherits(p, "numeric"))
  +   stopifnot(length(p) == length(x))
  + })
  > rm(e)
  > 
  > ################################################################################
  > # verify errors are thrown when inputs are of unequal length or otherwise not as
  > # expected
  > 
  > e <- new.env()
  > with(e, {
  +   d <- tryCatch(d_order_statistic(x = x, n = 4:6, j = 2, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(d, "simpleError"))
  +   stopifnot(d$message == "length(n) == 1 is not TRUE")
  +   p <- tryCatch(p_order_statistic(q = x, n = 4:6, j = 2, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(p, "simpleError"))
  +   stopifnot(p$message == "length(n) == 1 is not TRUE")
  + })
  > rm(e)
  > 
  > e <- new.env()
  > with(e, {
  +   d <- tryCatch(d_order_statistic(x = x, n = numeric(0), j = 2, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(d, "simpleError"))
  +   stopifnot(d$message == "length(n) == 1 is not TRUE")
  + 
  +   p <- tryCatch(p_order_statistic(q = x, n = numeric(0), j = 2, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(p, "simpleError"))
  +   stopifnot(p$message == "length(n) == 1 is not TRUE")
  + })
  > rm(e)
  > 
  > e <- new.env()
  > with(e, {
  +   d <- tryCatch(d_order_statistic(x = 0, n = NA_real_, j = 2, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(d, "simpleError"))
  +   stopifnot(d$message == "!is.na(n) is not TRUE")
  + 
  +   p <- tryCatch(p_order_statistic(q = 0, n = NA_real_, j = 2, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(p, "simpleError"))
  +   stopifnot(p$message == "!is.na(n) is not TRUE")
  + })
  > rm(e)
  > 
  > e <- new.env()
  > with(e, {
  +   d <- tryCatch(d_order_statistic(x = 0, n = 10, j = 11, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(d, "simpleError"))
  +   stopifnot(d$message == "n >= stats::na.omit(j) is not TRUE")
  + 
  +   p <- tryCatch(p_order_statistic(q = 0, n = 10, j = 11, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(p, "simpleError"))
  +   stopifnot(p$message == "n >= stats::na.omit(j) is not TRUE")
  + })
  > rm(e)
  > 
  > e <- new.env()
  > with(e, {
  +   d <- tryCatch(d_order_statistic(x = 0, n = 10, j = -1, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(d, "simpleError"))
  +   stopifnot(d$message == "stats::na.omit(j) >= 1 is not TRUE")
  + 
  +   p <- tryCatch(p_order_statistic(q = 0, n = 10, j = -1, distribution = "norm")
  +                 , error = function(e) e)
  +   stopifnot(inherits(p, "simpleError"))
  +   stopifnot(p$message == "stats::na.omit(j) >= 1 is not TRUE")
  + })
  > rm(e)
  > 
  > ################################################################################
  > e <- new.env()
  > with(e, {
  +   stopifnot(
  +     identical(
  +       p_order_statistic(numeric(0), n = 10, j = 3, distribution = "norm")
  +       , numeric(0)
  +       )
  +     )
  + })
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.34    0.09    0.42 
Running the tests in 'tests/test-plot.cpr_bs.R' failed.
Complete output:
  > 
  > proc.time()
     user  system elapsed 
     0.17    0.06    0.20 
Running the tests in 'tests/test-plot.cpr_cn.R' failed.
Complete output:
  > 
  > proc.time()
     user  system elapsed 
     0.15    0.06    0.21 
Running the tests in 'tests/test-plot.cpr_cp.R' failed.
Complete output:
  > 
  > proc.time()
     user  system elapsed 
     0.20    0.03    0.21 
Running the tests in 'tests/test-summary.cpr_cn.R' failed.
Complete output:
  > library(cpr)
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.29    0.04    0.34 
Running the tests in 'tests/test-summary.cpr_cpr.R' failed.
Complete output:
  > library(cpr)
  > 
  > set.seed(42)
  > x <- seq(0 + 1/5000, 6 - 1/5000, length.out = 100)
  > bmat <- bsplines(x, iknots = c(1, 1.5, 2.3, 4, 4.5), bknots = c(0, 6))
  > theta <- matrix(c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5), ncol = 1)
  > DF <- data.frame(x = x, truth = as.numeric(bmat %*% theta))
  > DF$y <- as.numeric(bmat %*% theta + rnorm(nrow(bmat), sd = 0.3))
  > 
  > initial_cp <-
  +   cp(y ~ bsplines(x, iknots = c(1, 1.5, 2.3, 3.0, 4, 4.5), bknots = c(0, 6))
  +      , data = DF
  +      , keep_fit = TRUE # default is FALSE
  +   )
  > 
  > cpr0 <- cpr(initial_cp)
  
    |                                                                            
    |                                                                      |   0%
    |                                                                            
    |==========                                                            |  14%
    |                                                                            
    |====================                                                  |  29%
    |                                                                            
    |==============================                                        |  43%
    |                                                                            
    |========================================                              |  57%
    |                                                                            
    |==================================================                    |  71%
    |                                                                            
    |============================================================          |  86%
    |                                                                            
    |======================================================================| 100%
  > s <- summary(cpr0)
  > stopifnot(identical(nrow(s), 7L))
  > stopifnot(identical(names(s), c("dfs", "n_iknots", "iknots", "loglik", "rss", "rse", "wiggle", "fdsc", "Pr(>w_(1))")))
  > 
  > stopifnot(isTRUE(
  +   all.equal(summary(cpr0[[1]])
  +             , s[1, c("dfs", "n_iknots", "iknots", "loglik", "rss", "rse", "wiggle", "fdsc")]
  +             , check.attributes = FALSE
  + )))
  > 
  > stopifnot(identical(
  +     attr(s, "elbow")
  +     ,
  +     structure(c(3, 3, 3, 3, 3, 3), dim = 2:3, dimnames = list(c("quadratic", "linear"), c("loglik", "rss", "rse")))
  +     )
  + )
  > 
  > # verify that print returns the object
  > printed <- print(s)
    dfs n_iknots       iknots    loglik       rss       rse    wiggle fdsc
  1   4        0              -74.49596 25.977023 0.5201865  46.88727    2
  2   5        1          4.5 -60.10762 19.481116 0.4528403  46.14794    2
  3   6        2       1, 4.5 -22.06652  9.103178 0.3111950 113.07351    4
  4   7        3    1, 3, 4.5 -19.96710  8.728863 0.3063637  99.28601    4
  5   8        4 1, 1.5, .... -19.95694  8.727089 0.3079930  95.78576    4
  6   9        5 1, 1.5, .... -19.79792  8.699377 0.3091885  92.01369    4
  7  10        6 1, 1.5, .... -19.67424  8.677886 0.3105172  87.35118    4
    Pr(>w_(1))
  1         NA
  2 1.8203e-08
  3 < 2.22e-16
  4 9.6222e-05
  5   0.630871
  6   0.072342
  7   0.066713
  
  -------
  Elbows (index of selected model):
            loglik rss rse
  quadratic      3   3   3
  linear         3   3   3
  > stopifnot(identical(printed, s))
  > 
  > # only the first and last row:
  > printed <- print(s, n = 1)
      dfs n_iknots       iknots            loglik              rss
  1     4        0              -74.4959645891524 25.9770228656933
  ---                                                             
  7    10        6 1, 1.5, .... -19.6742417303079 8.67788603091469
                    rse           wiggle fdsc Pr(>w_(1))
  1   0.520186493658097 46.8872701952964    2         NA
  ---                                                   
  7   0.310517239294459 87.3511819271172    4   0.066713
  
  -------
  Elbows (index of selected model):
            loglik rss rse
  quadratic      3   3   3
  linear         3   3   3
  > stopifnot(identical(printed, s))
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     1.70    0.25    1.93 
Running the tests in 'tests/test-wiggle.R' failed.
Complete output:
  > # verify wiggle returns the integeral of the squared second derivative
  > 
  > library(cpr)
  > 
  > ################################################################################
  > f <- function(x) {
  +   #(x + 2) * (x - 1) * (x - 3)
  +   x^3 - 2 * x^2 - 5 * x + 6
  + }
  > 
  > fprime <- function(x) { # first derivatives of f(x)
  +   3 * x^2 - 4 * x - 5
  + }
  > 
  > fdoubleprime <- function(x) { # second derivatives of f(x)
  +   6 * x - 4
  + }
  > 
  > fdoubleprime_squared <- function(x) {
  +   fdoubleprime(x)^2
  + }
  > 
  > bknots = c(-3, 5)
  > x     <- (runif(n = 100, min = -3, max = 5))
  > bmat  <- bsplines(x, bknots = bknots)
  > theta <- matrix(coef(lm(f(x) ~ bmat + 0)), ncol = 1)
  > 
  > bmatD1 <- bsplineD(x, bknots = bknots, derivative = 1L)
  > bmatD2 <- bsplineD(x, bknots = bknots, derivative = 2L)
  > 
  > stopifnot(isTRUE(all.equal(
  +   integrate(fdoubleprime_squared, lower = -3, upper = 5)$value
  +   ,
  +   wiggle(cp(bmat, theta), lower = -3, upper = 5)$value
  + )))
  > 
  > ################################################################################
  > # verify the result is the same if the bounds are not specified
  > stopifnot(isTRUE(all.equal(
  +   integrate(fdoubleprime_squared, lower = -3, upper = 5)$value
  +   ,
  +   wiggle(cp(bmat, theta))$value
  + )))
  > 
  > ################################################################################
  > # verify the wiggle result is the same if the bounds are set too wide
  > stopifnot(isTRUE(all.equal(
  +   integrate(fdoubleprime_squared, lower = -3, upper = 5)$value
  +   ,
  +   wiggle(cp(bmat, theta, lower = -10, upper = 20))$value
  + )))
  > 
  > ################################################################################
  > #                              Count Sign Changes                              #
  > 
  > # sign changes in fprime
  > # the data need to be sorted for fprime call
  > stopifnot(identical(
  +   sum(abs(diff(sign(fprime(sort(x)))))) / 2L
  +   ,
  +   sign_changes( object = cp(bmat, theta), lower = -3, upper = 5, derivative = 1)
  +   )
  + )
  > 
  > # sign_changes in fdoubleprime
  > stopifnot(identical(
  +   sum(abs(diff(sign(fdoubleprime(sort(x)))))) / 2L
  +   ,
  +   sign_changes( object = cp(bmat, theta), lower = -3, upper = 5, derivative = 2)
  +   )
  + )
  > 
  > ################################################################################
  > # verify the sign_changes result is the same if the bounds are set too wide
  > 
  > warn <- tryCatch(sign_changes(cp(bmat, theta), lower = -10, upper = 20, derivative = 1)
  +                  , warning = function(w) w)
  > stopifnot(inherits(warn, "warning"))
  > stopifnot(warn$message == "setting lower = min(object$bknots)")
  > 
  > 
  > stopifnot(isTRUE(all.equal(
  +   sign_changes( object = cp(bmat, theta), lower = -3, upper = 5, derivative = 1)
  +   ,
  +   suppressWarnings(sign_changes(cp(bmat, theta), lower = -10, upper = 20, derivative = 1))
  + )))
  > 
  > 
  > ################################################################################
  > #                                 End of File                                  #
  > ################################################################################
  > 
  > proc.time()
     user  system elapsed 
     0.26    0.07    0.32