
Sys.time()
## ---------------------------------------------
## R codes that are used or cited in the paper submitted to NEJSDS
## ---------------------------------------------
PathAll <- "/Users/farhad/Library/CloudStorage/Dropbox/Apps/Overleaf/NEJoDS_fmrs/Codes/"
PathData <- "/Users/farhad/Library/CloudStorage/Dropbox/Apps/Overleaf/NEJoDS_fmrs/Codes/PythonDataGen/"
#PathAll <- "..//Codes/"
#PathData <- "..//Codes/PythonDataGen/"
setwd(PathAll)
source(paste0(PathAll, "/fmrlasso.r"))
#source("..//fmrlasso.r")
library(tidyverse)
library(kableExtra)
library(knitr)
library(fmrs)

## ---------------------------------------------
## FMAFTR model with log-normal and component-wise tuning parameters
## ---------------------------------------------
library(fmrs); set.seed(1980)
myk <- 2; myd <- 10; myn <- 500; myrho <- 0.25
mydisp <- c(1,1); myu <- 40; mypi <- c(0.4,0.6)
coeff1 <- c(2,2,-1,-2,1,2,0,0,0,0,0)
coeff2 <- c(-1,-1,1,2,0,0,0,0,-1,2,-2)
## ---------------------------------------------
dat <- 
  list(y = as.numeric(unlist(read.csv(
    paste0(PathData, 'Ydata_1.csv'), header = TRUE))),
    x = as.matrix(read.csv(
      paste0(PathData, 'Xdata_1.csv'), header = TRUE)),
    delta = as.numeric(unlist(read.csv(
      paste0(PathData, 'Ddata_1.csv'), header = TRUE))))
## ---------------------------------------------
res.mle <- 
  fmrs.mle(y = dat$y, delta = dat$delta, 
           x = dat$x, nComp = myk, 
           disFamily = "lnorm", 
           initCoeff = rnorm(myk * myd + myk),
           initDispersion = rep(1, myk), 
           initmixProp = rep(1 / myk, myk))
summary(res.mle)
## ---------------------------------------------
res.lam <- 
  fmrs.tunsel(y = dat$y, delta = dat$delta, 
              x = dat$x, nComp = myk, 
              disFamily = "lnorm", 
              initCoeff = c(coefficients(res.mle)),
              initDispersion = dispersion(res.mle), 
              initmixProp = mixProp(res.mle), 
              penFamily = "adplasso", LambMin = 0.01, 
              LambMax = 1, nLamb = 100)
summary(res.lam)
## ---------------------------------------------
res.var <- 
  fmrs.varsel(y = dat$y, delta = dat$delta, 
              x = dat$x, nComp = myk, 
              disFamily = "lnorm", 
              initCoeff = c(coefficients(res.mle)),
              initDispersion = dispersion(res.mle), 
              initmixProp = mixProp(res.mle), 
              penFamily = "adplasso", 
              lambPen = slot(res.lam,"lambPen"))
summary(res.var)
## ---------------------------------------------
## FMAFTR model with log-normal and common tuning parameters
## ---------------------------------------------
set.seed(1980)
res.mle <- 
  fmrs.mle(y = dat$y, delta = dat$delta, 
           x = dat$x, nComp = myk, 
           disFamily = "lnorm", 
           initCoeff = rnorm(myk * myd + myk),
           initDispersion = rep(1, myk), 
           initmixProp = rep(1 / myk, myk))
summary(res.mle)
## ---------------------------------------------
BICmin <- 1e+22
optLam <- 0
Ctune <- seq(0.01, 1, length.out = 100)
for (j in seq_along(Ctune)) {
  holdfit <- 
    fmrs.varsel(y = dat$y, delta = dat$delta, 
                x = dat$x, nComp = myk, 
                disFamily = "lnorm", 
                initCoeff = c(coefficients(res.mle)),
                initDispersion = dispersion(res.mle), 
                initmixProp = mixProp(res.mle), 
                penFamily = "adplasso", 
                lambPen = rep(Ctune[j], myk))
  if (BIC(holdfit) < BICmin) {
    BICmin <- BIC(holdfit)
    optLam <- Ctune[j]
    res.var <- holdfit
  }
}
summary(res.var)
optLam

