Run_simulation/compCL_test/test.R

In jiji6454/Rpac_compReg: Compositional Data Regression

library(glmnet)
Rcpp::sourceCpp('C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/program/compReg/src/constraints.cpp')

n = 50
p = 30
rho = 0.2
sigma = 0.5
gamma = 0.5
add.on = 1:5
beta = c(1, -0.8, 0.6, 0, 0, -1.5, -0.5, 1.2)
beta = c( beta, rep(0, times = p - length(beta)) )
pos_group = c(1:3, 6:8)
neg_group = (1:p)[-pos_group]
intercept = FALSE

Comp_data = comp_simulation(n = n, p = p, rho = rho, sigma = sigma, gamma  = gamma, add.on = add.on,
                            beta = beta, intercept = intercept)




fit1 = glmnet(x = comp_cvm$compCL.fit$Z_log, y = Comp_data$y, family = "gaussian",
              #weights,
              alpha = 1,
              nlambda = 100,
              lambda=NULL,
              intercept = FALSE, thresh = 1e-5)




m1 <- compCL(y = Comp_data$y, Z = Comp_data$X.comp, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
             pf = rep(1, times = p),
             lam = NULL, nlam = 100,lambda.factor = ifelse(n < p, 0.05, 0.001),
             #dfmax = 20,
             #pfmax = min(dfmax * 1.5, p),
             #u = u,
             mu_ratio = 0, tol = 1e-10,
             outer_maxiter = 1e8, outer_eps = 1e-10,
             inner_maxiter = 1e3, inner_eps = 1e-6)







fit1 = glmnet(x = comp_cvm$compCL.fit$Z_log, y = Comp_data$y, family = "gaussian",
              #weights,
              alpha = 1,
              nlambda = 100,
              lambda=m1$lam,
              intercept = FALSE, thresh = 1e-5,
              standardize = FALSE)
y_use <- y / sqrt((sum((y - mean(y))^2) / n))
sum((y_use - mean(y_use))^2) / n
as.matrix(fit1$beta) - m1$beta[1:p, ]




m1 <- compCL(y = y_use, Z = Comp_data$X.comp, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
             pf = rep(1, times = p),
             lam = NULL, nlam = 100,lambda.factor = ifelse(n < p, 0.05, 0.001),
             #dfmax = 20,
             #pfmax = min(dfmax * 1.5, p),
             #u = u,
             mu_ratio = 0, tol = 1e-10,
             outer_maxiter = 1e8, outer_eps = 1e-10,
             inner_maxiter = 1e8, inner_eps = 1e-14)

fit1 = glmnet(x = comp_cvm$compCL.fit$Z_log, y = y_use, family = "gaussian",
              #weights,
              alpha = 1,
              nlambda = 100,
              lambda=m1$lam,
              intercept = FALSE, thresh = 1e-14,
              standardize = FALSE)
``



















#compCL

m1 <- compCL(y = Comp_data$y, Z = Comp_data$X.comp, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
             pf = rep(1, times = p),
             lam = NULL, nlam = 100,lambda.factor = ifelse(n < p, 0.05, 0.001),
             dfmax = 20,
             #pfmax = min(dfmax * 1.5, p),
             #u = u,
             mu_ratio = 1, tol = 1e-10,
             outer_maxiter = 1e8, outer_eps = 1e-10,
             inner_maxiter = 1e3, inner_eps = 1e-6)


m2 <- compCL2(y = Comp_data$y, Z = Comp_data$Z, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
             pf = rep(1, times = p),
             lam = NULL, nlam = 100,lambda.factor = ifelse(n < p, 0.05, 0.001),
             dfmax = 20,
             #pfmax = min(dfmax * 1.5, p),
             #u = u,
             mu_ratio = 1, tol = 1e-10,
             outer_maxiter = 1e8, outer_eps = 1e-10,
             inner_maxiter = 1e3, inner_eps = 1e-6)

all.equal(m1$lam, m2$lam)
all.equal(m1$beta, m2$beta)
all(m1$error == 0)
all(m2$error == 0)


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

comp_cvm2 <- cv.compCL2(y = Comp_data$y, Z = Comp_data$Z, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
                       pf = rep(1, times = p),
                       nlam = nlam, trim = trim,
                       #nfolds = nfolds,
                       foldid = foldid,
                       dfmax = dfmax,
                       mu_ratio = mu_ratio, tol = tol,
                       outer_maxiter = 1e8, outer_eps = 1e-10, inner_maxiter = 1e3, inner_eps = 1e-6)


all.equal(comp_cvm1$compCL.fit$lam, comp_cvm2$compCL.fit$lam)
all.equal(comp_cvm1$compCL.fit$beta, comp_cvm2$compCL.fit$beta)
all.equal(comp_cvm1$Ftrim$cvm, comp_cvm2$Ftrim$cvm)
all.equal(comp_cvm1$Ttrim$cvm, comp_cvm2$Ttrim$cvm)
all.equal(coef(comp_cvm1, s = "lam.min")$lam, coef(comp_cvm2, s = "lam.min")$lam)
all.equal(coef(comp_cvm1, s = "lam.min")$beta, coef(comp_cvm2, s = "lam.min")$beta)
cbind(coef(comp_cvm1, s = "lam.min")$beta, coef(comp_cvm2, s = "lam.min")$beta)

#GIC.compCL
comp_GIC1 <- GIC.compCL(y = Comp_data$y, Z = Comp_data$X.comp, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
                       pf = rep(1, times = p),
                       nlam = nlam, dfmax = dfmax,
                       mu_ratio = mu_ratio, tol = tol,
                       outer_maxiter = 1e8, outer_eps = 1e-10, inner_maxiter = 1e3, inner_eps = 1e-6)

comp_GIC2 <- GIC.compCL2(y = Comp_data$y, Z = Comp_data$Z, Zc = Comp_data$Zc, intercept = Comp_data$intercept,
                       pf = rep(1, times = p),
                       nlam = nlam, dfmax = dfmax,
                       mu_ratio = mu_ratio, tol = tol,
                       outer_maxiter = 1e8, outer_eps = 1e-10, inner_maxiter = 1e3, inner_eps = 1e-6)

all.equal(comp_GIC1$compCL.fit$lam, comp_GIC2$compCL.fit$lam)
all.equal(comp_GIC1$compCL.fit$beta, comp_GIC2$compCL.fit$beta)
all.equal(comp_GIC1$Ftrim$cvm, comp_GIC2$Ftrim$cvm)
all.equal(comp_GIC1$Ttrim$cvm, comp_GIC2$Ttrim$cvm)
all.equal(coef(comp_GIC1, s = "lam.min")$lam, coef(comp_GIC2, s = "lam.min")$lam)
all.equal(coef(comp_GIC1, s = "lam.min")$beta, coef(comp_GIC2, s = "lam.min")$beta)

cbind(coef(comp_GIC1, s = "lam.min")$beta, coef(comp_GIC2, s = "lam.min")$beta)