initialize.R

initialize = function(){
  

################################# GIBBS SAMPLING  ###################################################

Time <- cbind(time, censoring) 

K <<-  as.integer(N/5)

## HYPER PRIORS
## Hyper parameters of the DP
shape.alpha <- 2
rate.alpha <- 1
## Hyperparameters for the GMM
beta  = D+1
ro = 0.5


source('rchinese.R')
alpha  = rgamma(1, shape = shape.alpha, rate = rate.alpha )
c <-  rchinese(N,alpha)
f <- table(factor(c, levels = 1:max(c)))

## Empirical Bayes Estimate of the Hyperparameters
epsilon = as.vector(apply(Y,2,mean))
W = diag(diag(cov(Y)))


## Initialization of the parameters for Gaussian Mixture
mu = matrix(data = NA, nrow = K, ncol = D)
S = array(data = NA, dim =c(K,D,D))


#Sparsity controlling hyperparameter of the BAYESIAN LASSO MODEL
r =1
si = 1.78

lambda2 <- numeric(K)
tau2 = matrix(data = NA, nrow = K, ncol = D)
betahat = matrix(data = NA, nrow = K, ncol = D)
sigma2 <- rep(NA, K)
beta0 <- rep(NA, K)
That <-  numeric(N)

## Fitting a linear model to the whole model
Ysc <- scale(Y[1:N,1:D], center = TRUE, scale =TRUE)
lm.data <- lm(time ~ Ysc)
sig2.dat <-  var(lm.data$residuals)


## Set Some Initial Values for the Cluster Parameters

source('priorPARAMETERS.R')
disclass <- table(factor(c, levels = 1:K))
activeclass <- which(disclass!=0)
for ( j in 1:length(activeclass)){
  
  priorone <- priordraw(beta, W, epsilon, ro, r, si, N, D, sig2.dat)  
  mu[activeclass[j],] <- (priorone$mu) 
  S[activeclass[j],1:D,1:D]  <- priorone$Sigma  
  beta0[activeclass[j]] <- priorone$beta0 
  sigma2[activeclass[j]] <- priorone$sigma2
  betahat[activeclass[j],1:D] <- priorone$betahat 
  lambda2[activeclass[j]] <- priorone$lambda2 
  tau2[activeclass[j], 1:D] <- priorone$tau2
}

# The Time has to be initialized
source('posteriorCensoredTime.R')
ti <- updatetime(c, Y, Time,That, beta0, betahat, sigma2)
That <- ti$time


################# K-Means BLASSO INITIALIZATION ############################################
G <- F

k.data <- kmeans(Y,F,nstart =10)

c <- k.data$cluster

c.kmeans <<- c
#### Under special cases
###### c <- c.true


prior.numclust <- table(factor(c, levels = 1:K))
prior.activeclass<- which(prior.numclust!=0)

### The means are set using the k-means
for ( i in 1:length(prior.activeclass)){
  mu[prior.activeclass[i],1:D] <-  k.data$centers[i,1:D] 
  
  S[prior.activeclass[i],1:D,1:D] <-  priordraw(beta, W, epsilon, ro, r, si,N,D, sig2.dat)$Sigma
  
  lclust <- which(c == prior.activeclass[i])
  
  reg.blas <- 0
  
  sum <- c(0)
  
  coeff <- 0
  
  Ytemp <-  matrix(NA, nrow = length(lclust), ncol = D)
  
  Ytemp <- scale(Y[lclust,1:D], center = TRUE, scale = TRUE)
  
  
  ### Part where I use the MONOMVN PACKAGE
  
  Ttemp <- as.vector(That[lclust])
  
  ntemp <- length(lclust)
  
  reg.blas <- blasso(Ytemp, Ttemp, T = 300,thin =  50, RJ = TRUE, mprior = 0.0 ,normalize = TRUE, verb = 0)
  
  sum <- summary(reg.blas, burnin= 100)
  
  ## Selecting those features which are relevant
  
  coeff <- unlist(lapply(strsplit(sum$coef[3,], split = ":"), function(x) as.numeric(unlist(x)[2])))
  
  beta0[prior.activeclass[1]] <- coeff[1]
  
  indexplusone <- D+1
  
  ind <- 2:indexplusone
  
  betahat[prior.activeclass[i], ] <- coeff[ind]
  
  ta <- unlist(lapply(strsplit(sum$tau2i[3,], split = ":"), function(x) as.numeric(unlist(x)[2])))
  
  tau2[prior.activeclass[i],] <- ta
  
  sigma2[prior.activeclass[i]] <- sum$s2[3]
  
  lambda2[prior.activeclass[i]] <- sum$lambda2[3]
  
}

## Deleting those values which are no longer relevant
g <- table(factor(c, levels = 1:K))
inactive <- which(g==0)

for ( i in 1:length(inactive)){
  mu[inactive[i],1:D]  <- NA 
  S[inactive[i],1:D,1:D]  <- NA  
  beta0[inactive[i]] <- NA 
  sigma2[inactive[i]] <- NA
  betahat[inactive[i],1:D] <- NA 
  lambda2[inactive[i]] <- NA
  tau2[inactive[i], 1:D] <- NA
}





assign("Time", Time, envir = .GlobalEnv)
assign("r", r, envir = .GlobalEnv)
assign("si", si, envir = .GlobalEnv)
assign("shape.alpha", shape.alpha, envir = .GlobalEnv)
assign("rate.alpha", shape.alpha, envir = .GlobalEnv)
assign("alpha", alpha, envir = .GlobalEnv)
assign("beta", beta, envir = .GlobalEnv)
assign("ro", ro, envir = .GlobalEnv)
assign("That", That, envir = .GlobalEnv)
assign("c", c, envir = .GlobalEnv)
assign("epsilon", epsilon, envir = .GlobalEnv)
assign("W", W, envir = .GlobalEnv)
assign("mu", mu, envir = .GlobalEnv)
assign("S", S, envir = .GlobalEnv)
assign("beta0", beta0, envir = .GlobalEnv)
assign("betahat", betahat, envir = .GlobalEnv)
assign("sigma2", sigma2, envir = .GlobalEnv)
assign("lambda2", lambda2, envir = .GlobalEnv)
assign("tau2", tau2, envir = .GlobalEnv)
assign("lambda2", lambda2, envir = .GlobalEnv)
assign("sig2.dat", sig2.dat, envir = .GlobalEnv)

}