## ---------------------------------------------
## FMR model with normal and component-wise tuning parameters
## ---------------------------------------------
## ---------------------------------------------
set.seed(1980) 
dat <- 
  list(y = as.numeric(unlist(read.csv(
    paste0(PathData, 'Ydata_2.csv'), header = TRUE))),
    x = as.matrix(read.csv(
      paste0(PathData, 'Xdata_2.csv'), header = TRUE)))

## ---------------------------------------------
res.mle <- 
  fmrs.mle(y = dat$y, x = dat$x, nComp = myk, 
           disFamily = "norm",
           initCoeff = rnorm(myk * myd + myk), 
           initDispersion = rep(1,myk), 
           initmixProp = rep(1 / myk, myk))
summary(res.mle)
## ---------------------------------------------
res.lam <- 
  fmrs.tunsel(y = dat$y, x = dat$x, nComp = myk, 
              disFamily = "norm",
              initCoeff = c(coefficients(res.mle)), 
              initDispersion = dispersion(res.mle), 
              initmixProp = mixProp(res.mle), 
              penFamily = "adplasso", LambMin = 0.01, 
              LambMax = 1, nLamb = 100)
summary(res.lam)
## ---------------------------------------------
res.var <- 
  fmrs.varsel(y = dat$y, x = dat$x, nComp = myk, 
              disFamily = "norm", 
              initCoeff = c(coefficients(res.mle)),
              initDispersion = dispersion(res.mle),
              initmixProp = mixProp(res.mle), 
              penFamily = "adplasso",
              lambPen = slot(res.lam, "lambPen"))
summary(res.var)
## ---------------------------------------------
## FMR model with normal and common tuning parameters
## ---------------------------------------------
set.seed(1980)
res.mle <- 
  fmrs.mle(y = dat$y, x = dat$x, nComp = myk, 
           disFamily = "norm",
           initCoeff = rnorm(myk * myd + myk), 
           initDispersion = rep(1,myk), 
           initmixProp = rep(1 / myk, myk))
summary(res.mle)
## ---------------------------------------------
BICmin <- 1e+22
optLam <- 0
Ctune <- seq(0.01, 1, length.out = 100)
for (j in seq_along(Ctune)) {
  holdfit <- 
    fmrs.varsel(y = dat$y, x = dat$x, 
                nComp = myk, disFamily = "norm", 
                initCoeff = c(coefficients(res.mle)),
                initDispersion = dispersion(res.mle),
                initmixProp = mixProp(res.mle), 
                penFamily = "adplasso", 
                lambPen = rep(Ctune[j], myk))
  if (BIC(holdfit) < BICmin) {
    BICmin <- BIC(holdfit)
    optLam <- Ctune[j]
    res.var <- holdfit
  }
  
}
summary(res.var)
optLam

## ---------------------------------------------
## FMR model with large sample size
## ---------------------------------------------
library(fmrs); set.seed(1980)
myk <- 2; myd <- 10; myn <- 20000; myrho <- 0.25
mydisp <- c(1,1); myu <- 40; mypi <- c(0.4,0.6)
coeff1 <- c(2,2,-1,-2,1,2,0,0,0,0,0)
coeff2 <- c(-1,-1,1,2,0,0,0,0,-1,2,-2)
## ---------------------------------------------
dat <- 
  list(y = as.numeric(unlist(read.csv(
    paste0(PathData, 'Ydata_3.csv'), header = TRUE))),
    x = as.matrix(read.csv(
      paste0(PathData, 'Xdata_3.csv'), header = TRUE)))

