# Read in the maew library
library(maew)

## --------- Simulate a cohort with genotype, covariates and phenotype ---------- ## 
n.subj = 100;   # The cohort contains 100 subjects
n.snps = 10;    # Each subject corresponds to 10 SNPs
n.covs = 2;     # Each subject corresponds to 2 covariates
maf = 0.3;      # Minor allele frequency for each SNP is 0.3

## Genotype is a numeric matrix of the observed genotype, each row is for one subject, each column is for one feature
gtype.data <- rbinom(n.subj*n.snps,size=2,prob=maf)
Gtype <- matrix(gtype.data,
                nrow=n.subj,ncol=n.snps);
## Covariate is a numeric matrix of the covariates used in the model, each row is for one subject, each column is for one covariate
covars.data <- rnorm(n.subj*n.covs)
Covars <- matrix(covars.data,
                 nrow=n.subj,ncol=n.covs);
## Phenotype is a numeric vector of the phenoytpe values.
Ptype.cont <- rnorm(n.subj);                        # continuous phenotype
Ptype.bina <- rbinom(n.subj,size=1,prob=0.5)        # binary phenotype
Ptype.surv <- Surv(rexp(n.subj),                    # survival phenotype
                   rbinom(n.subj,size=1,prob=0.5))
  




## --------- For each genomic feature, compute AICs and MAEWs for each candidate model.

## Continuous phenotype
bic.cont <- compute.BIC(Gtype,Ptype.cont,Covars,model.string="glm");
maew.cont <- compute.maew(bic.cont)
# aic.cont
# maew.cont

## Binary phenotype
bic.bina <- compute.BIC(Gtype,Ptype.bina,Covars,model.string="glm",family=binomial);
maew.bina <- compute.maew(bic.bina)
# aic.bina
# maew.bina

## Survival phenotype
bic.surv <- compute.BIC(Gtype,Ptype.surv,Covars,model.string="coxph");
maew.surv <- compute.maew(bic.surv)
# aic.surv
# maew.surv


## ----------- function to compute EBP based on p-value histogram

## 10000 p-values under null hypothesis
p <- runif(10000,0,1)
res <- pval.ebp(p)
EBP <- res[,"EBP"]

## 9970 p-values under null hypothesis and 30 p-values under alternative hypothesis
p <- c(runif(30,0,1e-4),runif(9970,0,1))
res <- pval.ebp(p)
EBP <- res[,"EBP"]

## draw plot
layout(matrix(1:2,1,2))
plot(-log10(p),main="Scatter plot for -log10(p-value)",
     col=c(rep("red",30),rep("black",9970)))
legend("topright",legend = c("NULL","ALTER"),col=c("black","red"),pch=1)
plot(EBP,main="Scatter plot for EBP value",
     col=c(rep("red",30),rep("black",9970)))
legend("bottomright",legend = c("NULL","ALTER"),col=c("black","red"),pch=1)