tests/testthat/testExampleCode.R

library(lme4)

context("Run Examples")

set.seed(1)

test_that("gen_design example code runs without errors", {
  skip_on_cran()
  #'#Generating a basic 2 factor design:
  basicdesign <- expand.grid(x1 = c(-1, 1), x2 = c(-1, 1))
  set.seed(1)
  #'
  #'#This factorial design is used as an input in the optimal design generation for a
  #'#D-optimal design with 11 runs.
  expect_silent({design <- gen_design(candidateset = basicdesign, model = ~x1 + x2, trials = 11)})
  #'
  #'#We can also use the dot operator to automatically use all of the terms in the model:
  expect_silent({design <- gen_design(candidateset = basicdesign, model = ~., trials = 11)})
  #'
  #'#Here we add categorical factors, specified by using "as.factor" in expand.grid:
  expect_silent({
    categoricaldesign <- expand.grid(a = c(-1, 1), b = as.factor(c("A", "B")), c = as.factor(c("High", "Med", "Low")))})
  #'
  #'#This factorial design is used as an input in the optimal design generation.
  expect_silent({
    design2 = gen_design(candidateset = categoricaldesign, model = ~a + b + c, trials = 19)})
  #'
  #'#We can also increase the number of times the algorithm repeats the search to increase the probability
  #'#that the globally optimal design was found.
  expect_silent({
    design2 = gen_design(candidateset = categoricaldesign, model = ~a + b + c, trials = 19, repeats = 100)})
  #'
  #'#You can also use a higher order model when generating the design:
  expect_silent({
    design2 = gen_design(candidateset = categoricaldesign, model = ~a + b + c + a * b * c, trials = 12)
    })
  #'
  #'#To evaluate a response surface design, include center points in the candidate set and do not include
  #'#quadratic effects with categorical factors.
  #'
  expect_silent({designquad = expand.grid(a = c(1, 0, -1), b = c(-1, 0, 1), c = c("A", "B", "C"))})
  #'
  expect_silent({design = gen_design(designquad, ~a + b + I(a ^ 2) + I(b ^ 2) + a * b * c, 20,timer =FALSE)})
  #'
  #'#The optimality criterion can also be changed:
  expect_silent({design = gen_design(designquad, ~a + b + I(a ^ 2) + I(b ^ 2) + a * b * c, 20, optimality = "I",timer =FALSE)})
  expect_silent({design = gen_design(designquad, ~a + b + I(a ^ 2) + I(b ^ 2) + a * b * c, 20, optimality = "A",timer =FALSE)})
  expect_silent({design = gen_design(designquad, ~a + b + I(a ^ 2) + I(b ^ 2) + a * b * c, 20, optimality = "D",timer =FALSE)})
  expect_silent({design = gen_design(designquad, ~a + b + I(a ^ 2) + I(b ^ 2) + a * b * c, 20, optimality = "E",timer =FALSE)})
  expect_silent({design = gen_design(designquad, ~a + b + I(a ^ 2) + I(b ^ 2) + a * b * c, 20, optimality = "T",timer =FALSE)})
  expect_silent({design = gen_design(designquad, ~a + b + I(a ^ 2) + I(b ^ 2), 20, optimality = "G",timer =F)})
  expect_silent({design = gen_design(designquad, ~a + b + I(a ^ 2) + I(b ^ 2), 20, optimality = "Alias",timer =F)})
  #'
  #'#A split-plot design can be generated by first generating an optimal blocking design using the
  #'#hard-to-change factors and then using that as the input for the split-plots design.
  #'#This generates an optimal subplot design that accounts for the existing split-plot settings.
  #'#See the accompannying paper "___________" for details of the implementation.
  #'
  candlist1 = expand.grid(Altitude = c(-1, 1), Range = as.factor(c("Close", "Medium", "Far")), Power = c(1, -1))
  expect_silent({hardtochangedesign = gen_design(candidateset = candlist1, model = ~Altitude, trials = 11)})
  #'
  #'#Now we can use the D-optimal blocked design as an input to our full design.
  #'
  #'#Here, we specify the easy to change factors for the factorial design, and input the hard-to-change design
  #'#along with a vector listing the number of repetitions within each block for the blocked design. There should be
  #'#a size entry for every block and the number of runs specified in the trials argument needs to equal the
  #'#sum of all of the block sizes or else the program will throw an error.
  #'
  #'#Putting this all together:
  expect_silent({
    designsplitplot = gen_design(candlist1, ~Altitude + Range + Power, trials = 33, splitplotdesign = hardtochangedesign,
                                 blocksizes = 3, timer = FALSE)})

  expect_silent({design = gen_design(candlist1, ~Altitude + Range + Power, trials = 33,
                                     splitplotdesign = hardtochangedesign, blocksizes = 3, optimality = "I", timer = FALSE)})
  expect_silent({design = gen_design(candlist1, ~Altitude + Range + Power, trials = 33,
                                     splitplotdesign = hardtochangedesign, blocksizes = 3, optimality = "A", timer = FALSE)})
  expect_silent({design = gen_design(candlist1, ~Altitude + Range + Power, trials = 33,
                                     splitplotdesign = hardtochangedesign, blocksizes = 3, optimality = "D", timer = FALSE)})
  expect_silent({design = gen_design(candlist1, ~Altitude + Range + Power, trials = 33,
                                     splitplotdesign = hardtochangedesign, blocksizes = 3, optimality = "E", timer = FALSE)})
  expect_silent({design = gen_design(candlist1, ~Altitude + Range + Power, trials = 33,
                                     splitplotdesign = hardtochangedesign, blocksizes = 3, optimality = "T", timer = FALSE)})
  expect_silent({design = gen_design(candlist1, ~Altitude + Range + Power, trials = 33,
                                     splitplotdesign = hardtochangedesign, blocksizes = 3, optimality = "Alias", timer = FALSE)})
  expect_silent({design = gen_design(candlist1, ~Altitude + Range + Power, trials = 33,
                                     splitplotdesign = hardtochangedesign, blocksizes = 3, optimality = "G", timer = FALSE)})
  #testing custom
  # customOpt = function(currentDesign) {
  #   return(det(t(currentDesign) %*% currentDesign))
  # }
  #
  # customBlockedOpt = function(currentDesign, vInv) {
  #   return(det(t(currentDesign) %*% vInv %*% currentDesign))
  # }
  #
  # candlist1 = expand.grid(Altitude = c(-1, 0, 1), Range = as.factor(c("Close", "Medium", "Far")), Power = c(-1, 0, 1))
  # hardtochangedesign = gen_design(candidateset = candlist1, model = ~Altitude, trials = 11, optimality = "custom")
  #
  # splitplotblocksize = rep(3, 11)
  # design = gen_design(candlist1, ~Altitude + Range + Power, trials = 33, splitplotdesign = hardtochangedesign, blocksizes = splitplotblocksize, optimality = "custom")
  #'
  #'#The split-plot structure is encoded into the row names, with a period demarcating the blocking level. This process
  #'#can be repeated for arbitrary levels of blocking (i.e. a split-plot design can be entered in as the hard-to-change
  #'#to produce a split-split-plot design, which can be passed as another hard-to-change design to produce a
  #'#split-split-split plot design, etc).
  #'
  candlist3 = expand.grid(Location = as.character(c("East", "West")),
                          Climate = as.factor(c("Dry", "Wet", "Arid")),
                          Vineyard = as.factor(c("A", "B", "C", "D")),
                          Age = c(1, -1))
  #'
  expect_silent(gen_design(candlist3, ~Location, trials = 6) -> temp)
  expect_silent(gen_design(candlist3, ~Location + Climate, trials = 12, splitplotdesign = temp, blocksizes = rep(2, 6), timer = FALSE) -> temp)
  expect_silent(gen_design(candlist3, ~Location + Climate + Vineyard, 48, splitplotdesign = temp, blocksizes = rep(4, 12), timer = FALSE) -> temp)
  expect_silent(gen_design(candlist3, ~Location + Climate + Vineyard + Age, 192, splitplotdesign = temp, blocksizes = rep(4, 48), timer = FALSE) -> splitsplitsplitplotdesign)
  #'
  #'test user-supplied block detection
  externaldesign  = data.frame(Block1 = c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8), var1 = c(-1, 1))
  externaldesign2 = data.frame(Block1 = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4), Block2 = c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8), var1 = c(-1, 1))
  externaldesign3 = data.frame(Block2 = c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8), var1 = c(-1, 1), Block1 = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4))
  expect_warning(eval_design_mc(externaldesign,  ~., 0.2, nsim = 10), "ignoring blocking structure and removing blocking columns.")
  expect_warning(eval_design_mc(externaldesign2, ~., 0.2, nsim = 10), "ignoring blocking structure and removing blocking columns.")
  expect_warning(eval_design_mc(externaldesign3, ~., 0.2, nsim = 10), "ignoring blocking structure and removing blocking columns.")
  expect_warning(eval_design(externaldesign,  ~., 0.2), "ignoring blocking structure and removing blocking columns.")
  expect_warning(eval_design(externaldesign2, ~., 0.2), "ignoring blocking structure and removing blocking columns.")
  expect_warning(eval_design(externaldesign3, ~., 0.2), "ignoring blocking structure and removing blocking columns.")
  expect_silent(eval_design(externaldesign, ~., 0.2, blocking = TRUE))
  expect_silent(eval_design(externaldesign2, ~., 0.2, blocking = TRUE))
  expect_silent(eval_design(externaldesign3, ~., 0.2, blocking = TRUE))

  #test interaction between whole and subplots

  expect_silent(gen_design(candlist3, ~Location, trials = 6) -> temp)
  expect_silent(gen_design(candlist3, ~Location + Climate + Location:Climate, trials = 12, splitplotdesign = temp, blocksizes = rep(2, 6), timer = FALSE) -> temp2)
  expect_silent(gen_design(candlist3, ~Location + Climate + Location * Climate, trials = 12, splitplotdesign = temp, blocksizes = rep(2, 6), timer = FALSE) -> temp2)

  #test parallel features

  options(cores = 2)

  expect_silent(gen_design(candlist3, ~Location, trials = 6, parallel = TRUE) -> temp)
  expect_silent(gen_design(candlist3, ~Location + Climate, trials = 12, splitplotdesign = temp, blocksizes = rep(2, 6), parallel = TRUE, timer = FALSE) -> temp)
  expect_silent(gen_design(candlist3, ~Location, trials = 6, parallel = TRUE) -> temp)
  expect_silent(gen_design(candlist3, ~Location + Climate, trials = 12, splitplotdesign = temp, blocksizes = rep(2, 6), parallel = TRUE, timer = FALSE) -> temp2)
  expect_silent(eval_design_mc(temp, ~., 0.2, nsim = 10, parallel = TRUE))
  expect_silent(eval_design_mc(temp2, ~., 0.2, nsim = 10, parallel = TRUE))
  expect_silent(eval_design_mc(temp, ~., 0.2, nsim = 10, glmfamily = "poisson", effectsize = c(1, 10), parallel = TRUE))
  expect_silent(eval_design_mc(temp2, ~., 0.2, nsim = 10, glmfamily = "poisson", effectsize = c(1, 10), parallel = TRUE))
  expect_warning(eval_design_mc(temp, ~., 0.2, nsim = 10, glmfamily = "poisson"), " This can lead to unrealistic effect sizes")



  #A design's diagnostics can be accessed via the following attributes:

  expect_silent(attr(design, "D")) #D-Efficiency
  expect_silent(attr(design, "A")) #A-Efficiency
  expect_silent(attr(design, "I")) #The average prediction variance across the design space
  #'
  #'#The correlation matrix can be accessed via the "correlation.matrix" attribute:
  expect_silent({correlation.matrix = attr(design2, "correlation.matrix")})


})