## ---------------------------------------------
res.mle <- 
  fmrs.mle(y = dat$y, x = dat$x, nComp = myk, 
           disFamily = "norm",
           initCoeff = rnorm(myk * myd + myk), 
           initDispersion = rep(1,myk), 
           initmixProp = rep(1 / myk, myk))
summary(res.mle)
## ---------------------------------------------
res.lam <- 
  fmrs.tunsel(y = dat$y, x = dat$x, nComp = myk, 
              disFamily = "norm",
              initCoeff = c(coefficients(res.mle)), 
              initDispersion = dispersion(res.mle), 
              initmixProp = mixProp(res.mle), 
              penFamily = "adplasso", LambMin = 0.01, 
              LambMax = 1, nLamb = 100)
summary(res.lam)
## ---------------------------------------------
res.var <- 
  fmrs.varsel(y = dat$y, x = dat$x, nComp = myk, 
              disFamily = "norm", 
              initCoeff = c(coefficients(res.mle)),
              initDispersion = dispersion(res.mle),
              initmixProp = mixProp(res.mle), 
              penFamily = "adplasso",
              lambPen = slot(res.lam, "lambPen"))
summary(res.var)


## ---------------------------------------------
## Non-mixture AFT model with log-normal
## ---------------------------------------------
library(fmrs); set.seed(1980)
myk <- 1; myd <- 10; myn <- 500; myrho <- 0.25
mydisp <- c(1); myu <- 40; mypi <- c(1.0)
coeff1 <- c(2,2,-1,-2,1,2,0,0,0,0,0)
## ---------------------------------------------
dat <- 
  list(y = as.numeric(unlist(read.csv(
    paste0(PathData, 'Ydata_4.csv'), header = TRUE))),
    x = as.matrix(read.csv(
      paste0(PathData, 'Xdata_4.csv'), header = TRUE)),
    delta = as.numeric(unlist(read.csv(
      paste0(PathData, 'Ddata_4.csv'), header = TRUE))))
## ---------------------------------------------
res.mle <- 
  fmrs.mle(y = dat$y, delta = dat$delta, 
           x = dat$x, nComp = myk, 
           disFamily = "lnorm", 
           initCoeff = rnorm(myk * myd + myk),
           initDispersion = rep(1, myk), 
           initmixProp = rep(1 / myk, myk))
summary(res.mle)
## ---------------------------------------------
res.lam <- 
  fmrs.tunsel(y = dat$y, delta = dat$delta, 
              x = dat$x, nComp = myk, 
              disFamily = "lnorm", 
              initCoeff = c(coefficients(res.mle)),
              initDispersion = dispersion(res.mle), 
              initmixProp = mixProp(res.mle), 
              penFamily = "adplasso", LambMin = 0.01, 
              LambMax = 1, nLamb = 100)
summary(res.lam)
## ---------------------------------------------
res.var <- 
  fmrs.varsel(y = dat$y, delta = dat$delta, 
              x = dat$x, nComp = myk, 
              disFamily = "lnorm", 
              initCoeff = c(coefficients(res.mle)),
              initDispersion = dispersion(res.mle), 
              initmixProp = mixProp(res.mle), 
              penFamily = "adplasso", 
              lambPen = slot(res.lam,"lambPen"))
summary(res.var)

## ---------------------------------------------
## Simulation study by comparing fmrs and fmrlasso methods
## ---------------------------------------------
library(lubridate)
set.seed(1980); Rep <- 100; 
K <- 2; n <- 500; d <- 5; mpi <- c(0.4,0.6)
disp <- c(1, 1); mrho <- 0.25
coeff1 <- c(-1,2,0,0,-1,2)
coeff2 <- c(2,1,1,0,0,0) 
datalist <- list()
beta.fmrs <- matrix(NA, Rep, K * d)
sel.fmrs <- matrix(NA, Rep, K * d)
beta.fmrl <- matrix(NA, Rep, K * d)
sel.fmrl <- matrix(NA, Rep, K * d)

