require(MASS)

# Variables to change for each run - Build the simulation design matrix

subth <- c(.5, 1, 2)
prev  <- c(.15, .10, .01)
Maf   <- c(.1, .05, .01)
CORR  <- c( .05, .025, .1)

design <- cbind(
  rep(CORR , each = length(subth)*length(prev)*length(Maf) ),
  rep(Maf , each = length(subth)*length(prev) ),
  rep(prev, each = length(subth)), 
  rep(subth, length(prev) )
)

colnames(design) <- c("eff", "maf", "prev", "sub")

# Sample Size
N <- 400000

# Effect Size
R <- design[,1] 

# Major allele frequency
maf = design[,2]

#Prevalence for case definition
prevc <- design[,3]

#Prevalence for subthreshold
prevs <- design[,4]

# Number of iterations
iter <- 1000


# Create empty items to store data from simulations

ordinal_coef <- matrix(NA, iter, 3)
colnames(ordinal_coef) <- c("beta", "SE", "t-value")

super_coef   <- matrix(NA, iter, 4)
colnames(super_coef) <- c("beta", "SE", "Z", "p-value")

cc_coef      <- matrix(NA, iter, 4)
colnames(cc_coef) <- c("beta", "SE", "Z", "p-value")

Results <- list()

for(n in 1:nrow(design)){
  # Effect Size
  R <- design[n,1] 
  
  # Major allele frequency
  maf = design[n,2]
  
  #Prevalence for case definition
  prevc <- design[n,3]
  
  #Prevalence for subthreshold
  prevs <- design[n,4]*prevc
  
  params <- (c(N, R, maf, prevc, prevs, iter))
  names(params) <- c("N", "R", "MAF", "prevc", "prevs", "iterations")
  
  # Set up the loop to run the simulations
  for(i in 1:iter){
    
    sigma <- matrix(c(1,R,R,1),2,2)
    dat <- as.data.frame(mvrnorm(N, c(0,0), sigma, empirical = F))
    colnames(dat) <- c("Y", "X")
    
    # Want to break it up into general SNP and Phenotype categories
    #### Setting up the Phenotype
    
    #Case threshold
    threshc <- qnorm(1 - prevc) 
    
    #Not case threshold
    threshs <- qnorm(1 - (prevs + prevc)) 
    
    # Create the phenotype item in the data frame
    dat$dv <- cut(dat$Y, c(-Inf, threshc, threshs, Inf), labels = 0:2)
    
    
    #### Setting up the SNP
    
    # This part will ensure that everything is in Hardy Weinberg Equilibrium
    
    p = maf
    
    q <- 1-p
    
    homo_m <- p^2
    homo_M <- q^2
    het <- 2*p*q
    
    # Set the thresholds for Major and Minor alleles
    thresh_M <- qnorm(1 - homo_M) 
    thresh_m <- qnorm(homo_m) 
    
    # Create the SNP item in the dataframe (homozygous minor, heterozygous, homozygous major)
    dat$iv <- cut(dat$X, c(-Inf, thresh_m,  thresh_M, Inf), labels = 0:2)
    
    
    # Set up super-control
    dat$super <- dat$dv
    dat$super[dat$super == 1] <- NA
    dat$super[dat$super == 2] <- 1
    
    #Set up standard case-control
    dat$cc <- dat$dv
    dat$cc[dat$cc == 1] <- 0
    dat$cc[dat$cc == 2] <- 1
    
    
    # Run the probit models with the simulated data
    
    ord <- polr(dv ~ as.numeric(iv), method = "probit", data = dat)
    
    super <- glm(super ~ as.numeric(iv), family = binomial(link = "probit"), data = dat)
    
    cc <- glm(cc ~ as.numeric(iv), family = binomial(link = "probit"), data = dat)
    
    # Take the coefficients
    ordinal_coef[i,] <- coef(summary(ord))[1,]
    super_coef[i,]   <- coef(summary(super))[2,]
    cc_coef[i,]      <- coef(summary(cc))[2,]
    
     
    # Create the final list to store all of the results in with the effect size and maf specified
    simResults <- list(params, ordinal_coef, super_coef, cc_coef, correlationXY)
    names(simResults) <- c("parameters", "ordinal_coef", "super_coef", "cc_coef")
    
    
  }
  Results[[n]] <-simResults

  print(n)
  
}


###############################################################################

# The mean of each set of parameters was used in figures

## Create each item matrix to take the mean for each combination of parameters

ordinal_cfnts <- matrix(NA, nrow(design), 3)

colnames(ordinal_cfnts) <- c("Beta", "SE", "t")


super_cfnts <- matrix(NA, nrow(design), 3)

colnames(super_cfnts) <- c("Beta", "SE", "Z")


cc_cfnts <- matrix(NA, nrow(design), 3)

colnames(cc_cfnts) <- c("Beta", "SE", "Z")

## Loop through the list of results, taking the mean, and storing it in the newly created matrices

for(n in 1:nrow(design)){

  ordinal_cfnts[n, 1] <- mean(Results[[n]]$ordinal_coef[,1])
  
  ordinal_cfnts[n, 2] <- mean(Results[[n]]$ordinal_coef[,2])
  
  ordinal_cfnts[n, 3] <- mean(Results[[n]]$ordinal_coef[,3])
  
  
  super_cfnts[n, 1] <- mean(Results[[n]]$super_coef[,1])
  
  super_cfnts[n, 2] <- mean(Results[[n]]$super_coef[,2])
  
  super_cfnts[n, 3] <- mean(Results[[n]]$super_coef[,3])
  
  
  cc_cfnts[n, 1] <- mean(Results[[n]]$cc_coef[,1])
  
  cc_cfnts[n, 2] <- mean(Results[[n]]$cc_coef[,2])
  
  cc_cfnts[n, 3] <- mean(Results[[n]]$cc_coef[,3])
  
  
  coefs <- cbind(ordinal_cfnts, super_cfnts, cc_cfnts)
  colnames(coefs) <- c("ordinal_beta", "ordinal_SE", "ordinal_t", "super_beta", "super_SE", "super_Z",
                       "cc_beta", "cc_SE", "cc_Z")
  
}

# Bind the matrices for each control coding scheme into a single dataframe with the parameters

FullSimResults <- cbind(design, coefs)
colnames(FullSimResults) <- c("eff",  "maf",  "prev", "sub", "ordinal_beta", "ordinal_SE", "ordinal_t", "super_beta", "super_SE", "super_Z",
                     "cc_beta", "cc_SE", "cc_Z")

FullSimResults <- as.data.frame(FullSimResults)





# MAFs - .05, .1, .01, .5

# Case Prevalences
## 15% for depression (0.15)
## 10% for impulse control disorders (0.1) 
## 1% for Schizophrenia (0.01)

# Increase subthreshold prevalence
## half the case prevalence
## equal to case prevalence
## double case prevalence

# Impact of the correlation (R)
## .05  
## .025  
## .1