SimulationsTable3.R

# Simulation for the QREM+SEMMS paper
# see sims.R for simulation configurations.

# compare with mean-model (SEMMS)
# plots - after simulations, quantile functions

rm(list=ls())
library("SEMMS")
library("QREM")
source("sims.R")
nn <- 4
maxRep <- 40
plots <- FALSE
simsToRun <- 8 # see sims.R
res <- matrix(0,ncol=7,nrow=100*9*length(simsToRun))

qn <- 0.8
simno = 1
i=1

confmat <- function(allidx, actual, estimated) {
  TP <- which(estimated  %in% actual)
  FP <- which(estimated %in% setdiff(allidx,actual))
  TN <- which(setdiff(allidx,estimated) %in% setdiff(allidx,actual))
  FN <- which(setdiff(allidx,estimated) %in% actual)
  c(length(TP), length(FP),length(TN), length(FN))
}

runSEMMSandQREM <- function(filename, qn) {
  dataYXZ <- readInputFile(filename, ycol=1, Zcols=2:501)
  n <- dataYXZ$N
  K <- dataYXZ$K
  t0=Sys.time()
  cat("initializing...\n")
  pval <- rep(0, K)  ### need to speed up pval initialization!
  for (i in 1:K) {
    linmod <- as.formula(paste("Y ~", colnames(dataYXZ$Z)[i]))
    dframetmp <- data.frame(cbind(dataYXZ$Y, dataYXZ$Z[,i]))
    colnames(dframetmp) <- c("Y",colnames(dataYXZ$Z)[i])
    qremFit <- QREM(lm,linmod, dframetmp, qn, maxInvLambda = 1000)
    sm <- summary(qremFit$fitted.mod)
    pval[i] <- sm$coefficients[2,4]
  }
  nnset = order(pval)[1:nn]
  t1=Sys.time()
  cat(round(difftime(t1,t0, units = "secs")),"seconds. initial set", nnset,"\n")
  for (rep in 1:maxRep) {
    # create a subset of the selected columns and run QREM
    preds <- paste(colnames(dataYXZ$Z)[nnset], collapse = "+")
    linmod <- as.formula(paste("Y ~", preds))
    dframetmp <- data.frame(cbind(dataYXZ$Y, dataYXZ$Z[,nnset]))
    colnames(dframetmp) <- c("Y",colnames(dataYXZ$Z)[nnset])
    qremFit <- QREM(lm,linmod, dframetmp, qn, maxInvLambda = 1000)
    # apply the weights found by QREM and rerun SEMMS
    dataYXZtmp <- dataYXZ
    dataYXZtmp$Y <- (dataYXZ$Y-(1-2*qn)/qremFit$weights)
    fittedVSnew <- fitSEMMS(dataYXZtmp,distribution = 'N', mincor=0.8,rnd=F,
                            nnset=nnset, minchange = 1, maxst = 20)
    if (length(fittedVSnew$gam.out$nn) == 0) {
      return(fittedVSnew$gam.out$nn)
    }
    #foundSEMMSnew <- sort(union(which(fittedVSnew$gam.out$lockedOut != 0),
    #                            fittedVSnew$gam.out$nn))
    t2=Sys.time()
    cat(round(difftime(t2,t1,units = "secs")),"seconds.",rep,":\t",fittedVSnew$gam.out$nn,":\t", qremFit$empq,"\n")
    t1=t2
    if (length(fittedVSnew$gam.out$nn) == length(nnset)) {
      if (all(fittedVSnew$gam.out$nn == nnset)) {
        break
      }
    }
    nnset <- fittedVSnew$gam.out$nn
  }
  if (plots) {
    fittedGLM <- runLinearModel(dataYXZtmp,nnset, "N")
    print(summary(fittedGLM$mod))
    plotMDS(dataYXZ, fittedVSnew, fittedGLM, ttl="...")
    plotFit(fittedGLM)
    plot(dataYXZ$Y, col=(2+(qremFit$ui>0)), cex=0.7, pch=19)
    for (i in 2:ncol(dframetmp)) {
      qrdiag <- QRdiagnostics(dframetmp[,i],colnames(dframetmp)[i],qremFit$ui,qn)
    }
  }
  t3=Sys.time()
  cat(round(difftime(t3,t0,units = "secs")),"seconds (total). Done\n")
  nnset
}

cnt <- 1
for (simno in 1:length(sims)) {
  if (simno %in% simsToRun) {
    cat("\nSim. #",simno,"\n")
    n <- sims[[simno]]$n
    qns <- sims[[simno]]$qns
    reps <- sims[[simno]]$reps
    coefs <- sims[[simno]]$coefs
    lp <- as.list(attr(terms(sims[[simno]]$mod), "variables"))[-(1:2)]
    truepreds = rep(0,length(lp))
    for(i in 1:length(lp)) {
      truepreds[i] = gsub("[a-zA-Z]","",lp[[i]], perl=TRUE)
    }
    truepreds <- as.numeric(truepreds)
    
    xrng <- matrix(sims[[simno]]$xrng,  nrow=(length(coefs)-1), ncol=2)
    for (i in 1:reps) {
      set.seed(simno*100000 + sims[[simno]]$seed + i)
      X <- matrix(0,nrow=n, ncol=length(coefs))
      X[,1] <- rep(1,n)
      for (jj in 1:(length(coefs)-1)) {
        X[,jj+1] <- runif(n, min=xrng[jj,1], max=xrng[jj,2])
      }
      colnames(X) <- c("const", paste("X",1:(length(coefs)-1),sep=""))
      if (sims[[simno]]$errvar == 0) {
        errs <- sims[[simno]]$err(rep(0,n))
      } else {
        if (length(sims[[simno]]$errvar) == 2) {
          errs <- sims[[simno]]$err(X[,sims[[simno]]$errvar[1]+1], X[,sims[[simno]]$errvar[2]+1])
        } else {
          errs <- sims[[simno]]$err(X[,sims[[simno]]$errvar+1])
        }
      }
      y <- X%*%coefs + errs
      dframe <- data.frame(y,X[,-1], matrix(rnorm((500-ncol(X)+1)*n,0,0.1), 
                                            nrow=n, ncol=(500-ncol(X)+1)))
      datfn <- sprintf("data/sim%02d.RData",simno)
      save(dframe, file=datfn)
      for (qn in qns) {
        selected <- runSEMMSandQREM(datfn, qn)
        res[cnt,]  <- c(simno, i, qn, confmat(1:500, truepreds, selected))
        cat(c(simno, i, qn, confmat(1:500, truepreds, selected)),"\n")
        cnt <- cnt+1
      }
    }
    save(res,file=sprintf("results/resSim%02d.RData",simno))
  }
}