## ---------------------------------------------
for (r in seq_len(Rep)) {
  datalist[[r]] <-
    list(y = as.numeric(unlist(read.csv(
      paste0(PathData, 'YdataA_', r, '.csv'), header = TRUE))),
      x = as.matrix(read.csv(
        paste0(PathData, 'XdataA_', r, '.csv'), header = TRUE)))
}

sleep_for_a_minute <- function() {Sys.sleep(60)}
fmrs_str_time <- Sys.time()

for (r in seq_len(Rep)) {
  res.mle <- 
    fmrs.mle(y = datalist[[r]]$y, 
             x = datalist[[r]]$x, 
             nComp = K, disFamily = "norm", 
             initCoeff = rnorm(K * d + K),
             initDispersion = rep(1,K), 
             initmixProp = rep(1/K,K))
  res.lam <- 
    fmrs.tunsel(y = datalist[[r]]$y, 
                x = datalist[[r]]$x, 
                nComp = K, disFamily = "norm",
                initCoeff = c(coefficients(res.mle)), 
                initDispersion = dispersion(res.mle), 
                initmixProp = mixProp(res.mle), 
                penFamily = "lasso", LambMin = 0.01, 
                LambMax = 1, nLamb = 100)
  res.var <- 
    fmrs.varsel(y = datalist[[r]]$y, 
                x = datalist[[r]]$x, 
                nComp = K, disFamily = "norm", 
                initCoeff = c(coefficients(res.mle)),
                initDispersion = dispersion(res.mle),
                initmixProp = mixProp(res.mle), 
                penFamily = "lasso",
                lambPen = slot(res.lam, "lambPen"))
  coefficients(res.var)
  beta.fmrs[r,] <- 
    c(coefficients(res.var)[-1, order(
      coefficients(res.var)[1,])])
  sel.fmrs[r,] <- 
    c(res.var@selectedset[, order(
      coefficients(res.var)[1,])])
}

fmrs_end_time <- Sys.time()
fmrs.time <- fmrs_end_time - fmrs_str_time
fmrs.time <- 
  difftime(fmrs_end_time,fmrs_str_time,units = "secs")

fmrlassoBIC_str_time <- Sys.time()

for (r in seq_len(Rep)) {
  init <- ini.ex(k = K,n = n)
  y <- datalist[[r]]$y
  x <- cbind(1, datalist[[r]]$x)
  colnames(x)[1] <- "Int"
  la <- seq(0.01, lambdamax(y, x), length = 100)
  fit2 <- 
    fmrlassopath(x, y, k = K, lambda = la, 
                 ssd.ini = 1, ex.ini = init)
  fit2 <- 
    fmrlasso(x, y, k = 2, 
             lambda = la[which.min(fit2$bic)],
             ssd.ini = 1, ex.ini = init)
  beta.fmrl[r,] <- 
    c(fit2$coef[-1, order(fit2$coef[1,])])
  sel.fmrl[r,] <- 1 * (beta.fmrl[r,] != 0.0)
}

fmrlassoBIC_end_time <- Sys.time()
fmrlasso.time <- 
  fmrlassoBIC_end_time - fmrlassoBIC_str_time
fmrlasso.time <- 
  difftime(fmrlassoBIC_end_time,
           fmrlassoBIC_str_time,
           units = "secs")

IDv <- (c(coeff2[-1], coeff1[-1]) == 0)*1
sel.fmrs.f <- colMeans(sel.fmrs) * 100
sel.fmrs.f[IDv == 1] <- 100 - sel.fmrs.f[IDv == 1]
sel.fmrl.f <- colMeans(sel.fmrl) * 100
sel.fmrl.f[IDv == 1] <- 100 - sel.fmrl.f[IDv == 1]
names(sel.fmrs.f) <- 
  paste0(rep(paste0("beta",1:d), K), 
         rep(1:K ,each = d))
