GIC.FuncompCGL: GIC cirterion selection for FuncompCGL

In jiji6454/Rpac_compReg: Compositional Data Regression

Description

Calculate GIC for compCL, return value of lam.

Usage

GIC.FuncompCGL(y, X, Zc = NULL, ref = NULL, lam = NULL, nlam = 100,
  W = rep(1, times = p - length(ref)), k = 4:10, outer_maxiter = 1e+06,
  ...)

Arguments

y	a vector of response variable.
X	a data frame of longitudinal compositinal predictors with number p, subject ID and time variable. Order of subject ID should be the same as that of y. If df k is a scalar, X could be the matrix with dimension n*(k*p - length(ref)) after integral.
Zc	A design matrix for control variables, could be missing. Default is NULL. No penalty is imposed.
ref	reference variable. If ref is set to a scalar between [1,p], log-contract method is applied with the variable ref as baseline. If ref = NULL (default value), constrained group lasso method is applied
lam	a user supplied lambda sequence. Typically, by leaving this option unspecified users can have the program compute its own lam sequence based on nlam and lambda.factor If lam is provided but a scaler, lam sequence is also created starting from lam. Supplying a value of lambda overrides this. It is better to supply a decreasing sequence of lambda values, if not, the program will sort user-defined lambda sequence in decreasing order automatically.
nlam	the length of lam sequence - default is 100.
W	a vector in length of p (the total number of groups), matrix with dimension p1*p1 or character specifying function used to calculate inverted weight matrix for each group. If vector, works as penalty factor. Separate penalty weights can be applied to each group of beta'ss. to allow differential shrinkage. Can be 0 for some groups, which implies no shrinkage, and results in that group always being included in the model. If matrix, a block diagonal matrix. Diagonal elements are inverted weights matrics for each group. if character, user should provide the function for inverted weights matrics. Default value is rep(1, times = p).
k	a vector or scalar consists of df for basis - default is 4:10.
outer_maxiter	maximun munber of loops allowed for Augmented Lanrange method.
...	other arguments that could be passed to FuncompCL.

Value

an object of class GIC.FuncompCGL is returned.

Funcomp.CGL.fit	a list, length of k, of fitted FuncompCGL object for the full data. objects with S3 calss FuncompCGL
lam	the values of lam used in the fits
MSE	matrix of mean squared error with size k by nlam (the length of actually used lambda sequence, migth pre-stop by dfmax or pfmax). MSE is equivalent to likelihood under normal error model. Could be edited for other linkage.
Nzero	a k by nlam matrix for Nzero group cut-off by cut_off and lower_tri

Examples

df_beta = 5
p = 30 #30
beta_C_true = matrix(0, nrow = p, ncol = df_beta)
beta_C_true[3, ] <- c(-1, 0, 0, 0, -0.5) #c(-0.5, 0, 0, 0, -0.5)
beta_C_true[1, ] <- c(1, 0, 1 , 0, -0.5) #c(0.5, 0, 1 , 0, -0.5)
beta_C_true[2, ] <- c(0, 0,  -1,  0,  1)

nfolds = 10
k_list <- c(4,5,6)
n_train = 100 #100
n_test = 500

Data <- Model2(n = n_train, p = p, m = 0, intercept = TRUE,
               SNR = 3, sigma = 3,
               rho_X = 0, rho_W = 0,
               Corr_X = "CorrCS", Corr_W = "CorrAR",
               df_W = 5, df_beta = df_beta,
               ns = 20, obs_spar = 1, theta.add = FALSE, #c(0,0,0),
               beta_C = as.vector(t(beta_C_true)))
y <- drop(Data$data$y)
n <- length(y)
X <- Data$data$Comp
Zc <- Data$data$Zc
intercept <- Data$data$intercept
m <- ifelse(is.null(Zc), 0, dim(Zc)[2]) #+ as.integer(intercept)
m1 <- m + as.integer(intercept)
sseq <- Data$basis.info[,1]
beta_C.true <- matrix(Data$beta[1:(p*(df_beta))],
                      nrow = p, ncol = df_beta, byrow = TRUE)
beta_curve.true <- Data$basis.info[,-1] %*% t(beta_C.true)
Non_zero.true <- (1:p)[apply(beta_C.true, 1, function(x) max(abs(x)) > 0)]
foldid <- sample(rep(seq(nfolds), length = n))

arg_list <- as.list(Data$call)[-1]
arg_list$n <- n_test
Test <- do.call(Model2, arg_list)
y_test <- drop(Test$data$y)

GIC_arg <- list(basis_fun = "bs", degree = 3, sseq = Data$basis.info[, 1])

