Rpac_compReg: Compositional Data Regression

#try$Funcomp.CGL.fit[[1]]
input <- get(load("C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/Real Data/L4/graphcode/leaveoneout/data/input_K=4.Rda"))
apply(input$X / 25, 2, sd)
try <- get(load("C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/Real Data/L4/graphcode/leaveoneout/data/[0,1]10folds_Weighted_cv.FuncompCGL-0.1.Rda"))
names(try)
as.numeric(names(try[[1]]))
k_list <- 4:9
GIC <-as.list(1:length(k_list))
a <- k_list[1] - 1
n <- length(input$y)
p <- 22
MSE <- matrix(NA, nrow = length(k_list), ncol = length(try$lam))
N_zero <- matrix(NA, nrow = length(k_list), ncol = length(try$lam))
GIC_min <- vector()
# cat("K: ")
for(k in k_list){
  cat("K=", k, "; ")
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*%  try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  
  
  #part1 <- log(apply(pre, 2, function(x, y) sum((x- y)^2)/ (length(y) - 7), y = input$y))
  ##part1 <- log(apply(pre, 2, function(x, y) sum((x- y)^2)/ (length(y) - 7), y = input$y)) #- log (summary(lm(input$y~input$Z-1))$sigma^2) #log(apply(pre, 2, function(x, y) sum((x- y)^2)/length(y), y = input$y))
  ## part1 <- apply(pre, 2, function(x, y) sum((x- y)^2), y = input$y) / summary(lm(input$y~input$Z-1))$sigma^2 / n
  
  ##scaler <- log(log(n)) * log(max(p*k, n)) / n #log(n) / n * log(max(p, n)) #log(log(n)) / n * log(max(p, n))
  #scaler <- log(log(n)) * log(max(p * k, n)) / n
  
  #scaler <- log(log(n)) * log(max(p * k, n)) / n
  #scaler <- log(max(p * k, n)) / n
  scaler <- log(n) / n
  #scaler <- log(log(n)) * log(n) / n
  #scaler <- log(max(p*n)) / n
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2, function(x, p, k) {
    A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
    a = apply(A, 1, function(x) sum(abs(x)))
    return(length(which(a > 0)))
  },p = p, k = k)
  
  ##part2 <- ifelse(N_zero >=2, scaler*(N_zero*k - k + 7) , scaler*7 )
  #part2 <- scaler* ifelse(N_zeroN_zero[k-a, ] >= 2, 
  #                        N_zeroN_zero[k-a, ] - 1 + 7, 7)
  part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
                           (N_zero[k-a, ] - 1) * k + 7, 7)
  plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])
  
  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])
  
}
cat("\r\n")
plot.args_MSE = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = log(MSE)[1, ],
                     ylim = range(log(MSE)),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="log(MSE)",
                     type="n",
                     main = "MSE")
do.call("plot",plot.args_MSE)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = log(MSE[i, ]), col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       log(MSE[i, length(try$lam)]), labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
}


plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "coef", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        #type = "l"
)

# beta_W0.1 <- beta1
# beta_UW0.1 <- beta1
# # beta_choose <- list(k = k_se, beta = beta1)
# # save(beta_choose, file = "UWeighted.GIC_beta.Rda")
# 
# cvm <- get(load("C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/Real Data/L4/graphcode/leaveoneout/data/[0,1]cv.FuncompCGL.Rda"))
# beta_W0.1_cvm <- coef(cvm, s = "lam.min")
# beta_W0.1_cvm <- beta_W0.1_cvm[-length(beta_W0.1_cvm)]
# tail(beta_W0.1_cvm, 7)
# tail(beta_W0.1, 7)