SimulationAR1.R

rm(list=ls())
library("SEMMS")
library("QREM")
library("rqPen")
library("MASS")

nn <- 5
maxRep <- 40
plots <- FALSE


simulateAR1 <- function(P, N, corcoef=0.9) {
  Sigma <- toeplitz(1:P)
  for (i in 2:P) {
    Sigma[which(Sigma == i)] <- corcoef^(i-1)
  }
  mvrnorm(N, rep(0,P), Sigma)
}

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))
}

checkLoss <- function(ui,qn) {
  ui*(qn - as.numeric(ui < 0))
}

runSEMMSandQREM <- function(filename, qn, nnset=NULL) {
  dataYXZ <- readInputFile(filename, ycol=1, Zcols=2:1001)
  n <- dataYXZ$N
  K <- dataYXZ$K
  t0=Sys.time()
  if (is.null(nnset)) {
    cat("initializing...\n")
    zval <- rep(0, K)
    rnd <- sample(K, replace = FALSE)
    m <- 5
    for (i in 1:(1000/m)) {
      idx <- ((i-1)*m+1) : (i*m)
      preds <- paste0(colnames(dataYXZ$Z)[rnd[idx]], collapse = " + ")
      linmod <- as.formula(paste("Y ~", preds))
      dframetmp <- data.frame(cbind(dataYXZ$Y, dataYXZ$Z[,rnd[idx]]))
      colnames(dframetmp) <- c("Y",colnames(dataYXZ$Z)[rnd[idx]])
      qremFit <- QREM(lm,linmod, dframetmp, qn, maxInvLambda = 1000)
      zval[rnd[idx]] <- qremFit$coef$beta[-1]/sqrt(diag(bcov(qremFit,linmod,dframetmp,qn)))[-1]
    }
    nnset <- order(abs(zval),decreasing = TRUE)[1:nn]
  }
  
  t1=Sys.time()
  cat(round(difftime(t1,t0, units = "secs")),"seconds. initial set", nnset,"\n")
  ll <- 0
  for (repno 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)
    newll <- sum(checkLoss(qremFit$ui,qn))
    if (abs(ll - newll) < 1e-3) {
      nnset <- sort(c(fittedVSnew$gam.out$nn, which(fittedVSnew$gam.out$lockedOut == 1 )))
      cat(repno,"steps\n")
      return(nnset)
    }
    ll <- newll
    # 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)
    }
    t2=Sys.time()
    cat(round(difftime(t2,t1,units = "secs")),"seconds.",repno,":\t",fittedVSnew$gam.out$nn,":\t", qremFit$empq,"\n")
    t1=t2
    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)
    }
  }
  nnset <- sort(c(fittedVSnew$gam.out$nn, which(fittedVSnew$gam.out$lockedOut == 1 )))
  t3=Sys.time()
  cat(round(difftime(t3,t0,units = "secs")),"seconds (total). Done\n")
  nnset
}

res <- matrix(0,ncol=11,nrow=100*9)
cnt <- 1
simno = 9
qns = seq(0.1,0.9,by=0.1)
for (i in 1:100) {
  M0 = simulateAR1(980,100,0)
  M = simulateAR1(20,100, corcoef=0.95)
  truepreds = 1:20
  X = cbind(M,M0)
  y = rowSums(X[,1:20]+rnorm(100,0,0.1))
  dframe <- as.data.frame(cbind(y,X))
  datfn <- "simAR.RData"
  save(dframe, file=datfn)
  for (qn in qns) {
    L1fit <- rq.lasso.fit(X,y,lambda=0.13,tau = qn)
    selectedL1 <- which(L1fit$coefficients[-1]!=0)
    if (length(selectedL1) > 2) {
      selected <- runSEMMSandQREM(datfn, qn, nnset=selectedL1)
    } else {
      selected <- runSEMMSandQREM(datfn, qn)
    }
    res[cnt,]  <- c(simno, i, qn, confmat(1:1000, truepreds, selected),
                    confmat(1:1000, truepreds, selectedL1))
    cat(c(simno, i, qn, confmat(1:1000, truepreds, selected), 
          confmat(1:1000, truepreds, selectedL1)),"\n")
    cnt <- cnt+1
  }
}

save(res,file="resSim.RData")