Nothing
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)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.