L4 data, p = 22 (# composition covariates), m = 6 (# control covariates) and include intercept. Time ranges from PNA = 5 to PNA = 30 and are mapped to standard interval [0,1]. Z (integral of compositional data) is scaled to sd = 1 (column-wise). Select from k_list (degree freedom of cubic splines) in $[4, \dots, 9]$, and lambda.factor is set to 0.1.
library(compReg) #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]])) try2 <- get(load("C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/Real Data/L4/graphcode/leaveoneout/data/[0,1]cv.FuncompCGL.Rda")) # try2$lam[1] # length(try2$lam) # try2$lam[1]/tail(try2$lam, 1) # class(try2$Funcomp.CGL.fit$`4`$W) # dim(try2$Funcomp.CGL.fit$`4`$W) # head(try2$Funcomp.CGL.fit$`4`$Z) # head(try$Funcomp.CGL.fit$`4`$Z) k_list <- 4:9 GIC <-as.list(1:length(k_list)) a <- k_list[1] - 1 n <- length(input$y) p <- 22 # ## with n 34, p = 22, p * dk >= n, p < n # ## GIC1 = GIC3, GIC2 = GIC4, GIC1.1 = GIC3.1, GIC1.2 = GIC3.2, GIC2.1 = GIC4.1, GIC2.2 = GIC4.2, GIC1.3 = GIC3.3, GIC2.3 = GIC4.3 # ## GIC5 = BIC1, GIC6 = BIC2, GIC5.1 = BIC1.1, GIC6.1 = BIC2.1, GIC5.2 = BIC1.2, GIC6.2 = BIC2.2, GIC5.3 = BIC1.2, GIC6.3 = BIC2.3 # ## # ## alpha = switch(GIC_type, # ## "GIC1" = log(max(p*dk,n)) * dk, # ## "GIC2" = log(max(p*dk,n)) * dk, # ## "GIC3" = log(p * dk) * dk, # ## "GIC4" = log(p * dk) * dk, # ## "GIC1.1" = log(max(p*dk, n) * dk), # ## "GIC2.1" = log(max(p*dk, n) * dk), # ## "GIC3.1" = log(p*dk * dk), # ## "GIC4.1" = log(p*dk * dk), # ## "GIC1.2" = log(max(p*dk, n)) * log(dk), # ## "GIC2.2" = log(max(p*dk, n)) * log(dk), # ## "GIC3.2" = log(p*dk) * log(dk), # ## "GIC4.2" = log(p*dk) * log(dk), # ## "GIC1.3" = log(max(p * dk, n)), # ## "GIC2.3" = log(max(p * dk, n)), # ## "GIC3.3" = log(p * dk), # ## "GIC4.3" = log(p * dk), # ## "GIC5" = log(max(p, n)) * dk, # ## "GIC6" = log(max(p, n)) * dk, # ## "GIC7" = log(p) * dk, # ## "GIC8" = log(p) * dk, # ## "GIC5.1" = log(max(p, n) * dk), # ## "GIC6.1" = log(max(p, n) * dk), # ## "GIC5.2" = log(max(p, n)) * log(dk), # ## "GIC6.2" = log(max(p, n)) * log(dk), # ## "GIC7.2" = log(p) * log(dk), # ## "GIC8.2" = log(p) * log(dk), # ## "GIC5.3" = log(max(p, n)), # ## "GIC6.3" = log(max(p, n)), # ## "GIC7.3" = log(p), # ## "GIC8.3" = log(p), # ## "AIC1" = dk * 2, # ## "AIC2" = log(dk) * 2, # ## "BIC1" = dk * log(n), # ## "BIC2" = dk * log(n), # ## "BIC1.1" = log(dk * n), # ## "BIC2.1" = log(dk * n), # ## "BIC1.2" = log(dk) * log(n), # ## "BIC2.2" = log(dk) * log(n), # ## "BIC1.3" = log(n), # ## "BIC2.3" = log(n) # ## ) / n # ## even version of GICs and BICs multiplying log(log(n)) # ## odd version of GICs and BICs multiplying 1 GIC_typelist = c(paste0("GIC", 1:2), paste0("GIC1.", 1:3), paste0("GIC2.", 1:3), paste0("GIC", 7:8), paste0("GIC7.", 2:3), paste0("GIC8.", 2:3), paste0("BIC", 1:2), paste0("BIC1.", 1:3), paste0("BIC2.", 1:3), "AIC1", "AIC2") 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() scaler <- vector() knitr::opts_chunk$set(cache=TRUE)
\newpage
#try$Ftrim$cvm #plot(try, k_list = k_list) par(mfrow=c(1,2)) plot.args_CV = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam), y = try$Ftrim$cvm[1, ], ylim = range(try$Ftrim$cvm), xlab= "lambda Index",#"lambda", #"log(lambda)", ylab="MEAN-Squared Error", main = "10-folds Cross Vadilation", type="n") do.call("plot",plot.args_CV) for(i in 1:length(k_list)) { points(x = seq(length(try$lam)), #log(GIC_m1$lam), y = try$Ftrim$cvm[i, ], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i]) text(length(try$lam), #tail(log(GIC_m1$lam), 1), try$Ftrim$cvm[i, ][length(try$lam)], labels=paste(k_list[i]), cex= 1, pos= 4, col = rainbow(length(k_list))[i]) abline(v = which.min(try$Ftrim$cvm[i, ]), col = rainbow(length(k_list))[i]) } CV_min <- apply(try$Ftrim$cvm, 1, min) k_opt <- which.min(CV_min) + a lam_opt <- which(try$Ftrim$cvm[k_opt - a, ] == CV_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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste0(paste(Non.zero1, collapse = ","), ","), "df=", k_opt), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) text(tail(seq(0, 1, length.out = 20), 1), beta_curve1[dim(beta_curve1)[1], Non.zero1], labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
par(mfrow=c(1,2)) plot.args_CV = list(x = seq(length(try2$lam)), #GIC_m1$lam, #log(GIC_m1$lam), y = try2$Ftrim$cvm[1, ], ylim = range(try2$Ftrim$cvm), xlab= "lambda Index",#"lambda", #"log(lambda)", ylab="MEAN-Squared Error", main = "Leave one out Cross Vadilation", type="n") do.call("plot",plot.args_CV) for(i in 1:length(k_list)) { points(x = seq(length(try2$lam)), #log(GIC_m1$lam), y = try2$Ftrim$cvm[i, ], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i]) text(length(try2$lam), #tail(log(GIC_m1$lam), 1), try2$Ftrim$cvm[i, ][length(try$lam)], labels=paste(k_list[i]), cex= 1, pos= 4, col = rainbow(length(k_list))[i]) abline(v = which.min(try2$Ftrim$cvm[i, ]), col = rainbow(length(k_list))[i]) } CV_min <- apply(try2$Ftrim$cvm, 1, min) k_opt <- which.min(CV_min) + a lam_opt <- which(try2$Ftrim$cvm[k_opt - a, ] == CV_min[k_opt-a]) beta1 <- try2$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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste0(paste(Non.zero1, collapse = ","), ","), "df=", k_opt), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) text(tail(seq(0, 1, length.out = 20), 1), beta_curve1[dim(beta_curve1)[1], Non.zero1], labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{GIC}1(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \times df \times \log\big(\max(p\times df + m + 1,n)\big)+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)$
GIC_type = GIC_typelist[1] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste0(paste(Non.zero1, collapse = ","), ","), "df=", k_opt), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) text(tail(seq(0, 1, length.out = 20), 1), beta_curve1[dim(beta_curve1)[1], Non.zero1], labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{GIC}2(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \times df \times \log\big(\max(p\times df + m + 1,n)\big)\log(\log(n))+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)\log(\log(n))$
GIC_type = GIC_typelist[2] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{GIC}{1.1}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \times df \big)+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)$
GIC_type = GIC_typelist[3] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{GIC}{1.2}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \big) \log(df) + (m+1) \log\big(\max(p\times df + m + 1,n)\big)$
GIC_type = GIC_typelist[4] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{GIC}{1.3}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \big)+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)$
GIC_type = GIC_typelist[5] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{GIC}{2.1}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \times df \big)\log(\log(n))+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)\log(\log(n))$
GIC_type = GIC_typelist[6] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{GIC}{2.2}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \big) \log(df) \log(\log(n)) + (m+1) \log\big(\max(p\times df + m + 1,n)\big) \log(\log(n))$
GIC_type = GIC_typelist[7] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{GIC}{2.3}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \big)\log(\log(n))+ (m+1) \log\big(\max(p\times df + m + 1,n)\big) \log(\log(n))$
GIC_type = GIC_typelist[8] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
$\textrm{BIC}1(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n) \times df + (m+1) \log(n)$
GIC_type = GIC_typelist[15] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{BIC}2(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n) \times df \times \log(\log(n))+ (m+1) \log(n) \log(\log(n))$
GIC_type = GIC_typelist[16] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{BIC}{1.1}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n \times df) + (m+1) \log(n)$
GIC_type = GIC_typelist[17] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{BIC}{1.2}(\lambda) = n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n ) \times \log(df) + (m+1) \log(n)$
GIC_type = GIC_typelist[18] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{BIC}{1.3}(\lambda) = n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n ) + (m+1) \log(n)$
GIC_type = GIC_typelist[19] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{BIC}{2.1}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n \times df)\log(\log(n)) + (m+1) \log(n)\log(\log(n))$
GIC_type = GIC_typelist[20] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{BIC}{2.2}(\lambda) = n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n ) \times \log(df)\log(\log(n)) + (m+1) \log(n)\log(\log(n))$
GIC_type = GIC_typelist[21] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
\newpage
$\textrm{BIC}{2.3}(\lambda) = n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n ) \log(\log(n)) + (m+1) \log(n)\log(\log(n))$
GIC_type = GIC_typelist[23] for(k in k_list){ 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)) 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) # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) dk = k alpha = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) * dk, "GIC2" = log(max(p*dk + 7, n)) * dk, "GIC3" = log(p * dk + 7) * dk, "GIC4" = log(p * dk + 7) * dk, "GIC1.1" = log(max(p*dk + 7, n) * dk), "GIC2.1" = log(max(p*dk + 7, n) * dk), "GIC3.1" = log((p*dk+7) * dk), "GIC4.1" = log((p*dk+7) * dk), "GIC1.2" = log(max(p*dk + 7, n)) * log(dk), "GIC2.2" = log(max(p*dk + 7, n)) * log(dk), "GIC3.2" = log(p*dk + 7) * log(dk), "GIC4.2" = log(p*dk + 7) * log(dk), "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)) * dk, "GIC6" = log(max(p + 7, n)) * dk, "GIC7" = log(p + 7) * dk, "GIC8" = log(p + 7) * dk, "GIC5.1" = log(max(p + 7, n) * dk), "GIC6.1" = log(max(p + 7, n) * dk), "GIC5.2" = log(max(p + 7, n)) * log(dk), "GIC6.2" = log(max(p + 7, n)) * log(dk), "GIC7.2" = log(p + 7) * log(dk), "GIC8.2" = log(p + 7) * log(dk), "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = dk * 2, "AIC2" = log(dk) * 2, "BIC1" = dk * log(n), "BIC2" = dk * log(n), "BIC1.1" = log(dk * n), "BIC2.1" = log(dk * n), "BIC1.2" = log(dk) * log(n), "BIC2.2" = log(dk) * log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n alpha2 = switch(GIC_type, "GIC1" = log(max(p*dk + 7, n)) , "GIC2" = log(max(p*dk + 7, n)) , "GIC3" = log(p * dk + 7) , "GIC4" = log(p * dk + 7) , "GIC1.1" = log(max(p*dk + 7, n)), "GIC2.1" = log(max(p*dk + 7, n)), "GIC3.1" = log((p*dk+7) ), "GIC4.1" = log((p*dk+7) ), "GIC1.2" = log(max(p*dk + 7, n)) , "GIC2.2" = log(max(p*dk + 7, n)) , "GIC3.2" = log(p*dk + 7) , "GIC4.2" = log(p*dk + 7) , "GIC1.3" = log(max(p * dk + 7, n)), "GIC2.3" = log(max(p * dk + 7, n)), "GIC3.3" = log(p * dk + 7), "GIC4.3" = log(p * dk + 7), "GIC5" = log(max(p + 7, n)), "GIC6" = log(max(p + 7, n)), "GIC7" = log(p + 7), "GIC8" = log(p + 7), "GIC5.1" = log(max(p + 7, n)), "GIC6.1" = log(max(p + 7, n) ), "GIC5.2" = log(max(p + 7, n)), "GIC6.2" = log(max(p + 7, n)) , "GIC7.2" = log(p + 7) , "GIC8.2" = log(p + 7) , "GIC5.3" = log(max(p + 7, n)), "GIC6.3" = log(max(p + 7, n)), "GIC7.3" = log(p + 7), "GIC8.3" = log(p + 7), "AIC1" = 2, "AIC2" = 2, "BIC1" = log(n), "BIC2" = log(n), "BIC1.1" = log(n), "BIC2.1" = log(n), "BIC1.2" = log(n), "BIC2.2" = log(n), "BIC1.3" = log(n), "BIC2.3" = log(n) ) / n if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3", "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2", "GIC6.3", "GIC8.3", "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) { c_const = log(log(n)) } else { c_const = 1 } scaler[k-a] <- alpha * c_const part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 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") par(mfrow=c(1,2)) 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) Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))] 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 = "Estimated Curve", xlab = "TIME", #main = "ESTI", ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#, ,type = "l" ) abline(a = 0, b = 0, col = "grey", lwd = 2) title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5) text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1) # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1], # labels = paste(Non.zero1)) par(mfrow=c(1,1))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.