inst/doc/tutorial-lt-mapit.R

In mvMAPIT: Multivariate Genome Wide Marginal Epistasis Test

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(mvMAPIT)

## ----parameters, eval = F-----------------------------------------------------
# population_size <- 1e4
# n_snps <- 2e3
# pve <- 0.6
# rho <- 0.5
# disease_prevalence <- 0.3
# sample_size <- 500

## ----random_genotypes, eval = F-----------------------------------------------
# maf <- 0.05 + 0.45 * runif(n_snps)
# random_genotypes   <- (runif(population_size * n_snps) < maf) + (runif(population_size *
#                                                                          n_snps) < maf)
# random_genotypes   <- matrix(as.double(random_genotypes), population_size, n_snps, byrow = TRUE)

## ----simulate_traits, eval = F------------------------------------------------
# n_causal <- 100
# n_epistatic <- 10
# simulated_data <- simulate_traits(
#     random_genotypes,
#     n_causal = n_causal,
#     n_trait_specific = n_epistatic,
#     n_pleiotropic = 0,
#     d = 1,
#     H2 = pve,
#     rho = rho
# )

## ----case_control, eval = F---------------------------------------------------
# liabilities <- simulated_data$trait
# threshold <- qnorm(1 - disease_prevalence, mean = 0, sd = 1)
# case_control_trait <- rep(0, population_size)
# case_control_trait[liabilities > threshold] <- 1
# 
# # Subsample a particular number of cases and controls
# cases <- sample(which(liabilities > threshold), sample_size / 2, replace = FALSE)
# controls <- sample(which(liabilities <= threshold), sample_size / 2, replace = FALSE)
# y <- as.integer(case_control_trait[c(cases, controls)])
# X <- simulated_data$genotype[c(cases, controls), ]

## ----mapit, eval = F----------------------------------------------------------
# y_liabilities <- binary_to_liability(y, disease_prevalence)
# 
# lt_mapit <- mvmapit(
#         t(X),
#         t(y_liabilities),
#         test = "hybrid"
# )