Nothing
# test cases for significance tests for delay
test_that('Structure of test objects.', {
set.seed(123)
x <- rexp_delayed(n = 131L, delay1 = 5, rate1 = .1)
y <- rexp_delayed(n = 111L, delay1 = 5, rate1 = .3)
ted_d <- test_diff(x = x, y = y, distribution = 'expon', param = 'delay1',
R = 19, type = 'bootstrap', chiSqApprox = TRUE)
expect_s3_class(ted_d, class = 'incubate_test')
expect_identical(names(ted_d), c('t_obs', 'testDist', 'R', 'chisq_df_hat', 'param', 'P'))
expect_type(ted_d$P$bootstrap, type = 'double')
expect_lte( ted_d$P$bootstrap, expected = 1L)
expect_gte( ted_d$P$bootstrap, expected = 0L)
})
test_that('GOF-test on single-group exponentials', {
testthat::skip_on_cran()
testthat::skip(message = 'Too long to run every time!')
future::plan(future.callr::callr, workers = 5L)
# GOF-tests on true exponential data with varying sample size, delay and rate
# results in a matrix of dimension #scenarios x #replications
fitting_expos <- future.apply::future_replicate(n = 467L, simplify = FALSE, expr = {
scenarios <- expand.grid(n = c(10, 25, 50), delay1 = c(0, 5, 15), rate1 = c(.01, .2, .4, 1, 1.5, 4))
# fit exponential models with varying n, delay and rate
purrr::pmap(.l = scenarios,
.f = ~ delay_model(x = rexp_delayed(n = ..1, delay1 = ..2, rate1 = ..3),
distribution = 'exponential'))
}) %>% purrr::transpose() # get a list of scenarios, each containing its models of replicated data
# list: for each scenario, the vector of Moran's GOF-test p-value
GOF_pvals <- list(
moran = purrr::map(.x = fitting_expos, .f = ~ purrr::map_dbl(., function(fit) test_GOF(fit, method = 'moran')$p.value)),
pearson = purrr::map(.x = fitting_expos, .f = ~ purrr::map_dbl(., function(fit) test_GOF(fit, method = 'pearson')$p.value))
#ad = purrr::map(.x = fitting_expos, .f = ~ purrr::map_dbl(., function(fit) test_GOF(fit, method = 'ad')$p.value))
)
# # to visualize the GOF P-values:
# as_tibble(GOF_pvals[['moran']], .name_repair = 'unique') %>%
# tidyr::pivot_longer(cols = everything()) %>%
# ggplot(aes(x=value, col = name)) +
# geom_freqpoly(binwidth = .15) + xlim(0, 1)
# expect uniform P-values for GOF under valid H0
# go over the three types of GOF-tests, within test type, check each scenario
# use purrr::flatten(GOF_pvals) as data argument to walk to test *all* in one go
purrr::walk(GOF_pvals[['pearson']], .f = ~expect_gt(object = mean(.), expected = 0.25))
purrr::walk(GOF_pvals[['pearson']], .f = ~expect_lt(object = mean(.), expected = 0.75))
# purrr::walk(GOF_pvals[['ad']], .f = ~expect_gt(object = mean(.), expected = 0.25))
# purrr::walk(GOF_pvals[['ad']], .f = ~expect_lt(object = mean(.), expected = 0.75))
purrr::walk(GOF_pvals[['moran']], .f = ~expect_gt(object = mean(.), expected = 0.35))
purrr::walk(GOF_pvals[['moran']], .f = ~expect_lt(object = mean(.), expected = 0.65))
# Moran's test shows indeed P-values on average close to 0.5
purrr::walk(GOF_pvals[['moran']], .f = ~expect_equal(object = mean(.), expected = 0.5, tolerance = .22, info = 'moran'))
# Moran's test does not give P-values very close to 0, less than expected under uniformity
GOF_moran_pvals_KSpval <- purrr::map_dbl(GOF_pvals[['moran']], .f = ~suppressWarnings(ks.test(.[which(. > .085)],
y = 'punif', min = 0.085)$p.value))
# not too many small P-values
expect_lte(length(which(GOF_moran_pvals_KSpval < .05)) / length(GOF_moran_pvals_KSpval),
expected = .3)
# some high P-values
expect_gte(length(which(GOF_moran_pvals_KSpval > .6)) / length(GOF_moran_pvals_KSpval),
expected = .15)
future::plan(future::sequential)
})
test_that("Test difference in delay for two exponential fits", {
testthat::skip_on_cran()
testthat::skip(message = 'Too long to run every time!')
future::plan(future.callr::callr, workers = 3L)
set.seed(12345)
x <- rexp_delayed(13L, delay1 = 11, rate1 = .05)
y <- rexp_delayed(17L, delay1 = 11, rate1 = .08)
# increasing effect
te_diff_delays <- purrr::map(purrr::set_names(c(0, 9, 19)),
~ test_diff(x = x + .x, y = y, param = "delay1", type = 'bootstrap', R = 399))
te_diff_delays_P_bs <- purrr::map_dbl(te_diff_delays, ~ purrr::chuck(., "P", "bootstrap"))
# null model (no effect) has a high P-value
expect_gt(te_diff_delays_P_bs[['0']], expected = .1)
# the bigger the effect (=difference in delay) the smaller the P-value
expect_true(all(diff(te_diff_delays_P_bs) < 0L))
# negative correlation betw effect size and P-values
expect_lt(cor(x = as.integer(names(te_diff_delays)),
y = te_diff_delays_P_bs),
expected = -.67)
#test effect of sample size: increase n and power goes up.
set.seed(123456)
# data with difference in delay by 2.5 time units
#+but different sample sizes
xs <- purrr::map(purrr::set_names(c(9, 20, 32, 37)),
~ rexp_delayed(., delay1 = 6.5, rate1 = .07))
ys <- purrr::map(purrr::set_names(c(10, 19, 30, 38)),
~ rexp_delayed(., delay1 = 9, rate1 = .07))
te_diff_delays_n <- purrr::map2(.x = xs, .y = ys,
.f = ~ suppressWarnings(test_diff(x = .x, y = .y, param = "delay1", type = 'bootstrap', R = 399)))
expect_lt(cor(x = as.integer(names(te_diff_delays_n)),
y = purrr::map_dbl(te_diff_delays_n, ~ purrr::chuck(., "P", "bootstrap"))),
expected = -.67)
future::plan(future::sequential)
})
test_that("Test difference in delay when H0 is true (no difference in delay)", {
testthat::skip_on_cran()
testthat::skip(message = 'Too long to run every time!')
future::plan(future.callr::callr, workers = 3L)
set.seed(20210506)
testres_P_H0 <- future.apply::future_vapply(X = seq_len(21L), FUN.VALUE = double(1L),
FUN = function(dummy) {
x <- rexp_delayed(13, delay1 = 4, rate1 = .07)
y <- rexp_delayed(11, delay1 = 4, rate1 = .1)
# return P-value of bootstrap test
Pval <- NA_real_
try(Pval <- purrr::chuck(test_diff(x = x, y = y, param = "delay1",
distribution = "expon", R = 301),
"P", "bootstrap"),
silent = TRUE)
Pval
}, future.seed = TRUE)
testres_P_H0 <- testres_P_H0[is.finite(testres_P_H0)]
# KS-test does not reject H0: uniform distribution
expect_gt(suppressWarnings(stats::ks.test(x = testres_P_H0, y = "punif")$p.value), expected = .2)
future::plan(future::sequential)
})
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.