Coloc: code to simulate a small dataset

Chris Wallace

2021-04-06

vignette: > % % %

Simulate a small dataset

setClass("simdata",
         representation(df1="data.frame",df2="data.frame"))
setValidity("simdata", function(object) {
  n <- nrow(object@df1)
  if(nrow(object@df2)!=n)
    return("nrow of '@df1' should equal nrow of '@df2'")
})

## Class "simdata" [in ".GlobalEnv"]
## 
## Slots:
##                             
## Name:         df1        df2
## Class: data.frame data.frame

setMethod("show", signature="simdata", function(object) {
  cat("pair of simulated datasets, with",ncol(object@df1)-1,"SNPs and",nrow(object@df1),"samples.\n")
})

library(mvtnorm)
simx <- function(nsnps,nsamples,S,maf=0.1) {
    mu <- rep(0,nsnps)
    rawvars <- rmvnorm(n=nsamples, mean=mu, sigma=S)
    pvars <- pnorm(rawvars)
    x <- qbinom(1-pvars, 2, maf)
}


sim.data <- function(nsnps=50,nsamples=200,causals=floor(nsnps/2),nsim=1) {
  cat("Generate",nsim,"small sets of data\n")
  ntotal <- nsnps * nsamples * nsim
  S <- (1 - (abs(outer(1:nsnps,1:nsnps,`-`))/nsnps))^4
  X1 <- simx(nsnps,ntotal,S)
  X2 <- simx(nsnps,ntotal,S)
  Y1 <- rnorm(ntotal,rowSums(X1[,causals,drop=FALSE]/2),2)
  Y2 <- rnorm(ntotal,rowSums(X2[,causals,drop=FALSE]/2),2)
  colnames(X1) <- colnames(X2) <- paste("s",1:nsnps,sep="")
  df1 <- cbind(Y=Y1,X1)
  df2 <- cbind(Y=Y2,X2)
  if(nsim==1) {
    return(new("simdata",
               df1=as.data.frame(df1),
               df2=as.data.frame(df2)))
  } else {
    index <- split(1:(nsamples * nsim), rep(1:nsim, nsamples))
    objects <- lapply(index, function(i) new("simdata", df1=as.data.frame(df1[i,]),
                                             df2=as.data.frame(df2[i,])))
    return(objects)
  }
}

set.seed(46411)
data <- sim.data()

## Generate 1 small sets of data