data("statin46")
# 500 bootstrap iterations (nsim) in each example below
# are for quick demonstration only --
# we recommended setting nsim to 10000 (default) or bigger
# no grouping -- each drug forms its own class,
# default "rrr" (lambda) parametrization, possible zi,
# null distribution evaluated using parametric bootstrap
test1 <- pvlrt(statin46, nsim = 500)
test1
## extract the observed LRT statistic
extract_lrstat_matrix(test1)
## extract the estimated omegas
extract_zi_probability(test1)
# grouped drugs --
# group 1: drug 1, drug 2
# group 2: drug 3, drug 4
# drug 5, 6, in their own groups
## 7 is not tested, so excluded from test_drug_idx (default)
## if needed, include 7 in test_drug_idx
drug_groups <- list(c(1, 2), c(3, 4))
## 5, 6 not present in drug_groups, so each will form their own groups
set.seed(50)
##
test2 <- pvlrt(statin46, drug_class_idx = drug_groups, nsim = 500)
test2
# instead of column positions column names can also be supplied
# in drug_class_idx and/or test_drug_idx
## column name version of drug_groups
drug_groups_colnames <- lapply(drug_groups, function(i) colnames(statin46)[i])
test_drug_colnames <- head(colnames(statin46), -1)
set.seed(50)
test20 <- pvlrt(
statin46,
test_drug_idx = test_drug_colnames,
drug_class_idx = drug_groups_colnames,
nsim = 500
)
test20
all.equal(test2, test20)
# specify no zero inflation at the entirety of the last row and the last column
no_zi_idx <- list(c(nrow(statin46), 0), c(0, ncol(statin46)))
test3 <- pvlrt(statin46, no_zi_idx = no_zi_idx, nsim = 500)
test3
# \donttest{
# use non-parametric bootstrap to evaluate the null distribution
# takes longer, due to computational costs with non-parametric
# contigency table generation
test4 <- pvlrt(statin46, null_boot_type = "non-parametric", nsim = 500)
test4
# }
# test zi probabilities (omegas)
test5 <- pvlrt(statin46, test_omega = TRUE, nsim = 500)
test5
Run the code above in your browser using DataLab