test_that("eval_design example code runs without errors", {
  #'#this can also be generated with expand.grid as:
  #'
  factorialdes <- expand.grid(A = c(1, -1), B = c(1, -1), C = c(1, -1))
  expect_silent({optdesign = gen_design(candidateset = factorialdes, model = ~A + B + C, trials = 17, optimality = "D", repeats = 100)})
  #'
  #'#Now evaluating that design (with default anticipated coefficients and a effectsize of 2):
  expect_silent(eval_design(design = optdesign, model = ~A + B + C, alpha = 0.2))
  #'
  #'#Evaluating a subset of the design (changing the power due to a different number of
  #'#degrees of freedom)
  expect_silent(eval_design(design = optdesign, model = ~A + C, alpha = 0.2))
  #'
  #'#Halving the signal-to-noise ratio by setting a different effectsize (default is 2):
  expect_silent(eval_design(design = optdesign, model = ~A + B + C, alpha = 0.2, effectsize = 1))

  #'#Trying with ~.*. operator
  expect_silent(eval_design(design = optdesign, model = ~. * ., alpha = 0.2, effectsize = 1))

  #'
  #'#With 3+ level categorical factors, the choice of anticipated coefficients directly changes the
  #'#final power calculation. For the most conservative power calculation, that involves
  #'#setting all anticipated coefficients in a factor to zero except for one. We can specify this
  #'#option with the "conservative" argument.
  #'
  expect_silent({factorialcoffee = expand.grid(caffeine = c(1, -1),
                                               cost = c(1, 2),
                                type = as.factor(c("Kona", "Colombian", "Ethiopian", "Sumatra")),
                                size = as.factor(c("Short", "Grande", "Venti")))})
  #'
  expect_silent({designcoffee = gen_design(factorialcoffee, ~cost + size + type, trials = 29, optimality = "D", repeats = 100)})
  #'
  #'#Evaluate the design, with default anticipated coefficients (conservative is FALSE by default).
  expect_silent(eval_design(designcoffee, model = ~cost + size + type, alpha = 0.05))
  #'
  #'#Evaluate the design, with conservative anticipated coefficients:
  expect_silent(eval_design(designcoffee, model = ~cost + size + type, alpha = 0.05, conservative = TRUE))
  #'
  #'#which is the same as the following, but now explicitly entering the coefficients:
  expect_silent(eval_design(designcoffee, model = ~cost + size + type, alpha = 0.05, anticoef = c(1, 1, 1, 0, 0, 1, 0)))
  #'
  #'#If the first level in a factor is not the one that you want to set to one
  #'#in the conservative calculation, enter the anticipated coefficients in manually.
  expect_silent(eval_design(designcoffee, model = ~cost + size + type, alpha = 0.05, anticoef = c(1, 1, 0, 0, 1, 0, 1)))
  #'
  #'#You can also evaluate the design with higher order effects:
  expect_silent(eval_design(designcoffee, model = ~cost + size + type + cost * type, alpha = 0.05))
  #'
  #'#Blocked designs can also be evaluated by specifying the blocking model.
  #'
  #'#Generating blocked design
  expect_silent({coffeeblockdesign = gen_design(factorialcoffee, ~caffeine, trials = 12)})
  expect_silent({coffeefinaldesign = gen_design(factorialcoffee, model = ~caffeine + cost + size + type, trials = 36,
                                 splitplotdesign = coffeeblockdesign)})
  #'
  #'#Evaluating design
  expect_silent(eval_design(coffeefinaldesign,  model = ~cost + size + type + caffeine, alpha = 0.2, blocking = TRUE))
  #'
  #'#We can also evaluate the design with a custom ratio between the whole plot error to
  #'#the run-to-run error.
  expect_silent(eval_design(coffeefinaldesign,  model = ~cost + size + type + caffeine, alpha = 0.2, blocking = TRUE,
              varianceratio = 2))

  expect_silent(eval_design(coffeefinaldesign,  model = ~cost + size + type + caffeine, alpha = 0.2, blocking = TRUE,
                             varianceratio = c(2, 1)))
})


