R/alternatives_normal.R

In jean997/jadeSims:

#'Run alternative methods for normal simulations
#'@description Run alternative t-test based methods on normal data
#'@param file.prefix File prefixe
#'@return Nothing
#'@export
alternatives_normal <- function(file.prefix){
  data.file <- paste("data/", file.prefix, "_data.RData", sep="")
  out.file <- paste("alt/", file.prefix, "_altpvals.RData", sep="")
  #cat(out.file, "\n")
  R <- jadeTF:::getobj(data.file)

  p <- dim(R$Y)[1]
  K <- length(R$sample.size)
  stopifnot(dim(R$Y)[2] == sum(R$sample.size))
  stopifnot(K==2)
  sites <- 1:p

  strt <- c(1)
  stp <- c()
  for(i in 1:K){
    stp <- c(stp, strt[i] + R$sample.size[i] -1)
    if(i < K) strt <- c(strt, strt[i] + R$sample.size[i])
  }


  #No smoothing + T-test
  tt.pvals <- apply(R$Y, MARGIN=1, FUN=ttest.func, sample.size=R$sample.size)

  #Spline + T-Tests
  spline.Y <- matrix(nrow=p, ncol=sum(R$sample.size))
  for(i in 1:K){
    spline.Y[, strt[i]:stp[i]] <- apply(R$Y[, strt[i]:stp[i]], MARGIN=2,
                                        FUN=spline.func, sites=sites)
  }
  spline.pvals <- apply(spline.Y, MARGIN=1, FUN=ttest.func, sample.size=R$sample.size)
  spline.wt.ttests <- fit_withttest(R$Y, R$sample.size, sites, B=100, type="spline")

  #Locfit + T-test
  locfit.Y <- matrix(nrow=p, ncol=sum(R$sample.size))
  for(i in 1:K){
    locfit.Y[, strt[i]:stp[i]] <- apply(R$Y[, strt[i]:stp[i]], MARGIN=2,
                                        FUN=locfit.func, sites=sites)
  }
  locfit.pvals <- apply(locfit.Y, MARGIN=1, FUN=ttest.func, sample.size=R$sample.size)
  locfit.wt.ttests <- fit_withttest(R$Y, R$sample.size, sites, B=100, type="locfit")


  #Adjust all pvalues using SLIM
  tt.slim.pi0 <- SLIMfunc(tt.pvals)
  tt.slim <- QValuesfun(tt.pvals, tt.slim.pi0$pi0_Est)

  spline.tt.slim.pi0 <- SLIMfunc(spline.pvals)
  spline.slim <- QValuesfun(spline.pvals, spline.tt.slim.pi0$pi0_Est)

  locfit.tt.slim.pi0 <- SLIMfunc(locfit.pvals)
  locfit.slim <- QValuesfun(locfit.pvals, locfit.tt.slim.pi0$pi0_Est)

  #Pvals
  stats <- data.frame(tt.pvals, tt.slim,
                      spline.pvals, spline.slim, spline.wt.ttests,
                      locfit.pvals, locfit.slim, locfit.wt.ttests)
  save(stats, file=out.file)
}


locfit.gcv.func <- function(alpha, x, sites){
  z <- gcv(x~sites, alpha=alpha)
  return(z["gcv"])
}

locfit.func <- function(x, sites){
  alpha.opt <- optimize(f=locfit.gcv.func, interval=c(0, 0.1), x=x, sites=sites, maximum=FALSE)
  alpha <- alpha.opt$minimum
  fit <- locfit(x~sites, alpha=alpha)
  pp <- preplot(fit, where="data")
  return(pp$fit)
}

fit_withttest <- function(Y, sample.size, sites, B=100, type=c("spline", "locfit")){
  type <- match.arg(type)
  p <- dim(Y)[1]
  avg <- matrix(0, p , 2)
  var <- matrix(0, p , 2)
  strt <- 1
  for(j in 1:2){
    stp <- strt + sample.size[j] -1
    y <- rowMeans(Y[, strt:stp])
    sd_hat <- apply(Y[, strt:stp], MARGIN=1, FUN=sd)
    fits <- matrix(0, p, B)
    i <- 1
    while(i <=B){
      myy <- rnorm(n=p, mean=y, sd=sd_hat)
      if(type=="spline"){
        fits[,i] <- spline.func(myy, sites)
      }else if(type=="locfit"){
        fits[,i] <- locfit.func(myy, sites)
      }
      i <- i+1
    }
    avg[,j] <- rowMeans(fits)
    var[,j]<- (1/(B-1))*rowSums(apply(fits, MARGIN=2, FUN=function(x, avg){
      return((x-avg)^2)
    }, avg=avg[,j]))
  }
  ttests <- unlist(lapply(1:p, FUN=function(idx, avg, var){
    t <- (avg[idx,1]-avg[idx,2])/sqrt(var[idx,1] + var[idx,2])
    return(t)}, avg=avg, var=var))
  return(ttests)

}