Run_simulation/compCL_test/classo/comp_simulation.R

In jiji6454/Rpac_compReg: Compositional Data Regression

devtools::use_build_ignore("Run_simulation")
library("compReg")

#simulation 1

#parameter set up
n = 50
p = 30
rho = 0.2
sigma = 0.5
gamma = 0.5
add.on = 1:5
beta = c(1, -0.8, 0.6, 0, 0, -1.5, -0.5, 1.2)
beta = c( beta, rep(0, times = p - length(beta)) )
pos_group = c(1:3, 6:8)
neg_group = (1:p)[-pos_group]
intercept = FALSE
# Comp_data = comp_simulation(n = n, p = p, rho = rho, sigma = sigma, gamma  = gamma, add.on = add.on,
#                             beta = beta, intercept = intercept)
#
#
# m1 <- compCL(y = Comp_data$y, Z = Comp_data$X.comp, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
#              pf = rep(1, times = p),
#              lam = NULL, nlam = 100,lambda.factor = ifelse(n < p, 0.05, 0.001),
#              dfmax = 20,
#              #pfmax = min(dfmax * 1.5, p),
#              #u = u,
#              mu_ratio = 1, tol = 1e-10,
#              outer_maxiter = 1e8, outer_eps = 1e-10,
#              inner_maxiter = 1e3, inner_eps = 1e-6)

nfolds = 10
nlam = 100
trim= 0.05
mu_ratio = 1
tol = 1e-10
dfmax = 15

#table set up
rep <- 100
result_minlam <- matrix(NA, nrow = 100, ncol = 6)
colnames(result_minlam) <- c("PE", "L1", "L2", "Linif", "FP", "FN")
result_min1se <- matrix(NA, nrow = 100, ncol = 6)
colnames(result_min1se) <- c("PE", "L1", "L2", "Linif", "FP", "FN")
result_GIC <- matrix(NA, nrow = 100, ncol = 6)
colnames(result_GIC) <- c("PE", "L1", "L2", "Linif", "FP", "FN")
B_minlam <- matrix(NA, nrow = 100, ncol = p + 1)
B_GIC <- matrix(NA, nrow = 100, ncol = p + 1)
B_min1se <- matrix(NA, nrow = 100, ncol = p + 1)

#100 simulation
for(l in 1:rep) {
  set.seed(l )
  cat("replication", l, "\r\n")
  Comp_data = comp_simulation(n = n, p = p, rho = rho, sigma = sigma, gamma  = gamma, add.on = add.on,
                              beta = beta, intercept = intercept)


  #cv.compCL
  foldid <- sample(rep(seq(nfolds), length = n))
  comp_cvm <- cv.compCL(y = Comp_data$y, Z = Comp_data$X.comp, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
                         pf = rep(1, times = p),
                         nlam = nlam, trim = trim,
                         #nfolds = nfolds,
                         foldid = foldid,
                         dfmax = dfmax, 
                         mu_ratio = mu_ratio, tol = tol,
                         outer_maxiter = 1e8, outer_eps = 1e-10, inner_maxiter = 1e3, inner_eps = 1e-6)

  B_minlam[l, ] <- coef(comp_cvm, trim = FALSE, s = "lam.min")$beta
  B_min1se[l, ] <- coef(comp_cvm, trim = FALSE, s = "lam.1se")$beta

  pos_select <- which(abs(B_minlam[l, ]) > 0 )
  FP_select <- pos_select[which( pos_select %in% neg_group )]
  FN_select <- pos_group[which( !(pos_group %in% pos_select) )]
  result_minlam[l, 5] <- length(FP_select)
  result_minlam[l, 6] <- length(FN_select)
  result_minlam[l, 1] <- sum( (Comp_data$y - cbind2(log(Comp_data$X.comp), 1) %*% B_minlam[l, ])^2) / n

  diff <- B_minlam[l, 1:p] - beta
  result_minlam[l, 2] <- sum(abs(diff))
  result_minlam[l, 3] <- sqrt(sum(diff^2))
  result_minlam[l, 4] <- max(abs(diff))
  result_minlam[l, 7] <- sum(diff^2)

  B_min1se[l, ] <- coef(comp_cvm, trim = FALSE, s = "lam.1se")$beta

  pos_select <- which(abs(B_min1se[l, ]) > 0 )
  FP_select <- pos_select[which( pos_select %in% neg_group )]
  FN_select <- pos_group[which( !(pos_group %in% pos_select) )]
  result_min1se[l, 5] <- length(FP_select)
  result_min1se[l, 6] <- length(FN_select)
  result_min1se[l, 1] <- sum( (Comp_data$y - cbind2(log(Comp_data$X.comp), 1) %*% B_min1se[l, ])^2) / n

  diff <- B_min1se[l, 1:p] - beta
  result_min1se[l, 2] <- sum(abs(diff))
  result_min1se[l, 3] <- sqrt(sum(diff^2))
  result_min1se[l, 4] <- max(abs(diff))
  result_min1se[l, 7] <- sum(diff^2)
  comp_GIC <- GIC.compCL(y = Comp_data$y, Z = Comp_data$X.comp, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
                         pf = rep(1, times = p),
                         nlam = nlam, dfmax = dfmax,
                         mu_ratio = mu_ratio, tol = tol,
                         outer_maxiter = 1e8, outer_eps = 1e-12, inner_maxiter = 1e3, inner_eps = 1e-10)

  B_GIC[l, ] <- coef(comp_GIC)$beta

  pos_select <- which(abs(B_GIC[l, ]) > 0 )
  FP_select <- pos_select[which( pos_select %in% neg_group )]
  FN_select <- pos_group[which( !(pos_group %in% pos_select) )]
  result_GIC[l, 5] <- length(FP_select)
  result_GIC[l, 6] <- length(FN_select)
  result_GIC[l, 1] <- sum( (Comp_data$y - cbind2(log(Comp_data$X.comp), 1) %*%  B_GIC[l, ])^2) / n

  diff <-  B_GIC[l, 1:p] - beta
  result_GIC[l, 2] <- sum(abs(diff))
  result_GIC[l, 3] <- sqrt(sum(diff^2))
  result_GIC[l, 4] <- max(abs(diff))
  result_GIC[l, 7] <- sum(diff^2)
}


