# Load and process the data to run SOHPIE
data("simulated_data")
n_genes = ncol(simulated_data[, grep("^(M_)", colnames(simulated_data))])
MedZIsc.res <- Magics(data.name = simulated_data, n_genes = n_genes, covariate.names = c("Z1", "Z2", "Z3"))

# Create a function to check the sign of p-values
# and adjusted p-values:
sign_check <- function(pval) {
        if (sum(pval < 0) > 0) {
                return("ERROR: Negative p-values")
        } else {
                return("No negative p-values")
        }
}


# A function to check if results have any negative p-values or q-values:
test_that("Check any negative p-values after Magics", {
        expect_equal(sign_check(MedZIsc.res$adj_pvals_M), "No negative p-values")
        expect_equal(sign_check(MedZIsc.res$adj_pvals_F), "No negative p-values")
})

# Use testthat to check if results return any missing or NA:
# Note: coefficient estimates and p-values
#       (and q-values) should not have any missing or NA.
test_that("Check the missingness after Magics", {
        expect_false(any(is.na(MedZIsc.res$coef_outcome)))
        expect_false(any(is.na(MedZIsc.res$se_outcome)))
})