set.seed(42)
# first load the 46 statin data
data(statin46)
# zero inflation probabilities
omega_tru <- c(rep(0.15, ncol(statin46) - 1), 0)
# signal matrix
signal_mat <- matrix(1, nrow(statin46), ncol(statin46))
# "positive" signal at the (1, 1) entry
# the first column
signal_mat[1, 1] <- 10
# Now simulate data with the same row/column marginals
# as in statin46, with embedded signals as specified in
# the above signal_mat
# no zero inflation at (1, 1)
# (where signal is elicited) and the last row ("Other PT")
# and at the last column ("Other drugs") of the simulated matrix
sim_statin <- r_contin_table_zip(
n = 1,
row_marginals = rowSums(statin46),
col_marginals = colSums(statin46),
signal_mat = signal_mat,
omega_vec = omega_tru,
no_zi_idx = list(
c(1, 1),
c(nrow(statin46), 0), # the entire last row
c(0, ncol(statin46)) # the entire last column
)
)[[1]]
# now analyze the above simulated data
# using a pseudo LRT with a ZIP model
test1 <- pvlrt(
contin_table = sim_statin,
nsim = 500
# set to 500 for demonstration purposes only,
# we recommend the default 10000 or bigger
)
extract_zi_probability(test1)
extract_significant_pairs(test1)
# LRT with a poisson model
test2 <- lrt_poisson(
contin_table = sim_statin,
nsim = 500
# set to 500 for demonstration purposes only,
# we recommend the default 10000 or bigger
)
extract_zi_probability(test2)
extract_significant_pairs(test2)
# LRT with true omega supplied
test3 <- pvlrt(
contin_table = sim_statin,
zi_prob = omega_tru,
nsim = 500
# set to 500 for demonstration purposes only,
# we recommend the default 10000 or bigger
)
extract_zi_probability(test3)
extract_significant_pairs(test3)
Run the code above in your browser using DataLab