#formulate tables
table_GIC <- matrix(NA, nrow = 2, ncol = 7)
table_minlam <- matrix(NA, nrow = 2, ncol = 7)
table_min1se <- matrix(NA, nrow = 2, ncol = 7)
colnames(table_GIC) <- c("MSE", "L1", "L2", "Linif", "FP", "FN", "L2^2")
colnames(table_minlam) <- c("MSE", "L1", "L2", "Linif", "FP", "FN", "L2^2")
colnames(table_min1se) <- c("MSE", "L1", "L2", "Linif", "FP", "FN", "L2^2")
rownames(table_GIC) <- c("MEAN", "SD")
rownames(table_minlam) <- c("MEAN", "SD")
rownames(table_min1se) <- c("MEAN", "SD")

table_GIC[1, ] <- apply(result_GIC, 2, mean)
table_GIC[2, ] <- apply(result_GIC, 2, sd) / sqrt(rep)
table_minlam[1, ] <- apply(result_minlam, 2, mean)
table_minlam[2, ] <- apply(result_minlam, 2, sd) / sqrt(rep)
table_min1se[1, ] <- apply(result_min1se, 2, mean)
table_min1se[2, ] <- apply(result_min1se, 2, sd) / sqrt(rep)



# result_GIC2 <- matrix(NA, nrow = 100, ncol = 6)
# colnames(result_GIC2) <- c("PE", "L1", "L2", "Linif", "FP", "FN")
# B_GIC2 <- matrix(NA, nrow = 100, ncol = p + 1)


# for(l in 1:rep) {
#   set.seed(l+100)
#   Comp_data = comp_simulation(n = n, p = p, rho = rho, sigma = sigma, gamma  = gamma, add.on = add.on,
#                               beta = beta, intercept = intercept)
#
#
#   #cv.compCL
#   foldid <- sample(rep(seq(nfolds), length = n))
#
#   comp_GIC <- GIC.compCL(y = Comp_data$y, Z = Comp_data$X.comp, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
#                          pf = rep(1, times = p),
#                          nlam = nlam, dfmax = dfmax,
#                          mu_ratio = mu_ratio, tol = tol,
#                          outer_maxiter = 1e8, outer_eps = 1e-10, inner_maxiter = 1e3, inner_eps = 1e-6)
#
#   B_GIC2[l, ] <- coef(comp_GIC)$beta
#
#   pos_select <- which(abs(B_GIC[l, ]) > 0 )
#   FP_select <- pos_select[which( pos_select %in% neg_group )]
#   FN_select <- pos_group[which( !(pos_group %in% pos_select) )]
#   result_GIC2[l, 5] <- length(FP_select)
#   result_GIC2[l, 6] <- length(FN_select)
#   result_GIC2[l, 1] <- sum( (Comp_data$y - cbind2(log(Comp_data$X.comp), 1) %*%  B_GIC2[l, ])^2) / n
#
#   diff <-  B_GIC[l, 1:p] - beta
#   result_GIC2[l, 2] <- sum(abs(diff))
#   result_GIC2[l, 3] <- sqrt(sum(diff^2))
#   result_GIC2[l, 4] <- max(abs(diff))
# }


PE <- matrix(NA, nrow = 100, ncol = 3)
for(l in 1:rep) {
  set.seed(l + 200)
  Comp_data = comp_simulation(n = n, p = p, rho = rho, sigma = sigma, gamma  = gamma, add.on = add.on,
                              beta = beta, intercept = intercept)

  PE[l, 1] <- sum((Comp_data$y - log(Comp_data$X.com) %*% B_minlam[l, -(p+1)])^2)/n
  PE[l, 2] <- sum((Comp_data$y - log(Comp_data$X.com) %*% B_min1se[l, -(p+1)])^2)/n
  PE[l, 3] <- sum((Comp_data$y - log(Comp_data$X.com) %*% B_GIC[l, -(p+1)])^2)/n
}


table_GIC <- cbind(table_GIC, c(MEAN = mean(PE[, 3]), SD = sd(PE[, 3])/ sqrt(rep) ) )
table_minlam <- cbind(table_minlam, c(MEAN = mean(PE[, 1]), SD = sd(PE[, 1])/ sqrt(rep) ))
table_min1se <- cbind(table_min1se, c(MEAN = mean(PE[, 2]), SD = sd(PE[, 2])/ sqrt(rep) ))

colnames(table_GIC)[8] <- "PE"
colnames(table_minlam)[8] <- "PE"
colnames(table_min1se)[8] <- "PE"

result <- list(minlam = table_minlam,
               min1se = table_min1se,
               GIC = table_GIC)
save(result, file = "C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/program/compReg/Run_simulation/compCL_test/classo/Lin_setting1-1_compReg_muratio1.Rda")
mean(apply(B_min1se[, -(p + 1)], 1, function(beta, x) sum((x - beta)^2), beta = beta))

a <- vector()
b <- vector()
for(l in 1:100) {
  pos_select <- which(abs(B_GIC[l, ]) > 0.01 )
  a[l] <- length(pos_select[which( pos_select %in% neg_group )])
  b[l] <- length(pos_group[which( !(pos_group %in% pos_select) )])
}