R/gSim.R

In sambiostat/WAPL: What the Package Does (One Line, Title Case)

#'simulate examples for WAPL
#'This is the generative function for WAPL
#'@param scenario: takes value 1-4, represent different scenarios
#'@param N: number of subjects
#'@param sigma: value of noise
#'@export
#'@return a list contains generated data and related information

gSim<-function(scenario=c(1,2,3,4), N=200, sigma=0.2){
  #p=50 predictors in 20 blocks of 3; the later 40 have all coef of zero
  #within group the data are uniform with correlation 0.2


  A=sample(c(1,-1), N, replace=T, prob=c(0.5,0.5))

  if(scenario==1){
    X = Reduce(cbind, lapply(1:50, function(...){
      ##a block of uniform variable, within group correlation is 0.2
      matrix(runif(N*2, 0, 1), nrow=N) + matrix(rep(runif(N, 0, 1), 2), nrow=N, byrow=F)
    }))
    colnames(X)=paste("X",c(1:100),sep="")
    Ym= 2*X[,3]+ 2*pi*(sin(pi*X[,4]))
    Yc= 3*pi*(sin(pi*X[,1]*2)) - 0.5*(X[,2]-0.5)^2
  }else if(scenario==2){
    X = Reduce(cbind, lapply(1:10, function(...){
      ##a block of uniform variable, within group correlation is 0.2
      matrix(runif(N*2, 0, 1), nrow=N) + matrix(rep(runif(N, 0, 1), 2), nrow=N, byrow=F)
    }))
    colnames(X)=paste("X",c(1:20),sep="")
    Ym= 2*X[,3]+ 2*pi*(sin(pi*X[,4]))
    Yc= 3*pi*(sin(pi*X[,1]*2)) - 0.5*(X[,2]-0.5)^2
  }else if(scenario==3){
    X = Reduce(cbind, lapply(1:50, function(...){
      ##a block of uniform variable, within group correlation is 0.2
      matrix(runif(N*2, 0, 1), nrow=N) + matrix(rep(runif(N, 0, 1), 2), nrow=N, byrow=F)
    }))

    colnames(X)=paste("X",c(1:100),sep="")
    Ym= 4*X[,5]+ 4*pi*(sin(pi*X[,6]))
    Yc= 5*pi*(sin(pi*X[,1]*2)) - 1*(X[,2]-0.5)^2 +  10*(sin(pi*X[,3])) + 2*cos(pi*X[,4]*2)
  } else if(scenario==4){
    X = Reduce(cbind, lapply(1:10, function(...){
      ##a block of uniform variable, within group correlation is 0.2
      matrix(runif(N*2, 0, 1), nrow=N) + matrix(rep(runif(N, 0, 1), 2), nrow=N, byrow=F)
    }))

    colnames(X)=paste("X",c(1:20),sep="")
    Ym= 4*X[,5]+ 4*pi*(sin(pi*X[,6]))
    Yc= 5*pi*(sin(pi*X[,1]*2)) - 1*(X[,2]-0.5)^2 +  10*(sin(pi*X[,3])) + 2*cos(pi*X[,4]*2)
  } else {break;}

  noise=rnorm(N, sd=sigma)
  s2n=var(Yc) /(var(Ym))   #+ var(noise))
  Y=Ym + A*Yc + noise

  optA={Yc > 0}*2-1
  Y0=Ym - Yc #Everyone follows A=-1
  Y1=Ym + Yc #Everyone follows A=1
  optY=Ym + abs(Yc)

  tmp=table(optA)
  out=list("X"=X, "A"=A, "Y"=Y, "optA"=optA,  "Y1"=Y1, "Y0"=Y0, "optY"=optY,
           "s2n"=s2n, "PrOptA1"= tmp[2]/N, "EY"=mean(optY), "EYopt0"=mean(optY[optA==-1]),
           "EYopt1"=mean(optY[optA==1]), "EYAll0"=mean(Y0), "EYAll1"=mean(Y1),
           "MedY"=median(optY), "MedYopt0"=median(optY[optA==-1]),
           "MedYopt1"=median(optY[optA==1]), "MedYAll0"=median(Y0), "MedYAll1"=median(Y1))
  #return(unlist(out[8:17]))
  return(out)
}