# y = y
# X = X
# Zc = Zc
# intercept = intercept
# W = rep(1, p)
# k = k_list
# nfolds = 10
# trim = 0
# tol = 0
# inner_eps = 1e-6
# inner_maxiter = 1E3
# dfmax = 20
# lambda.factor = 1e-20
# mu_ratio = 1
# outer_eps = 1e-6
# keep = TRUE
# Trange = c(0,1)



GIC_m1 <- GIC.FuncompCGL( y = y, X = X, Zc = Zc, ref = NULL,
                          inner_eps = 1e-8, outer_eps = 1e-8, tol = 1e-8,
                          k = k_list)

temp <- get.GIC(p = p, df_list = k_list, lower_tri = 0.01,
                GIC_obj = GIC_m1, GIC_arg = GIC_arg,
                cut_type = "Strict", GIC_type = "GIC1",
                method_type = "cgl", refit = FALSE)
GIC_curve <- temp$GIC_curve
k_opt <- temp$k_opt
beta_GIC <- temp$beta

plot.args = list(x = seq(length(GIC_m1$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                 y = GIC_curve[1, ],
                 ylim = range(GIC_curve),
                 xlab= "lambda Index",#"lambda", #"log(lambda)",
                 ylab="GIC",
                 type="n")
#do.call("plot",plot.args)
# for(i in 1:length(k_list)) {
#
#   points(x = seq(length(GIC_m1$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
#          y = GIC_curve[i, ], col = rainbow(length(k_list))[i])
#   text(length(GIC_m1$lam), #tail(log(GIC_m1$lam), 1),
#        GIC_curve[i, length(GIC_m1$lam)], labels=paste(k_list[i]),
#        cex= 1, pos= 4, col = rainbow(length(k_list))[i])
# }
# axis(3, at = pretty(seq(length(GIC_m1$lam))), labels = rev(pretty(GIC_m1$lam)))
# loc  = which(GIC_curve == min(GIC_curve), arr.ind = TRUE)




beta_C <- matrix(beta_GIC[1:(p*k_opt)], byrow = TRUE, nrow = p)
cat("colSums:", colSums(beta_C) , "\r\n")
#Non.zero <- which(abs(beta_C[,1]) > 0)
Non.zero <- apply(beta_C, 1, function(x) ifelse(max(abs(x)) >0, TRUE, FALSE))
Non.zero <- (1:p)[Non.zero]
cat("None zero groups:", Non.zero)
#vet(beta, p = p, k = k_opt)

par(mfrow=c(1,4))
do.call("plot",plot.args)
for(i in 1:length(k_list)) {
  points(x = seq(length(GIC_m1$lam)), #log(GIC_m1$lam),
         y = GIC_curve[i, ], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(GIC_m1$lam), #tail(log(GIC_m1$lam), 1),
       GIC_curve[i, length(GIC_m1$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
}
#axis(3, at = pretty(seq(length(GIC_m1$lam))), labels = rev(pretty(GIC_m1$lam)))

matplot(sseq, beta_curve.true,
        ylab = "coeffcients curve", xlab = "TIME", #main = "TRUE",
        ylim = range(Data$beta[1:(p*df_beta)]),
        type = "l")
abline(a = 0, b = 0, col = "grey", lwd = 2)
title("TRUE", line = 0.5)
text(0, beta_curve.true[1, Non_zero.true], labels = paste(Non_zero.true))

B <- splines::bs(Data$basis.info[,1], df = k_opt, intercept = TRUE)
beta_curve <- B %*% t(beta_C)
matplot(sseq, beta_curve,
        ylab = "coef", xlab = "TIME", #main = "ESTI",
        ylim = range(Data$beta[1:(p*df_beta)])#,
        #type = "l"
)
abline(a = 0, b = 0, col = "grey", lwd = 2)
title("Estimate", line = 0.5)
text(0, beta_curve[1, Non.zero], labels = paste(Non.zero))
text(tail(sseq, 1), beta_curve[dim(beta_curve)[1], Non.zero], labels = paste(Non.zero))
plot(apply(abs(beta_C),1,sum))
text(seq(length(GIC_m1$lam))[which(apply(abs(beta_C),1,sum) > 0)], #tail(log(GIC_m1$lam), 1),
     apply(abs(beta_C),1,sum)[which(apply(abs(beta_C),1,sum) > 0)],
     labels=paste(seq(length(GIC_m1$lam))[which(apply(abs(beta_C),1,sum) > 0)]),
     cex= 1, pos= 4)

title(paste0("k=", k_opt), line = 0.5)
title(paste0("Method cgl"), outer=TRUE, line = -2)
par(mfrow=c(1,1))

##set a cutoff when you compute nonzeros
Non.zero <- apply(beta_C, 1, function(x)
  ifelse(sqrt(sum(x^2)) > sqrt(sum(beta_C^2))/100, TRUE, FALSE))
Non.zero <- (1:p)[Non.zero]
Non.zero

cgl_GIC <- list()
MSE <- temp$MSE[temp$k_opt - k_list[1] + 1, temp$lam_loc]
R_sqr <- 1 - MSE * length(y) / crossprod(y -  mean(y))

obj <- FuncompCGL(y = y_test, X = Test$data$Comp, k = k_opt, nlam = 1, outer_maxiter = 0)
X_test <- cbind2(cbind(obj$Z, Test$data$Zc), 1)
PE <- sum((y_test - X_test %*% beta_GIC)^2) / length(y_test)
cgl_GIC$pred_error <- c(MSE = MSE, PE = PE, Rsqr_train = R_sqr)

cgl_GIC$Non.zero_cut <- Non.zero
cgl_GIC <- c(cgl_GIC,
             ERROR_fun(beta_fit = beta_GIC, beta_true = Data$beta,
                       basis_fit = B, basis_true = Data$basis.info[,-1],
                       sseq = Data$basis.info[, 1],
                       m = m, p = p, Nzero_group = length(Non_zero.true), tol = 0),
             k = k_opt)
cgl_GIC$coef <- list(beta_C = beta_C, beta_c = tail(beta_GIC, m1))

## Not run: 

GIC_m2 <- GIC.FuncompCGL( y = y, X = X, Zc = Zc, ref = NULL,
                          outer_eps = 1e-8, mu_ratio = 0, tol = 1e-8,
                          k = k_list)
temp <- get.GIC(p = p, df_list = k_list, lower_tri = 0.01,
                GIC_obj = GIC_m1, GIC_arg = GIC_arg,
                cut_type = "Strict", GIC_type = "GIC1",
                method_type = "naive", refit = FALSE)
GIC_curve <- temp$GIC_curve
k_opt <- temp$k_opt
beta_GIC <- temp$beta

plot.args = list(x = seq(length(GIC_m2$lam)), #GIC_m2$lam, #log(GIC_m2$lam),
                 y = GIC_curve[1, ],
                 ylim = range(GIC_curve),
                 xlab= "lambda Index",#"lambda", #"log(lambda)",
                 ylab="GIC",
                 type="n")
# do.call("plot",plot.args)
#
# for(i in 1:length(k_list)) {
#
#   points(x = seq(length(GIC_m2$lam)), #GIC_m2$lam, #log(GIC_m2$lam),
#          y = GIC_curve[i, ], col = rainbow(length(k_list))[i])
#   text(length(GIC_m2$lam), #tail(log(GIC_m2$lam), 1),
#        GIC_curve[i, length(GIC_m2$lam)], labels=paste(k_list[i]),
#        cex= 1, pos= 4, col = rainbow(length(k_list))[i])
# }
# axis(3, at = pretty(seq(length(GIC_m2$lam))), labels = rev(pretty(GIC_m2$lam)))

beta_C <- matrix(beta_GIC[1:(p*k_opt)], byrow = TRUE, nrow = p)
cat("colSums:", colSums(beta_C), "\r\n")
#Non.zero <- which(abs(beta_C[,1]) > 0)
Non.zero <- apply(beta_C, 1, function(x) ifelse(max(abs(x)) >0, TRUE, FALSE))
Non.zero <- (1:p)[Non.zero]
cat("None zero groups:", Non.zero)
#vet(beta, p = p, k = k_opt)

par(mfrow=c(1,4))

do.call("plot",plot.args)
for(i in 1:length(k_list)) {
  points(x = seq(length(GIC_m2$lam)), #GIC_m2$lam, #log(GIC_m2$lam),
         y = GIC_curve[i, ], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(GIC_m2$lam), #tail(log(GIC_m2$lam), 1),
       GIC_curve[i, length(GIC_m2$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
}
#axis(3, at = pretty(seq(length(GIC_m2$lam))), labels = rev(pretty(GIC_m2$lam)))

matplot(sseq, beta_curve.true,
        ylab = "coeffcients curve", xlab = "TIME", #main = "TRUE",
        ylim = range(Data$beta[1:(p*df_beta)]),
        type = "l")
abline(a = 0, b = 0, col = "grey", lwd = 2)
title("TRUE", line = 0.5)
text(0, beta_curve.true[1, Non_zero.true], labels = paste(Non_zero.true))

B <- splines::bs(Data$basis.info[,1], df = k_opt, intercept = TRUE)
beta_curve <- B %*% t(beta_C)
matplot(sseq, beta_curve,
        ylab = "coef", xlab = "TIME", #main = "ESTI",
        ylim = range(Data$beta[1:(p*df_beta)])#,
        #type = "l"
)
abline(a = 0, b = 0, col = "grey", lwd = 2)
title("Estimate", line = 0.5)
text(0, beta_curve[1, Non.zero], labels = paste(Non.zero))
text(tail(sseq, 1), beta_curve[dim(beta_curve)[1], Non.zero], labels = paste(Non.zero))


plot(apply(abs(beta_C),1,sum))
text(seq(length(GIC_m2$lam))[which(apply(abs(beta_C),1,sum) > 0)], #tail(log(GIC_m2$lam), 1),
     apply(abs(beta_C),1,sum)[which(apply(abs(beta_C),1,sum) > 0)],
     labels=paste(seq(length(GIC_m2$lam))[which(apply(abs(beta_C),1,sum) > 0)]),
     cex= 1, pos= 4)

title(paste0("k=", k_opt), line = 0.5)
title(paste0("Method naive"), outer=TRUE, line = -2)

par(mfrow=c(1,1))


##set a cutoff when you compute nonzeros
Non.zero <- apply(beta_C, 1, function(x)
  ifelse(sqrt(sum(x^2)) > sqrt(sum(beta_C^2))/100, TRUE, FALSE))
Non.zero <- (1:p)[Non.zero]
Non.zero

naive_GIC <- list()
MSE <- temp$MSE[temp$k_opt - k_list[1] + 1, temp$lam_loc]
R_sqr <- 1 - MSE * length(y) / crossprod(y -  mean(y))
obj <- FuncompCGL(y = y_test, X = Test$data$Comp, k = k_opt, nlam = 1, outer_maxiter = 0)
X_test <- cbind2(cbind(obj$Z, Test$data$Zc), 1)
PE <- sum((y_test - X_test %*% beta_GIC)^2) / length(y_test)
naive_GIC$pred_error <- c(MSE = MSE, PE = PE, Rsqr_train = R_sqr)
naive_GIC$Non.zero_cut <- Non.zero
naive_GIC <- c(naive_GIC,
               ERROR_fun(beta_fit = beta_GIC, beta_true = Data$beta,
                         basis_fit = B, basis_true = Data$basis.info[,-1],
                         sseq = Data$basis.info[, 1],
                         m = m, p = p, Nzero_group = length(Non_zero.true), tol = 0),
               k = k_opt)
naive_GIC$coef <- list(beta_C = beta_C, beta_c = tail(beta_GIC, m1))



GIC_m3 <- GIC.FuncompCGL(y = y, X = X, Zc = Zc, ref = sample(4:p, 1), #sample(1:p, 1),
                         outer_eps = 1e-8, mu_ratio = 0, tol = 1e-8,
                         k = k_list)
temp <- get.GIC(p = p, df_list = k_list, lower_tri = 0.01,
                GIC_obj = GIC_m1, GIC_arg = GIC_arg,
                cut_type = "Strict", GIC_type = "GIC1",
                method_type = "base", refit = FALSE)
GIC_curve <- temp$GIC_curve
k_opt <- temp$k_opt
beta_GIC <- temp$beta
ref = GIC_m3$Funcomp.CGL.fit[[1]]$ref
beta_GIC <-  c(beta_GIC1[ifelse(ref==1, 0, 1):((ref-1)*k_opt)],
               -colSums(matrix(beta_GIC1[1:((p-1)*k_opt)], byrow = TRUE, ncol = k_opt)),
               beta_GIC1[((ref-1)*k_opt+1):(length(beta_GIC1))])

plot.args = list(x = seq(length(GIC_m3$lam)), #GIC_m3$lam, #log(GIC_m3$lam),
                 y = GIC_curve[1, ],
                 ylim = range(GIC_curve),
                 xlab= "lambda index",#"lambda", #"log(lambda)",
                 ylab="GIC",
                 type="n")
# do.call("plot",plot.args)
#
# for(i in 1:length(k_list)) {
#
#   points(x = seq(length(GIC_m3$lam)), #GIC_m3$lam, #log(GIC_m3$lam),
#          y = GIC_curve[i, ], col = rainbow(length(k_list))[i])
#   text(length(GIC_m3$lam), #tail(log(GIC_m3$lam), 1),
#        GIC_curve[i, length(GIC_m3$lam)], labels=paste(k_list[i]),
#        cex= 1, pos= 4, col = rainbow(length(k_list))[i])
#
# }
# axis(3, at = seq(1, length(GIC_m3$lam), length.out = 10),
#      labels = rev(GIC_m3$lam[seq(1, length(GIC_m3$lam), length.out = 10)]) )
# axis(3, at = pretty(seq(length(GIC_m3$lam))),
#         labels = rev(pretty(GIC_m3$lam, n = length(pretty(seq(length(GIC_m3$lam))))))
#      )



beta_C <- matrix(beta_GIC[1:(p*k_opt)], byrow = TRUE, nrow = p)
cat("colSums:", colSums(beta_C), "\r\n")
#Non.zero <- which(abs(beta_C[,1]) > 0)
Non.zero <- apply(beta_C, 1, function(x) ifelse(max(abs(x)) >0, TRUE, FALSE))
Non.zero <- (1:p)[Non.zero]
cat("None zero groups:", Non.zero)
#vet(beta, p = p, k = k_opt)

par(mfrow=c(1,4))
do.call("plot",plot.args)
for(i in 1:length(k_list)) {
  points(x = seq(length(GIC_m3$lam)), #log(GIC_m3$lam),
         y = GIC_curve[i, ], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(GIC_m3$lam), #tail(log(GIC_m3$lam), 1),
       GIC_curve[i, length(GIC_m3$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
}
#axis(3, at = pretty(seq(length(GIC_m3$lam))), labels = rev(pretty(GIC_m3$lam)))
matplot(sseq, beta_curve.true,
        ylab = "coeffcients curve", xlab = "TIME", #main = "TRUE",
        ylim = range(Data$beta[1:(p*df_beta)]),
        type = "l")
abline(a = 0, b = 0, col = "grey", lwd = 2)
title("TRUE", line = 0.5)
text(0, beta_curve.true[1, Non_zero.true], labels = paste(Non_zero.true))

B <- splines::bs(Data$basis.info[,1], df = k_opt, intercept = TRUE)
beta_curve <- B %*% t(beta_C)
matplot(sseq, beta_curve,
        ylab = "coef", xlab = "TIME", #main = "ESTI",
        ylim = range(Data$beta[1:(p*df_beta)])
        #type = "l"
)
abline(a = 0, b = 0, col = "grey", lwd = 2)
title("Estimate", line = 0.5)
text(0, beta_curve[1, Non.zero], labels = paste(Non.zero))
text(tail(sseq, 1), beta_curve[dim(beta_curve)[1], Non.zero], labels = paste(Non.zero))
plot(apply(abs(beta_C),1,sum))
text(seq(length(GIC_m3$lam))[which(apply(abs(beta_C),1,sum) > 0)], #tail(log(GIC_m3$lam), 1),
     apply(abs(beta_C),1,sum)[which(apply(abs(beta_C),1,sum) > 0)],
     labels=paste(seq(length(GIC_m3$lam))[which(apply(abs(beta_C),1,sum) > 0)]),
     cex= 1, pos= 4)
title(paste0("k=", k_opt), line = 0.5)
title(paste0("Method base",  ", ref=", ref), outer=TRUE, line = -2)
par(mfrow=c(1,1))
##set a cutoff when you compute nonzeros
Non.zero <- apply(beta_C, 1, function(x)
  ifelse(sqrt(sum(x^2)) > sqrt(sum(beta_C^2))/100, TRUE, FALSE))
Non.zero <- (1:p)[Non.zero]
Non.zero


base_GIC <- list()
MSE <- temp$MSE[temp$k_opt - k_list[1] + 1, temp$lam_loc]
R_sqr <- 1 - MSE * length(y) / crossprod(y -  mean(y))
obj <- FuncompCGL(y = y_test, X = Test$data$Comp, k = k_opt, nlam = 1, outer_maxiter = 0)
X_test <- cbind2(cbind(obj$Z, Test$data$Zc), 1)
PE <- sum((y_test - X_test %*% beta_GIC)^2) / length(y_test)
base_GIC$pred_error <- c(MSE = MSE, PE = PE, Rsqr_train = R_sqr)
base_GIC$Non.zero_cut <- Non.zero
base_GIC <- c(base_GIC,
              ERROR_fun(beta_fit = beta_GIC, beta_true = Data$beta,
                        basis_fit = B, basis_true = Data$basis.info[,-1],
                        sseq = Data$basis.info[, 1],
                        m = m, p = p, Nzero_group = length(Non_zero.true), tol = 0),
              k = k_opt)
base_GIC$coef <- list(beta_C = beta_C, beta_c = tail(beta_GIC, m1))


## End(Not run)

jiji6454/Rpac_compReg documentation built on May 31, 2019, 5:01 a.m.