names(sel.fmrl.f) <- 
  paste0(rep(paste0("beta",1:d), K), 
         rep(1:K ,each = d))
sel <- 
  data.frame(fmrs = sel.fmrs.f, 
             fmrlasso = sel.fmrl.f)
sel
L2fmrs <- apply(beta.fmrs, 1, function(x){
  sqrt(sum((x[1:d] - coeff2[-1])^2)) + 
    sqrt(sum((x[(d + 1):(K*d)] - coeff1[-1])^2))
})
L2fmrs <- sum(L2fmrs)
L2fmrs
L2fmrl <- apply(beta.fmrl, 1, function(x){
  sqrt(sum((x[1:d] - coeff2[-1])^2)) + 
    sqrt(sum((x[(d + 1):(K * d)] - coeff1[-1])^2))
})
L2fmrl <- sum(L2fmrl)
L2fmrl

## ---------------------------------------------
## Lung cancer data analysis using fmrs 
## ---------------------------------------------
library(survival)

lung.r <- lung[,-c(1,9)]
r.lung <- lung.r[-which(is.na(rowSums(lung.r))), ]
d <- 6
r.lung$sex <- as.numeric(r.lung$sex - 1)
Days <- r.lung$time
Status <- r.lung$status - 1 
Features <- as.matrix(r.lung[, -c(1, 2)])
Tol <- 0.001
Final.Model <- list()

for (K in 2:5) {
  set.seed(1980)
  res.mle <- 
    fmrs.mle(y = Days, x = Features, 
             delta = Status,
             nComp = K, disFamily = 'lnorm',
             initCoeff = rnorm(K * d + K),
             initDispersion = rep(1, K),
             initmixProp = rep(1 / K, K))
  
  res.lam <- 
    fmrs.tunsel(y = Days, x = Features, 
                delta = Status,
                nComp = ncomp(res.mle), 
                disFamily = 'lnorm',
                initCoeff = c(coefficients(res.mle)),
                initDispersion = dispersion(res.mle),
                initmixProp = mixProp(res.mle),
                penFamily = 'adplasso',
                cutpoint = Tol, LambMin = 0.001, 
                LambMax = 1, nLamb = 1000)
  
  res.var <- 
    fmrs.varsel(y = Days, x = Features, 
                delta = Status,
                nComp = ncomp(res.mle), 
                disFamily = 'lnorm',
                initCoeff = c(coefficients(res.mle)),
                initDispersion = dispersion(res.mle),
                initmixProp = mixProp(res.mle),
                penFamily = 'adplasso',
                lambPen = slot(res.lam, 'lambPen'),
                cutpoint = Tol)
  if (K == 2) {
    Final.Model <- 
      list(res.mle = res.mle, res.lam = res.lam, 
           res.var = res.var)
    Max_BIC <- BIC(res.var)
  } else {
    if (BIC(res.var) <= Max_BIC) {
      Final.Model <- 
        list(res.mle = res.mle, res.lam = res.lam, 
             res.var =  res.var)
      Max_BIC <- BIC(res.var)
    }
  }
}
summary(Final.Model$res.var)

## ---------------------------------------------
lungM <- Final.Model$res.var
cls <- 1 + (lungM@weights[,1] >= 0.5)

PostData <- data.frame( 
  fitted = sapply(1:length(Days), function(x){
    lungM@fitted[x,cls[x]]
  }),
  Component <- as.factor(cls))
library(ggplot2)
ggplot(PostData, aes(x = fitted, group = Component, 
                     color = Component)) + 
  geom_density(linewidth = 1, show.legend = FALSE) +
  stat_density(aes(x = fitted, group = Component, 
                   colour = Component), 
               geom = "line", 
               position = "identity") +
  theme(legend.position = c(0.7, 0.7), 
        legend.title = element_text(size = 7)) + 
  labs(x = "Fitted values")

## ---------------------------------------------
## ------------ session Information ------------
## ---------------------------------------------

sessionInfo()