test_that("eval_design_mc example code runs without errors", {
  set.seed(123)
  #'@examples #We first generate a full factorial design using expand.grid:
  factorialcoffee = expand.grid(cost = c(-1, 1),
                                type = as.factor(c("Kona", "Colombian", "Ethiopian", "Sumatra")),
                                size = as.factor(c("Short", "Grande", "Venti")))
  expect_silent({
    designcoffee = gen_design(factorialcoffee, model = ~cost + type + size, trials = 21, optimality = "D")
  })
  expect_silent(eval_design(design = designcoffee, model = ~cost + type + size, 0.05))
  expect_silent(eval_design_mc(design = designcoffee, model = ~cost + type + size, alpha = 0.05, nsim = 100,
                               glmfamily = "gaussian"))
  expect_silent(eval_design_mc(design = designcoffee, model = ~cost + type + size, alpha = 0.05,
                               nsim = 100, glmfamily = "gaussian", effectsize = 1))
  expect_silent(eval_design_mc(design = designcoffee, model = ~cost + type, alpha = 0.05,
                 nsim = 100, glmfamily = "gaussian"))
  expect_silent(eval_design_mc(design = designcoffee, model = ~cost + type + size, 0.05,
                 nsim = 100, glmfamily = "gaussian"))
  expect_silent(eval_design_mc(design = designcoffee, model = ~cost + type + size + cost * type, 0.05,
                 nsim = 100, glmfamily = "gaussian"))
  #'\dontrun{eval_design_mc(design = designcoffee, model = ~cost + type + size, 0.05,
  #'               nsim = 10000, glmfamily = "gaussian", parallel = TRUE)}
  #'
  #'
  #'#Trying with ~.*. operator
  expect_silent(eval_design_mc(design = designcoffee, model = ~. * ., 0.05,
                               nsim = 100, glmfamily = "gaussian"))

  factorialcoffee = expand.grid(Temp = c(1, -1),
                                Store = as.factor(c("A", "B")),
                                cost = c(-1, 1),
                                type = as.factor(c("Kona", "Colombian", "Ethiopian", "Sumatra")),
                                size = as.factor(c("Short", "Grande", "Venti")))

  vhtcdesign = gen_design(candidateset = factorialcoffee, model = ~Store, trials = 8)
  htcdesign = gen_design(candidateset = factorialcoffee, model = ~Store + Temp, trials = 24, splitplotdesign = vhtcdesign, blocksizes = 3)
  expect_silent({
    splitplotdesign = gen_design(candidateset = factorialcoffee, model = ~Store + Temp + cost + type + size, trials = 96,
                               splitplotdesign = htcdesign, blocksizes = 4, varianceratio = 3)
  })
  expect_silent(eval_design_mc(design = splitplotdesign, model = ~Store + Temp + cost + type + size, alpha = 0.05, blocking = TRUE,
                                nsim = 1, glmfamily = "gaussian", varianceratios = c(5, 4)))
  expect_silent(eval_design_mc(design = splitplotdesign, model = ~Store + Temp + cost + type + size, alpha = 0.05, blocking = TRUE,
                               nsim = 1, glmfamily = "gaussian", varianceratios = c(5, 4)))
  expect_error(eval_design_mc(design = splitplotdesign, model = ~Store + Temp + cost + type + size, alpha = 0.05, blocking = TRUE,
                                nsim = 1, glmfamily = "gaussian", varianceratios = c(5, 4, 2, 2)),"Wrong number of variance ratios specified.")
  expect_silent(eval_design_mc(design = splitplotdesign, model = ~Store + Temp + cost + type + size, alpha = 0.05, blocking = TRUE,
                                nsim = 1, glmfamily = "gaussian", varianceratios = c(5, 4, 2)))
  expect_silent(eval_design_mc(design = splitplotdesign, model = ~Store + Temp + cost + type + size, alpha = 0.05, blocking = TRUE,
                                nsim = 1, glmfamily = "gaussian", varianceratios = c(5, 4)))

  factorialbinom = expand.grid(a = c(-1, 1), b = c(-1, 1))
  expect_silent({
    designbinom = gen_design(factorialbinom, model = ~a + b, trials = 90, optimality = "D", repeats = 100)
  })
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "gaussian"))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "binomial"))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "poisson"))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "exponential"))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "binomial", advancedoptions = list(anovatest = "LR")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "poisson", advancedoptions = list(anovatest = "LR")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "exponential", advancedoptions = list(anovatest = "LR")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "binomial", advancedoptions = list(anovatest = "F")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "poisson", advancedoptions = list(anovatest = "F")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "exponential", advancedoptions = list(anovatest = "F")))

  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "gaussian", advancedoptions = list(anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "binomial", advancedoptions = list(anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "poisson", advancedoptions = list(anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "exponential", advancedoptions = list(anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "binomial", advancedoptions = list(anovatest = "LR", anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "poisson", advancedoptions = list(anovatest = "LR", anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "exponential", advancedoptions = list(anovatest = "LR", anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "binomial", advancedoptions = list(anovatest = "F", anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "poisson", advancedoptions = list(anovatest = "F", anovatype = "II")))
  expect_silent(eval_design_mc(designbinom, ~a + b, alpha = 0.2, nsim = 100, anticoef = c(1.5, 0.7, 0.7),
                               glmfamily = "exponential", advancedoptions = list(anovatest = "F", anovatype = "II")))

  factorialpois = expand.grid(a = as.numeric(c(-1, 0, 1)), b = c(-1, 0, 1))
  designpois = gen_design(factorialpois, ~a + b, trials = 90, optimality = "D", repeats = 100)

})


test_that("eval_design_survival_mc example code runs without errors", {
  basicdesign = expand.grid(a = c(-1, 1))
  expect_silent({
    design = gen_design(candidateset = basicdesign, model = ~a, trials = 100,
                            optimality = "D", repeats = 100)
  })
  rsurvival = function(X, b) {
    Y = rexp(n = nrow(X), rate = exp(-(X %*% b)))
    censored = Y > 1
    Y[censored] = 1
    return(survival::Surv(time = Y, event = !censored, type = "right"))
  }
  expect_warning({
    a = eval_design_survival_mc(design = design, model = ~a, alpha = 0.05, nsim = 100,
                          distribution = "exponential", rfunctionsurv = rsurvival, effectsize = 1)
  },
  "default or length 1 delta used with distribution == 'exponential'. This can lead to unrealistic effect sizes - make sure the generated anticipated coeffcients are appropriate.")

  rlognorm = function(X, b) {
    Y = rlnorm(n = nrow(X), meanlog = X %*% b, sdlog = 0.4)
    censored = Y > 1.2
    Y[censored] = 1.2
    return(survival::Surv(time = Y, event = !censored, type = "right"))
  }
  expect_silent(
    eval_design_survival_mc(design = design, model = ~a, alpha = 0.2, nsim = 100,
                          distribution = "lognormal", rfunctionsurv = rlognorm,
                          anticoef = c(0.184, 0.101), scale = 0.4)
  )
  #testing parallel
  options(cores = 2)
  eval_design_survival_mc(design = design, model = ~a, alpha = 0.05, effectsize = c(1, 2),
                          nsim = 100, distribution = "exponential",
                           censorpoint = 5, censortype = "right", parallel = TRUE)
  options(cores = NULL)
})

test_that("eval_design_custom_mc example code runs without errors", {
  #'@examples #To demonstrate how a user can use their own libraries for Monte Carlo power generation,
  #'#We will recreate eval_design_survival_mc using the eval_design_custom_mc framework.
  #'
  #'#To begin, first let us generate the same design and random generation function shown in the
  #'#eval_design_survival_mc examples:
  #'
  basicdesign = expand.grid(a = c(-1, 1))
  expect_silent(design <- gen_design(candidateset = basicdesign, model = ~a, trials = 100,
                            optimality = "D", repeats = 100))
  #'
  #'#Random number generating function
  #'
  rsurvival = function(X, b) {
    Y = rexp(n = nrow(X), rate = exp(-(X %*% b)))
    censored = Y > 1
    Y[censored] = 1
    return(survival::Surv(time = Y, event = !censored, type = "right"))
  }
  #'
  #'#We now need to tell the package how we want to fit our data,
  #'#given the formula and the model matrix X (and, if needed, the list of contrasts).
  #'#If the contrasts aren't required, "contrastlist" should be set to NULL.
  #'#This should return some type of fit object.
  #'
  fitsurv = function(formula, X, contrastlist = NULL) {
    return(survival::survreg(formula, data = X, dist = "exponential"))
  }
  #'
  #'
  #'#We now need to tell the package how to extract the p-values from the fit object returned
  #'#from the fit function. This is how to extract the p-values from the survreg fit object:
  #'
  pvalsurv = function(fit) {
    return(summary(fit)$table[, 4])
  }
  #'
  #'#And now we evaluate the design, passing the fitting function and p-value extracting function
  #'#in along with the standard inputs for eval_design_mc.
  #'
  expect_silent(eval_design_custom_mc(design = design, model = ~a, alpha = 0.05, nsim = 100,
                        fitfunction = fitsurv, pvalfunction = pvalsurv, rfunction = rsurvival, effectsize = 1))

  #trying with ~. operator
  expect_silent(eval_design_custom_mc(design = design, model = ~., alpha = 0.05, nsim = 100,
                                      fitfunction = fitsurv, pvalfunction = pvalsurv, rfunction = rsurvival, effectsize = 1))

  #testing parallel

  options(cores = c("localhost", "localhost"))
  basicdesign = expand.grid(a = c(-1, 1))
  design = gen_design(candidateset = basicdesign, model = ~a, trials = 100,
                      optimality = "D", repeats = 100)

  rsurvival = function(X, b) {
    Y = rexp(n = nrow(X), rate = exp(-(X %*% b)))
    censored = Y > 1
    Y[censored] = 1
    return(survival::Surv(time = Y, event = !censored, type = "right"))
  }

  fitsurv = function(formula, X, contrastslist = NULL) {
    return(survival::survreg(formula, data = X, dist = "exponential"))
  }

  pvalsurv = function(fit) {
    return(summary(fit)$table[, 4])
  }

  expect_silent({d = eval_design_custom_mc(design = design, model = ~a, alpha = 0.05, nsim = 100,
                          fitfunction = fitsurv, pvalfunction = pvalsurv, rfunction = rsurvival, effectsize = 1, parallel = TRUE)})
  options(cores = NULL)
})

Try the skpr package in your browser

Any scripts or data that you put into this service are public.

skpr documentation built on July 9, 2023, 7:23 p.m.