Install the package
## on CRAN
## install.packages("Compack")
library(Compack)
Regression with functional compostional predictors
Generate data
library(Compack)
df_beta = 5
p = 30
beta_C_true = matrix(0, nrow = p, ncol = df_beta)
beta_C_true[1, ] <- c(-0.5, -0.5, -0.5 , -1, -1)
beta_C_true[2, ] <- c(0.8, 0.8, 0.7, 0.6, 0.6)
beta_C_true[3, ] <- c(-0.8, -0.8 , 0.4 , 1 , 1)
beta_C_true[4, ] <- c(0.5, 0.5, -0.6 ,-0.6, -0.6)
n_train = 50
n_test = 30
Data <- Fcomp_Model(n = n_train, p = p, m = 0, intercept = TRUE,
SNR = 4, sigma = 3, rho_X = 0.6, rho_T = 0.2,
df_beta = df_beta, n_T = 20, obs_spar = 1, theta.add = FALSE,
beta_C = as.vector(t(beta_C_true)))
arg_list <- as.list(Data$call)[-1]
arg_list$n <- n_test
Test <- do.call(Fcomp_Model, arg_list)
sparse log-contrast regression with functional compositonal predictors
m1 <- FuncompCGL(y = Data$data$y, X = Data$data$Comp,
Zc = Data$data$Zc, intercept = Data$data$intercept,
k = df_beta)
plot(m1, xlab = "-log")
![plot of chunk FuncompCL](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAABiVBMVEUAAAAAADoAAGYAAP8AK2sAOjoAOmYAOpAAZpAAZrYAzQAA//8gAAAge9soAAAogf8qAC0qidsra7A6AAA6ABc6ACE6ADo6AEo6AF46AGY6AIQ6Kjo6K5M6OmY6OpA6ZrY6kNs6nf9Krv9LAC1MAABTa7BVp/9mAABmADpmAGZmOgBmOjpmZmZmtupmtv9m2/9pkNtrKwB+0v+BKACDVQCQOgCQOjqQOmaQkGaQtpCQ29uQ2/+Tud+jvoe2ZgC2Zjq2voe2//+4//+5//+7vsm8vtS+o5O+vqK+vqO+vr6+vsG+vsi+vsm+vtS+wdS+yd++1OrBvr7BvsDD2vTJvr7JvsnJvtTJyd/J3+rJ3/TJ4v/Nvr7P6//Q6v/Uvr7UvsnUvtTU6v/V9P/bkDrb///fyb7fycnfydTf6t/f9P/qtmbqyL7q1L7q///038n0////AAD/AP//tmb/25D/28X/4sn/6rv/6s3/6tT/9M3/9NT/9NX/9N///5z//7b//9v//+r///T///+qEIEJAAAACXBIWXMAAAsSAAALEgHS3X78AAAPUUlEQVR4nO2diZvcNhXAvSxQdssdKIQjXc6kpV2gsBvIwpZwdkuhAUJLOdKmmxYCOSY7HAGWnUR/OdZh67A0I9mSped5v+/Lt19mLEvj3zw92SPZFUFAUeVuABIGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgPGSmFnz1wlZxer6olbYzQHWYWPMOqMzJ/ulkVSMFzY/MItHmmhZZEeDBZ2cePKZRphFzp9IgpLwVBhhCz2qi1y8p5OgKGwJAwXlqQs4gKFlcW2FWWDWMKUQYfneGY9sQuxqrGDETYWnkJWgcKS4xs7fqCwhEQ1JUBhqeihauZC2Wb4ibMYYXRPxNZVmLsLdArpenEyOMIWezu9y04Qi6kwIasY3iWeXTroXXZKWOIqniYJ5rAYaKriRpQJChuGFlfpNElQ2ABaVQkjygSF9YXbGs+UAIX1Ybu1NXrVKCwMmbNy2CIoLAgtZ2VqAwrzRbWVsRkozIftQmwRFOaBzFnZbREUtoq8IwwLKGwJTUdYjC2CwtwoaStvQ3RQmJVyBhkmKMxCsbYICrPQJq7ybBEU1qG9pJu7IQ5QmAbVVa4sCgpT4Lpyt0KjM4cahbWUomv58j0UJsirK9oSyzURRk+7MqQuX0lamUFvJys7KiPb6qNJKT3o7WRlR4TpGqeqGKuw1l5Yel3DIqqzt0FvJys7FtvbKXWlWNe41sJqXWl2nHAJ6voKi75yi5F8sfC6Cotva6R13espbBZT17hL8AcLm59zrOcrV9gsVnjluFnCUGGLvV3296R7O7dChc0ihFfG+1oMXjIr7s8B5eZgs4Hhlf0OJOsVYfVJV29d2VWJZgx6m7Sr0iHkMKqrh6+y7uuzNqNEdnk3XFc5pgRrIoxdgAoLr6LiShJLWNE3BwvWVVTrddYgwoJ0lfNV21RQXp68MH413k9XblWbDkkqExcmdPn4ymprqSMNP2EnrqH7gF2PgPita4WurBl3eThZ8BLmvDtRyTcHa6dpLPWVL65CTQn8hHWvO7UUenMw+UPyEl95bAUHlYZfl3jskkIKvTmY/CXZ5StTL9jflMBTGKwcpoSX3VcOWYPiSjI0h/XfdTpWhNfYtoZ1gSaDc1jvXadieXiNaiuqKcHwHNZ314lYFl6j2YobVBp+Eeb8CWXArpPgDq+RBhnpTAkmdaVj5gyvEWQlV8WZTg7T1jPovpLrGkUVx0vYYr/4UaI241rvDpPqGimuJBPJYdqUa+1hQMl0ja6KM4kcNnP4SmYriyrOFITpKxoaXyltpdmxF37CFntV+JN/RxI2s/pKYitTL6jhN+hgV+SPA42NI8zQxX2l0JVdFSdgWB86uB9DmC28EvSF+QOrBXaEGevxqK8ktmLvcgCgc1jXV3RdZdkisEeJHV+xdRVni0AWZqav7ci6SrRF4M6aMhf/b0etraBBhonnpand+LsehKGrirkEtlxZFJBX643wqiKuMC/bFvHtEg+L+sXZ0BUvvIq3RSBerdfDqx5pRPIFwRYBOEo0dUXxVfAgwwSaMO0Z1LSGCL7AyKIAE6b44qddw32B0gVNmPQlzpIH+wKmC5YwZbgh9j3QF5zMJfERNj9XbRzk/3mlE17DfEG0RbyE0Xun0B9YrMLoi3TUb7mUH1eYDK/2muEQXzBtES9hXNThlksYe33+dOiuw+jqGuILrC7vCKs5ft95u7D5hVvJ7zXV+FJ09fcFWJdfDju7yK5MHduudJxd3LhymUbYhYT3mmq7Q7nLNdUVY5S42Ku2yEnKe011w6u/L+C6goTlGiUKX+rvk319gdcFQRj3pf2c3NPXBHTFE5bsXlOzTjt6CpuEruIjjA83quG+JqKrdGE2XX18TUaX34mz+2Y3JO1dtamvTs8a7mtCuoYP61Pe83dmmxga7GtSuoYLS3dXbZq+ursI9zWsFcVRaoTNrOEV7Gti4UUiXOlIc1dtu65wXwOaUChl/oA5s0+TD/Q1vfAihQqb2cuG+ZqkrjKFzRwXSUKETVRXkcJmjoJBvnrWXT7lCXP4CvnBcrLhRQoU5ngqJYaXoDRhw31NObxIccLsvoK6wx61QqIoYZXDl/8eJh5epChhlf0Z2BheGuUIw/DyohRhld1XQHitha5ihFX24Qbq6lCEsMruKyS8vLeETgnCKutwHgcbVgoQZvUVNHV+jXwVIMzuK6Siift68OCB8r/swmy+wsJrur4eMPTXcguz+gqpZZq6bKo4mYVZfDnC60EGBnz4vqysObuwji/bZq6PkDi8xrbnU0teYV6+3B9j1O4wafj57zKnMMv5ctfXko+RN31F60jdZZ5sUF7LKKybv7rpa5muUocbMbKn+zuQT1jVvQmsuQlEXd7Y46obUsYNRbIJ63SHpq6lnQtwXbZesGuK0rn9Sy5hpi8zvJanAri+LGnOElQU/tyQGUN5OZsw/YFEQbrA+rKrkv+ftbjvq5RJmOYrUBdQX6otLaikJuLxuNM8wlRfwbqg+dJGetzUTEVuaZpKMqzvs2RW86W/tfI0Bpiu2pTwoqrqbNbG1ZMqm5vdgMuxoE/xpYeXx1knAF8zA8egghh22q1aTZvbDUqZDEtmW1/bpq5VlZXpyxREeC9osWG3I6SYmioFpbbxI6zx1cldq6oiZfnidioDaaN9adsJ7+yEJ3NPhFhvUTP+klnuK3SoQck73OiaEILqIy81SRttQdOEoI0ozYjprMPoo0Tmq1xdWgRo/ZcaMWJjbRTnd+VXJqj2pZWONEYWRnt4I3V5/QqUTJfZS5k5xnIotfc5pirz4q3UZL7jdqS2Snk5lrAVNwdbdinaIzsXxWY47jTmxn6qlvsXZySUcYU9PDr62Y967q4Zj8bl9LW7hNw7Ojr6xV3xCr0SwGqa2+5xbJS+Lsvp8L0cs55m16+IybF1LDfqlY7T1/5affzayl3aOa4SCHvIjt6dm/KVs0sHZP7UAX0q4cqnOT168za59+p9yzvtXswTnkc3bjqKdDi0ft5xz8MON777uZ/4fb9M5p9/Nr6wO6/8pY6wx+/eli+dbBF2rOqmroyw09fvk0dv2T5Psxeuzizy5m1LEZPF/oHt5bGvdMy/0i/CFvtXknWJj27UXaISZOIgr+4S3REm93K8pdfnLYyd4HY/8rgRRvvPLyr/f8hzR33EVvUSxzvRchgT0WQqKuz0t7e1KPv3D77B2vTTH+ot32V/1I7yv786OvqTo/38aeUiwIT8xd7ztEt8xRBW5zX6hdFSKZl/8oNV9dkDmVM5o17pqFs//9RHX5ANrY/WvWvk8Ts36Z9lzC/ciiWM5aY2x7BBB0XmsX8e/fwubZPaGYjt6R/lgfH/+fUfyT9eftvafnGzf/5FbvJhnYdrOb97VxdGI45+adRUKqr88BeUbMgYdZRY9+zz89/+qvZafcRoGmgPnB0+3NoZ0JoGnpvaHNMV9reXvke7ybfuzpV0K7anj1NQOv+//7Lu3/7we1v75+f414s9PbTJhyIPm2nv4TVW/2PdI6vyWx+R2ZAzqrCHr96ff/pLL2iv1V9Nv349bpdIRG9FD3TdLPL4z+Ig1keadZOv359/5sdam+rtjQijYfGvl9+2tL/x1YwcWJ00Dz9P056lN6nLG6mUxub/PvH1tqWCcc/D7h299M2PqS+cXq/7c3rAxhfGc0xzHtamlTqW3/v9L+/S75YujG1fJwA1Wdc5+Og3tvY3J2DNoWZ1sjxcF7Hk68fvvNFJpXQfH5ItFYx8pUP07Ar1J/UfOUWBC+u2pEVEmN4mtj1NJtpTSx7deMNr5EfrXJKH2W4Y2ikhb6Le0nGFNT2Fyp2bHjksJmKU6I7XU57DtDbx7Wk3pfZPp9fp8fVoP63TnYf5bhiKMNFEo6WjCjOPkuhLaH+wYpQYE/ZtX+KLHXyjTc3B0yNMHGiP9rfdsK1asRstlbp8jSusc2lNJA+P87CIyG+7QxqLFr1NzfYnFXu6pICeONFxwur2LxXW7EZNpW2VZktTCkNSkE5YnJ1hoaBCKAxYIRQGrBAKA1YIhQErhMKAFUJhwArhzCdgoDBgoDBgoDBgoDBgoDBgoDBgoDBgoDBgoDBgRBEmJxvbV9BYOWm3NKaW+RUKqEn+ch9QkywUUJPcNKAmWcinphjClMnG9hU0NtjEC7aOgB4YY0XB6kIBNZF2JVNATbJQSE3tpiE1yf371BRDmJxs7FhB40JMIHzmqs/aOr1QUE3NSqagmppCATXJTQNqkoW8aoqVw/hkPscKGhf8O8jms1/yPiq8UEhN7UqmkJraQgE1yU0DapKFvGqKJExMNqYdo/c3cn5uo13B6C2sKRRSU7uSKaSmtlBATXLTgJpkIa+a4gjTJhv79/liXV1YhCmb+tUkZ1AH1GRMuw7NY4GfKSCPRRolqrWEfrjAHBaYpNWVTAE1GdOuy/lMMYRJX7QnWFz2aqfsNGh36jmikoUCaiLtCD2gJlko8DPxTQM/Ey/kVVMMYWKysZghveHZEfBNaaGwcxZRKKAmfuwDa5KFyvpMeKUDGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGCgMGOCFOSaTHe+4Z5nJmWRy6iGdSAiCiQqrj3+YsNpw7JalYTLCzi6ylTr1n/d/7YA6ad7ga6EW+9+pqp0TNjf08Dk2hZ1u+uxus1bK8qiSIpmMsMMdNhGT/tk4UCbesrVQ568u9rZqNVvs1cMnWPjRTatd8X7g+qV8TEUYVUQfFFD/Wewf0IykdIn1q3SRAf1H/VA1h7ts0712GRKUPnEqwugfYar+Q2OtFXZIlzWqwl48YE/taMLqkC97RGHjsCrCzi6ydGZEWB1aIsLE+yhsLFblMHYT+6cONGEilbEcJt7HHDYW/ImPW8oo8QP7fJQo3jiu6GBQE/YcW3Ow2GOjRP4+jhLzIR6DElgKSI84NWF11PAVO6HHH690IGlAYcBAYcBAYcBAYcBAYcBAYcBAYcBAYcBAYcBAYcBAYcBAYcBAYcBAYcBAYcD4P1aF1Le/E254AAAAAElFTkSuQmCC)
betamat <- coef(m1)
predmat <- predict(m1, Data$data$Comp, Data$data$Zc)
Cross-Validation tuned model
linearly contrained group lasso creterion: CGL
k_list = c(4,5,6)
nfolds = 10
foldid <- sample(rep(seq(nfolds), length = n_train))
## cv_cgl: Constrained group lasso
cv_cgl <- cv.FuncompCGL(y = Data$data$y, X = Data$data$Comp,
Zc = Data$data$Zc, intercept = Data$data$intercept,
k = k_list, foldid = foldid,
keep = TRUE)
plot(cv_cgl,k = k_list)
![plot of chunk cv_cgl](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAA1VBMVEUAAAAAADoAAGYAAP8AOjoAOmYAOpAAOv8AZpAAZrYAZv8yzTI6AAA6ADo6AGY6AP86OmY6OpA6Ov86ZrY6kNs6kP9mAABmADpmAGZmAP9mOgBmOjpmOv9mZmZmtv+QOgCQOjqQOmaQOv+QZv+QkGaQtv+Q29uQ2/+2ZgC2Zjq2Zv+2///T09PbkDrbkP/b25Db/7bb////AAD/ADr/AGb/OgD/Ojr/Omb/OpD/ZgD/Zrb/kDr/kLb/kNv/tmb/tpD/tv//25D/29v/2////7b//9v///9ja24oAAAACXBIWXMAAAsSAAALEgHS3X78AAASlklEQVR4nO2dD3vbthHG1WrL5rFN1s5xm66V02yrtZjOFjZdm3XRrEXC9/9IIwCCxJ+jeCBBiife+yixDQkExR/vcDwA5EqwSGl17h1gxYmBERMDIyYGRkwMjJgYGDExMGJiYMTEwIiJgRETAyMmBkZMDIyYGBgxMTBiYmDExMCIiYEREwMjJgZGTAyMmBgYMTEwYmJgxMTAiCkFsP07+V9RvAPf+bWQemyp470BbAVXBLQCNXz8qfjxo+gqUlXf9SmSm/O/6i78FFCEVQJgO9n4ofzSv/4CvaN+Ad+RdfbOwTqUx9fbCrBhsC2olbDo/SOmSJb5whVJil6pOmn/3VWE1nBg73/8b9n+vvzOh3fQO0IfYeCdndxn51urrbgfBjYMtQW1EhZBtYBT/fivYEu4onLffg4sTFT+pLMIpVQuET7rq70Cziazv+4hBYDhLQxoxS/a//Le939AUWnogS/FFZUQ/wcBQ+waVsn6sHL3w2OosQCnvgF2/OmDXapcom8V4YbBthAGVnZ+j+qkOF2k/LRnP7gi8euHsA8rmwh2DShCKhUwufe7wMo1FqBvMSiLD15xUfwnPCzehuG2ED1Yq8uFjl6ffqzcGAAMdzIhlQrYDvzeGsv7D3AddXYH8joVYMNwW0ArQdEBcLkpgenAFPi62JClWxNYGNw5v4N4qd7Ji6iQFoaLCt6H/g8okmfE8eceRQII69FnHE7J+rDy0iLcBd27QbvWcoW2C3snYMNQWzivcwgvgIAiuWd+XVwRdB2GrYkTZzqIiYEREwMjJgZGTCmAfXx2L/9br++hd96spTZtdTZecbAVXBHQCtjwq/X6TnQWqar3PYrk5vyv+hh+CijCKgGwR9X4izvx5ndvoXfkL/A7ZZ2Pz5yD9WITbAXYMNQW1EpY9HqDKZJlvnBFkqJXqk7aP3QVoTUc2Ov138r2P35ZfucX98A7Ui/uoDqPcp+db6234nwY2DDUFtQKUATVAk71V+GWcEXlvn0LYPz4DHQ+fZTMJYJnfbVXwNlU769zSAFgeAsDWvGLPn75F9//AUVlA6EvxRWVEAOXiNs1rJIBK3c/PIYVFuDUN8BePXeKpUv0Dx+wYagtjIF9fLbRJ8XJIu2nXfvBFYk3z8M+THa5/q4BRUilAib3/jHoSDUWoG+pUT73i9fPw8PibRhuC9GDtbtc4Oj16cfkxiALw5xMSKUCpo5NsA8ay+vncB11dodyOxVgw3BbQCth0bdhTaBIV+4BTAemwNfFhizdmsLCwM5Z1Qn2WfZOfkSFtDBUVPA69H9AkTojXoXnTWeRatT7UugzDqdkfZi8tAh2oerdWoABF0rlVoArAH/DYFs4r/MivAACit6EDeCKBHAdhq6JEmc6iImBERMDIyYGRkwMjJgYGDExMGJiYMTEwIiJgRETAyOmIcBWrDE0IrABdVltmhxYMWCLnVqPufF5iIERE7tEYhoK7PD1vTjcrFZPwnkZDKyvcqOikXkvATDJTOz/hKzLLvG0cp+WcI9ZAmD7L95qS0PVZWAnBNBKDuzmkx++lxb2ReAT2SX2UV7+s2mlBibE8XZ1JXafhjOLGVicGvPKDapRgEXWZZd4Qrm2MAY2fwWB4djArKADmfViecprCyvmZGGsNuV1tDErYOwSQ3nRxrgWtv9MuT90lMjAQNX+cGxgx1s9V3wX5qbYJUZoMmAm2EBnOligmg5MzMzC2CW6Ai6XhwMDzKeWTNVzHzZMtj9MAuz4ss9qJXaJWCUH1m5F3XVZ3bI7sNHD+ti67BItwfleBjZnVcnexMCOtytwFgCmLuukRgF2vL0u/3+IJMbAMBoFmA7rTwX37XV9sUs0qjsvaIR5PhbGwCy1jzBzHzZLtY8wO6+zRoksSzCwLNOvrCj0yyY2ODXVVdcXu0Slxhd6HZgilSlawDTSM6SmGJiRFXQ0L4CWtLKmFqemzqZwDkeW+bxU/5VZvLgPO588YFmAa5zhFaX9U+AD7BLbBHVgPq+xhlcqbwn5SwZ2Qv6km4aWjgrHG1453JTFERbGUvLmcNjWJcae5na4efIPBhYpB5jrDicYXtl/Bhkfu0RIud+BiaYDq6xrKDAZcsjVRJz8TSQ/g2gyG6bzmhcwlj+iUpmXKCYFxosh8HKBZTUwkWhOB7vEZAKGwDITc6QbD2NgKeV1YJnfgaUA1n5tLOIXQyxd7ohKFnRgaQYw28WLIdCqTcseAnMD+gmARS+GWLRLrC6XTQem/KFw/eHsLGzZwOwhsNodOv5wbGDRiyEWLXdEJQzopwA2Tt2Lk92BaWD1zA0w2TtelHi6bqjlukT3tg7GunSAmBRYqa2al3gVt4MMzBWcoC/GAJZy5u8ClXv+sKjyvXV+I7FL5Ln1CZQ7q8CMO2wf/xroEtXM30iPyC7RlpOQ8vIbIwDrJQYmBflDL78xb2ALlElF2fmNZgKH8ADJSFK/wLm/aJf45J9fRU7/ZWCVcsgfwiPMipRFK9wYNujYf/EWyD4h6vpalEtsUlGhPxTBgKUxrXZaAh/Wl8B4PKyPmux8rvMbWcXMy85rWkL4LtBXhIVxWN9HuZPgMOlecMDSMDspXtA3msLxr/pyWXgDlljzEhzWj6vcdolFltkJjhaX2CVe0DemAHdYLa2Exr86rUtq8IK+/Wer6y2Ph/mChlPs8LAB1lx3oTY8dDGEHHHefnrPDxoIldsu0VlOJNoGLBFKMadjd81zOgBB/jDz8r3VdRfWH4p0czrQFrYAYOBwis0LSkVhN44D9nBqfZgk9sB9mCvbH7rLYcN8L94fCmwf9tXd7irRiPOlK7cF+0M3ku9IRfnChvXmOWERWq5LFL55OXdzCBL0UccEG9aXL3BdrBF69crigGWBP3STG2IEYHIN8261uo7b7wW6RMcfamCZCPyh6xJj/KHgAczkslJRVfdl8oeFDywu2qg0GBg/ysOVPZxS5za89GHuJXtjtp8i0yHFc+uD7HzedF7Aeoe4q69GERb2APVh/CgPW7k9/mVFG+buNm5yo9fZGwEMDOt5fZgtezgF4OUlN+L9oYgCBj2YlB/lUcnPzju4Wmdj91BMH7ZJsekLBSYsf2i5Q+Fl5+uea1xg/bQ0l2ilD+1wI0z2xmajHDGw4fKz85njDxtgfZMbjmJcYuQKsUW5xCY+tHOHQXa+V3LDEXJ45Urw+rCT8syrCJK9lnkNOgS8PmygPH9opTbsq69gfUNvxawP4/EwWOq6SwviZc+zGegPBa8PGyZv7rzhlDLZ64snkg5Ukzu0Lr2EvjdAYYXxA5IbjjisH6g6PDSDKKF5xc5kOykMsMPNlRxFibzrw8UDc7PzFS89YCkMsGbaRhLzEjhg22uV4+WwPpB/rVwNWApRwDe3SSDMjVW+1vc550k4gYw/zBx/aJK9VmyYyLwEGphcG3ZyEk78pmnLufrStPL6VkRwdj6NUC5xoy7EtqBLLDu3T+4WOYDZ+EPRRPOiunWUbV6TAzvcrJ68lZEHINm5SZxLmltvJTcqfygy69a9ReHloYZeKztKsRhCGt+SgEl5AylNUqolO59MiRZDPPwm7OAu0yXmYe9VJzdUuJEXzpyNpOYlIoCBU3CE9JfqDWA1xGUCE82sKBMd1mG9/ayAUcxLJAAWu+kLcIlB7rC+OYDuwKzMRlrrkmJgMbKSvc3FV52NstKHI5mXSJdLXMyjPLxpbFbcUZuXsbD05iVwF85JbyFLV07uMPPcYcWrcCeKjiC8hT3wNDcTzNep+Tq50ZiTsbBRzEvggaknJ0JayGIIKJivfpPJXjc7L0b8mkhgD22LwxY0VduYV+MDjV+Uud28/cGwSYWc5tbafS1hMQR4rWz95mbn018rO8IA251Ye7mU5UaueSkoxryUdfnZ+dE0OEq89MUQtnlltXX5ucNiGvMSPKfjtNxZUXYAX/2eu8Hh2OYlGFingvVe1pQAM448ZAlsrGKuwxb27BVokk1mL3aor5XFVOYl0MB2KzWunGDTdICJcJKNijZcb5h6GluXMMDULIBtZJ6DuEu0zSuzzCuz0XmpqPH9oUBPwpETOxJver4C1zdkTV6+tK5cXSgXTmZjAvMS6Dkdq00qYERcIkSrXghrZq3ZmY1JzEvgg47tUvow6Lor88J6gyfFispY4aPE4+3lD6+A112ZTu7W5lUYYGNN2zgpvg7zFN5nw1m2HGY2JrQuKexU7QU82t4JDG1aViQPZDamNC/BwIwcWm7nVf+Rh5mNic1LnAHYDNUextt33MhN12XfK2pa65JaPLA2Wpl9O1FR56HszMb0tMTSXaIbFVq25ZpXs1yomPxC2dfks6bmA8xN7lpdVeb9bYUWzXXXWcxLLDesdx2hfbeh4O/zXnf5Grx6pd38Zgss92l5KY2ssGl5zvBsvrDSYAvTN13B1z23Swy6LT8q9C+SNazCCwzPBg0BTN1Vpb0PO7Q9fXaGwNq6rSCGt+aGFsVZw3hfGAt7iL5lPWbTkypvVHnBLlpunDHdAGWX0BNJ46HNBtiJPkvR8W5xKDIoKpyFeYmoPmx7Kqyf56M8bLvyOig7yHB5FdbwyfTjXV1CAtvG3/L3rH1YhwvM5NVWEGS4UeFMwnhfyD4serjZ3nRzko+t3FEFxHeBldfzza0xLHei/JxoiQRRIm71St6iIbvesvHmWjjg1OYFXcOyb/k6L1hSie4i0DW3vrlJXWGtskoovfF1i0VVLtC8XNmmVZw9udulNPfp6Fy9kgnjhQphpkgIa94f8PfJN+2/s2bja9CmdHoQtDgHleUIZ0pLTGVh8GGsz3gbaGYDwPzdpWrjACo4kzG/bsvR4NQUZvVK91FNrDzLO1yg5wWDIf/Y4ziZpsnW24drNEbla53b5gqpCOT3WXOmJaYBVkcG8pgW8v+8emU5fFixgEwoUfdh7cbsW1Xu3INy9oZVa6LxMFwcYfVROQikpU87oaKAXGAYudOgJaayMNsoACPR5tYSdBgLwQYZNabCu7qyXaBlVFRIVZrEwpxjiTvkvVSs/Ysq16oKmkbliHTQoY55/fygUusTpOjEFSc1BTB9VNVNBauX/bfDsnBhtEu0GJL9vn1j6wtApTWNhZmDVzivjr/7yMIkAqgDvst8NEnQ0e/gRwCqaMvhlfAjl6VJLGwYinZzDPzd+iIZOZqtS7SzRX4YDvFdiiYKOlJpwO5ciCaf+XvueYnUxcCIaalz68mKgRETu0RiGhMYqAIuTqP1mBufiUYExpWmrsTAiFViYMQqMTBilRgYsUoMjFglBkasEmcryImBERMDIyYGRkwMjJgYGDExMGJiYMTEwIiJgRFTP2DyzlR60XPzW0Slw80qXOLeJrMqPqKlplJES7t68xEtNZUiWpJ3PunxnbT6AWue0xLxxJb6o/JoPlxhq5nHS0c9G6aqFNHS/um9+Si+paZSREvydob7z+/iWqrVC9jx5V3wW0QldVPop8gTa//HbzaRLTWVoloS1UejWjKVIlraSa4KVWRLSr2AqTtBbNzfIiqpO3i33TXT0/HlD7exLTWVYloqpU0kpqW6UmRL+qORLSn1AiYNWp8dzW8RleQ9WrBf7uG66o4iWmoqxbRUdi36+U4xLdWVoloyt+WNaqlSNLDtauV7+m5P7FZCno2yUvlREz/gW2oqRbQkf1ofxX8nXSmqpcONdS/RyH6sf1jfP/CI8PcPaqZe/fVwLTWVIvuwab5TaZX2xqcAJu3/+P2981tEJekQ8FGi5d2QLTWVIlpqfFrkd9KVIlpqeMV9J63e12Gl61ZBLf7ptE2l6OuwyJaaShEtTfadtAPYRH8nLc50EBMDIyYGRkwMjJgYGDExMGJiYMTEwIiJgRETAyMmBkZMDIyYGBgxMTBiYmDExMCIiYEREwMjJgZGTAyMmMgDa5lb9nDdPumsmVjWzHmU8wpJ6EKBlcc/DlhJOPWejaOLAXa4UQt3yh+//fOdZGLe0It7ji//ulpd79T80u13aka7/Og3G7P4B3go6yx1McC212pepvzxyZ01D1ct7nl6f7y9KtFcqdLtE2V+8qOrTfV+r6U/59ClAJOIymMvfxxf3skeyXKJZalccyD/ST4SzXajPnpbr0qi4hMvBZj8UZEqf0hbq4Ft5SpHG9jf7ySw2qy2ehUkA5tGXRZ2uFHdmWdhpWlVFla9z8CmUlcfJn/sP79zgFVdmerDqve5D5tKMsiTi66aKPH3L3WUWL3xsJLBoAPsO7UE4XirokT9PkeJ55OymejrYCIe8dKAlVajF/DEHn/OdLDGEQMjJgZGTAyMmBgYMTEwYmJgxMTAiImBERMDIyYGRkwMjJgYGDExMGJiYMTEwIjp/7HKnKaVIcF/AAAAAElFTkSuQmCC)
cv_cgl$Ftrim[c("lam.min", "lam.1se")]
## $lam.min
## lam df
## 0.005003941 4.000000000
##
## $lam.1se
## lam df
## 0.01203433 4.00000000
beta <- coef(cv_cgl, trim = FALSE, s = "lam.min")
k_opt <- cv_cgl$Ftrim$lam.min['df']
plot(cv_cgl$FuncompCGL.fit[[as.character(k_opt)]])
![plot of chunk cv_cgl](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAB6VBMVEUAAAAAAA0AADAAADoAAFgAAGYAAP8AK2sANHcAOjoAOmYAOpAATJcAZpAAZrYAzQAA//8XAAAhAAAoADookNsqAC0qidsrAAAra7AwAC0yAAA2AAA2kNs6AAA6ADo6AEk6AEo6AF46AGY6K5M6OmY6OpA6ZrY6kLY6kNs6nf9JADpJtv9Krv9LAC1MAC1MADRRADJRvP9Ta7BVh8lVp/9etv9mAABmADpmAGZmOgBmOjpmOpBmZmZmp7tmtv9rKwBrLQBro6N+yd9+0v+BZgCDLSuDOmaDVQCDtr6D2/+J2/+QOgCQOjqQOmaQTACQZgCQkGaQtpCQ0t+Q29uQ2/+Tud+jvnSjvr6nZgC2ZgC2Zjq29NS2//+5//+8vtS+rX6+vqK+vqO+vr6+vsi+vsm+vs2+vtS+yd++1OrBvsHEvtTJvrbJvr7JvsXJvsnJvtTJydTJyd/J3/TQ1MnUvr7UvsPUvsnUvtTUyb7U0ZDU4t/U6v/V9P/bfCHbkDrb/+rb///fyb7fycnfydTf6t/f9PTf9P/qtmbq1Ifq1L7q25Dq///02sP035D038n0////AAD/AP//nDr/nTr/tl7/tmb/22b/23z/24P/24n/25D/6tT/9J//9N///7b//9v//+r///T////0IyzyAAAACXBIWXMAAAsSAAALEgHS3X78AAAPRElEQVR4nO2dh58cNxWAB0xCMdyaXm3AYHDoLXYC4Sihn0NJlrJgSo4AAQwGEwLBsIkpBsfg3AFHWXN39vylqExfaUdPbfRm3/f7Jbu+HY00+vZJGq1mJssJVGRDF4CAQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQQcKQ4SIsI0IQUJhDWkIHCUMGCUMGCUNGbGETHQ4ZrRWRhPV7IXtmpNkkkjotaQqrIW0dUhcmoXCrwCGsgLQhEyZZb2cIheVr3URGErZT4rC/ZdZRW+wI21mBZT7r5SylJtFe4ho1kSkJUwEytw7OUhdWYRh0o3eGRlhNn7pxO0MorEarbcR9GmphEn20jdHZCIRJ9LEWtRjBGY0wjjrWxuVsVMIk43Y2QmF5NZBs/mksg5BxCivoxtoYnI1aWD5CZ2MXli8PRXAriyTMcGV4QJrSMCuLHWHmy/oDUDvDq2zIJnEIdZUzrKPGNPqwqOqacYbPWRrCGomimEPcNqYmrE4d2Fs1CEGmLFlh5U5CaiucoVKWurBiV8G8lcrQSMMhrNxjCG31sNHzjsPgLGxxXNTiHect0trhfUiJSpmrsMOtTfG6d/QpcFpnvGkrRyAIjLkKO7jnfOsVktYTnrQJZ+kHGeoIa+XmwRt3lroy5z7s4EzsPmwlrtqYs7SVoRolmmft4i1tZeMUVpbAUlvKynwJaww6hv/tq41NSdJVNuoIq7D5Bu1M/F7M5on1EFYA1DZJUVkkYXe2cdirOwBnE8+XjHpgmAgb2p5xEzlJTpnzTMeZ4uCXT8RMZQ9nzsTZRLWSeECcI+xw67R12g5DeOt3xoeLCTlzbxIPzk6t0ypJzpkY4KeiLM1RYuxY6+nTJgkpS1NYQewmUu8sHWVJC5Ok4SwVZWbC9nQz8g67hhC1idQ5k3NVQxszEqYdV7js2oJo2tQdWhJBZiZs+edk911bEynUVM6KIBtSmVmTONedazns2oVITeSys2IKf8DzMkNhA/dhGgZwVv3oMpQyVH2Yggih1nZW/0w2jDGEfZiCmM4mwwYZzj5MQehAayqr3g1gzCzCtCujHHYdgLDtY62sYSy6MgQzHVDCOataxslwysbRh3UJ6Uy+NpboxFVmJOzwXLKjRC3BGsdC2WQgZaPqw7oEclYqa/wpnrIR9mEtwoxD5Fxja+ViLGVjFyYJEWpMWXutaRxjZsIOt9hXavn6FKddRyZAoGVZx1gMZWaDDrHQZg40lpYwjndn2SR6kAGG9dDBfXrCcv/OukHmdefKDE0+HkuESTw7m7SOM7ixtenDWnh1FtfYeowSFXgchExav78EHnqsrTCOL2WT9m9mQY3hWDUVDE/K+GCxoSxkkBlOTW3633Ui+FKWR1I2ztl6EH46MzG8b7aLgZSZNYmz9H9xdsLHCEQaC96VjXq2HoSzM3kK3WwX3fanZq1HiV0clRWTHrWyEMZIWAs3ZeU0VaUsgDES1sFJWTWxWCrzP/RwFcYHkLyLU8xb4RTm1pnVU8GVMtfidPAgTAz6F2+Bp00Xe2eNyfswxkyELY5nR6aaszH2x8Wppwa9X2IYbJU1jcka8GvMQBi/JSL/gUUt7MyRhx7gEXZqqU3ELcxaWfMHsgDGDIRJUbMNzXzH4Va2ke+lcr9Er9gpaxkTleDTmGGEMebPPznqE2clVspaP0JLY/6UmfRhB2fEzNR83DMdGmyUtY35DTJf52EJ3y/RFQtlk3BBBhA2ikU4Vtgoa/7DZ5CRMCPgytqrqYQyL8achcV/MsQwQJV170DLlfkw5iosmfvWh8fRWO7HmKuwBJ4MEQ1gkC0byzwYMzlx1t/Dcq0iLIcG2fJ9uT0Ycx7WJ/ZkiMDAgkxhLHM1Rr+HAXE0lrsaI2FQQEHm3xgJgwMwpny8hJMxEmYBIMjUxhyqhoRZYa5MbcxeGQmzxFiZ+qE7O7Zz4yTMGkdjlvVDwuwxDTKtMZsaImEuOBqzqSIS5oRhkKmfHreTWwQZCXPEJcj4CRm0lkiYKy5BZmGMhLnjEGTCGKiiSJgHHI2BaoqE+cC0WVT8DWqMhPnBzRgJi49ZkKmGHjugE7JIwnY5DvvCgFOQGVd0zAjb3R23N4cgM6/pAZrEEWuzDjJzY4P1YeMMN+uTaGNjQw86xubM9JRsSZmpsaGF5WNrIo1/JesqMzSWgDDJaJpI65+izYwlI6zAm7YsJKuztg0yI2OpCZOYh5tlnTriS/CSMXzC/H7Xh6LbLGpLunQb9d4DiiTMRcBuhUNhYrPcLK58OnRJv7FoEbYjubMDaI+7u3jkKQ5NrawdZL3GBmgSd0r4PzzIcyhiUBTHY6Csz9iQk787LXVteQB7yapTHUO/sh5jaUz+NrUV2NhLzZyy3Bpl9dvV67hTmvzdUXiT1OomHVZlMbw69TetP8hWVLuzMLO7CBwrMahw5qy7VU036BSbtHbWHKcMYE9jrE/ZCmOuwqDXOFfmqlKqyFeFm0TfWOodDmBO05qrT2Aaz/TWVryrMOu7CBxroq/DHm25WTfXjb2IQact2GplWmORIqynzznWDbwOfeGWr4q4laUIr05fpJXKdAfr3IeZ3UWgmslYpY2z0lyPM47ZqHJlb9eXGMqK0mj6MtknAJIYf2yctjvBIdB4kdWmCbr+UGtk17eZ+tvjPehWGdNUMC+W+igjCetUO6sMFm1Ki2U112fVCnUGTSTHPtyqYnrRtqoM2vlrViTVAfoS1nO/xB0AGo9FSMZnNzCwKfHUfl4heogk7Mb29vb3r6u2Yifed5zfv3i92OZpfmNhfjJePrFs/8L29uU8v/XY9g+e6SZzKFyLIg/25qKyjG6UA2kQN7e3v/NN1QeRZjquXdYk58+Smz/vq1wm32aP7+jg7DRfnJiKz289fjXf/+nV21cu5zcebSRjm0CfequlyIPXkvpL5cY8gwvbv/jn7FWPqj6Jcx52+/dXNen5/e7/8q6Ps28232Z25EsswvY22Aczud+bvNTXLt/61fXW158ns3sMnYIij/zaI38MEGGLN99rEWGsJt70LVVh4sx0sAataHS6iJo/+yCrKLnNorjXOg+hKvXjV/d//owIhCpbnxFW5JEHaRIPzz1k0ySymnnfgBHGGxxVlPHnQ9/F+y6mim/z39dk2dt5I5Jlz+HieEvI/ncpv/nt12bZ+6/ORfPLc2Tn6/2+RP13ez8lPI98WVjRmbp0mPPTuj7spuzW1eXjOb9DlSrSTAdH0Y+xdu/vWfbWcw/+RJT54Oxnnjjx0mk++/D2d0XN3WBheesxVpf/e8/n8v+8+n7+N/HV4F2c6mkULUSX1O39lIg88iVhRWfqEs6sBdEI43mxgqnLx7JcvOHl9yuSRRzWLwsr2tGDe77xw9/yN3sb1544+YnNw09+7Mui5vZ/9uvL+xd4sn++kPVhn787L9tKbq3ZaqqQXVK391Mh88iXhBWdqUuHKdsE3SNEWX7q8rGcFyc//QFFkjjCbrIAuv2HpVpbnPpC9qz35vONv37verHN069/2/RfL2bHeDcfqvzpym9kXfIm62/v5oLmGyKlSYSJiuj2fqrNSl8qs+xr4dhhrhjWs9BSl49VxuKN7xwwwlg/9chypS2Ob+b/yLKjT/7yR9eLbf79kpflv8uy537llXflNy7dvvIL2cGxRuvhj74urwcje5l4RNZqeP3z70GPsKITzVXC+GOdzDpMLVph+xdYnWjKd2P74Y+8QpVo0JkOOUScikJLise8sJr72ib77olBR+sTxeBGj2GErUDkahbONijGv62cFQwq7OBDUti1S8UfeMgJuDD5xb/U/gT0ROl9wz5Mi8zVpMO0ZPkMs5WzgmHnEtkokUVZNeAvSskqaP/Hn6qG9Y1P2FkNpN54RfDer3eUqKHINUyEFW2hunxaXwMLY50Dq4iqRShPs9jrsz/L/1AJKz+BfdHNz8PUlLkadZhgin5dWb76hLNLSGFECAIKo0SxE5EwZIlIGLJEJAxZIhKGLBEJQ5aIhCFLRCuf0EHCkEHCkEHCkEHCkEHCkEHCkEHCkEHCkEHCkGEnjP/ALddA1O8AiUCryspFZICc6kSAnPaq3QNyqhMBcqovtgIdk8RO2Gxz+Z15Il6bcrmoCeXFPICcqkSAnPhFGsWm5jnViQA5NS62Ah2TxEpYvbQJsMip3pSvii4vYumlvJgHtJyqTATKKS82hS3cKhIBcqovtgLmJLASJq6P2Gy/AyQql5eaUF3MA8ipTgTJKS9XiUNyqhIBc5KbAnMSWAnjAS2/HfU7QCLIus3qYh5ATnUi0ArRxXG51A2SU5UIthZVrA8H5lQAFjbLsm5L398StxMZfht5os7FPGY51YkAOfHXxqbmxyQTgXJqrdIG9mP2w3r7gQegve9czGOWU50I2IfFOabOqt8Ywnj8Hz5wvvUOkIg3COajxEbrZphTnQiQU92mAY9JJgLkVPuCHZPE+jyMNd1iUGu+yLlOBD4PA+ZUJwLkFO2YilXa4GOS0EwHMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMkgYMtAL06wtm5/WLzqrF5bVax75ukIUjFQYq3+YMGbYd8nCMBph8t6n/OUFH5yKW9cWH8iLew7PfTHLTu+J9aWz+8SKdr7pvZvlxT/eHvQSmNEIm50W6zL5y5FpYx2uuLjn5PnDrQ2mZkP8dXZUhB/fNNssPre69GcIxiKMK+KPDGAvh+emvEdqNInsr/yaA/4f98PVzDbFplvVVUlY2sSxCOMvhSn2wmOtEjbjVzk2hX19Kp7fUYbVTF4FScLi0Bdh/IFzvElsRxgLrSLCis9JWCz6+jD+sjgxbQkrujLRhxWfUx8WC/nEx43GKPFF5+QosfhgnvHBYEvYfeIShMMtMUqUn9MocThEzIDPg5G0iGMTxp8IKC7ggdY/zXQQYSBhyCBhyCBhyCBhyCBhyCBhyCBhyCBhyCBhyCBhyCBhyCBhyCBhyCBhyCBhyPg/eR7h5JrMnzYAAAAASUVORK5CYII=)
m1 <- ifelse(is.null(ncol(Data$data$Zc)), 0, ncol(Data$data$Zc))
m1 <- m1 + Data$data$intercept
if(k_opt == df_beta) {
plot(Data$beta, col = "red", pch = 19,
ylim = range(c(range(Data$beta), range(beta))))
abline(v= seq(from = 0, to = (p*df_beta), by = df_beta ))
abline(h = 0)
points(beta)
if(m1 > 0) points(p*df_beta + 1:m1, tail(Data$beta, m1),
col = "blue", pch = 19)
} else {
plot(beta, ylim = range(c(range(Data$beta), range(beta))) )
abline(v= seq(from = 0, to = (p*k_opt), by = k_opt ))
abline(h = 0, col = "red")
if(m1 > 0) points(p*k_opt + 1:m1, tail(Data$beta, m1),
col = "blue", pch = 19)
}
![plot of chunk cv_cgl](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAe1BMVEUAAAAAADoAAGYAAP8AOmYAOpAAZrY6AAA6ADo6AGY6OmY6OpA6ZrY6kNtmAABmADpmAGZmOgBmOpBmZmZmtv+QOgCQOjqQOmaQZgCQtpCQ29uQ2/+2ZgC2/7a2///bkDrb/7bb/9vb////AAD/tmb/25D//7b//9v////tKSRpAAAACXBIWXMAAAsSAAALEgHS3X78AAALIklEQVR4nO2dD2ObyBXEcdNLLk2d3LV22p7dU5taFt//E1YgxJ/HwL6FBe3oZu4c84B5WvbHsiCzUJQSlYpbF0CKk4CRScDIJGBkEjAyCRiZBIxMAkYmASOTgJFJwMgkYGQSMDIJGJkEjEwCRiYBI5OAkUnAyCRgZBIwMgkYmQSMTAJGJgEjk4CRScDIJGBkEjAyCRiZBIxMAkYmASOTgJFJwMgkYGQSMDIJGJkEjEwCRiYBI5OAkUnAyCRgZBIwMgkYmQSMTAJGJgEjk4CRScDIJGBkEjAyCRiZBIxMAkYmASNTGNjxU1HpT7/vUBopqCCw0/NT/fvtpx8jr7SFVgJ7/+X3wW/oRVOLZu5nyjETCq1WtbBVpROw8SQIrcJ92PvXuqGCPkzAEmUaHAdXA/N4BWyNqfjw4QNcgiRgt8/0oRJaguQG1jvpGJ/PCNga0zbAZr0CtsYkYGyZ+rwEjCgTCq0ELKNMKLQKAmsuw9CFmIAlzoRCq3ALOz0/Br0ClsYEQivHIfH920vIK2BpTCC0cgBzeOcLEnHFxljNdwesqP/zrBlenmM1C1jsTAFzSsASZ0Kh1Q7A1Ie5TSC0Sg+sxbPfdgqYTxBY0YYCFm0CodX2wIqJv6cK2HgShFabAyvK7bdTwHyCwGwfJmC5AzOTAkYGTH0YG7ComYzVzA+saVICFm0CodUGwIqy+5kpnYCNJ0FoJWAZZUKhlYBllAmFVhsAUx+22ARCqy2ArZjJWM0Ctospx0wotNoOWO96WcCcJhBaBYFVgyCqexPH4/nmgfX/zixgThMIrTzA6oErx7/OeAUsjQmEVh5gxy8/Ysc4V8dDAYs2gdAqDOzrw2//qFrYl4gxztW/6sPiTSC0CgKr7tUuPpZvUWOcL8DmSydg40kQWjmAObwClsYEQqttgDV9GPgQAZszgdDKDUxjnLfPhEKrjVrY0pmM1Sxgu5hyzIRCqzCw6ae5CVjiTCi0CgJL86wp7/3bjNWcGbAkT3Mr2lDA5kwgtNqnhQmYzwRCqyCwJE9zEzCfCYRWYWAeb7B06sNcJhBa7QTMO5OxmgVsF1OOmVBoJWAZZUKhlYBllAmFVgKWUSYUWglYRplQaCVgGWVCoZWAZZQJhVYCllEmFFoJWEaZUGglYBllQqGVgGWUCYVWApZRJhRaCdgOmdx/qwChlYBtn6low3TADlN/pHR4wwURsNTA3r+9vH0sDx/n18VeR0EELDmwehAYuM/G4XUUZMl2bv0Y0z0zpe/DTt9fzv8fP+cDbPMHBeeYCYVWzeIzq7eimHoFxLx3YekEbDwJQqvA4vJWd/6uBrbfeMJbALv0XrAPu9VbZmEfhjs2lHTHEbv7A2tfYASeFJDuPc5F72xpusizyyeaHQEw38OBUGjVb2FQqVpY8xOzne0YtGhgu46Jd2TqVcCsCYRWgcVlsvc4TwEDQwObtnip8/5yL7DePh1YM3J5YCeZvBRJD+z0XPz0n6nXTgW84YIMgBVG1+bQxWWz9WWzrBt+a5zjWV2L7MXjz3SpK75z/ebH+IvrJplMqMpAiCv99Px4/PIDHPUc3qkpO7PotZve5vWD66/rzwhYs6AzjeJe3bRr25XsDyxO0ZZ6ULaRaSZTs+3X7bgUrNtDJwh5+7AzsNlvOsAY5/9JW8gD7NLCDpu2sCvsYQtrd+Hejln0zyd7xw/TWGC76Pbe3t48bsuuFgYzTLWw0YZ0ZR00pvFWLu3D8Fm9wzvJpiy7UnUzC6Oyd+hqtvIS2Q8adUe9kxI7syzLuZW86kB4Vh4Vrz2UthsyWP26rVO1Gqj0Ka35pmO8D02s2W5Vx3fQLLsaGZq6k4zrmh2wnjdYULC8j3p+zalM11Y6kX4jYKuuw9zAuskBul6mmY0fJEUUIz6+V5Dr8WF5pqKPu+idO13Stw2u7CvikAj/HLbqm44FwKYyTRwnPUkXAxu19RV73vXQ367ZABxyHOZAuiy+vKsZnnSs+6YDHMgm1gwsHx8+tgZWgra+YkO62hjOXAZs5svfVN90+GZOLMd92G4fnyKTYbMSWPlatzDdIrBhpmJwht/NXNCHtd/W6yacTTMBYGYShFaBxU6vgLmWhy4QUGglYBllQqGVgGWUCYVWApZRJhRaCVhGmVBoJWAZZUKhlYBllAmFVgKWUSYUWglYRplQaCVgGWVCoVUiYBN/xRp/iIDNmUBola6FFdMFETCnCYRWApZRJhRaCVhGmVBolQiY+rAkJhBapWthMwURMKcJhFYCllEmFFoFgS19y6yALTCB0MoDbNFbZgVsgQmEVh5gC94yu3QmYzXnBmzJW2bxTHPvKywEYTVnBmzRW2bhzOpfeJs7eTVnB8zhFbA0JhBaCVhGmVBo5QaW4KRDfRhXC/PMZKzm/IH1xxx2WWIHqgjYeBKEVmFgzhGYRffLVzoBG0+C0CoIzDs+TMASmEBoFQTmHYEpYAlMILRK1sLUhyUwgdAqCCyLEZgbmHLMhEKrMDCPt5maGjUuYE4TCK1SAtMhMR9gB/Q4YAFLnAmFVgKWUSYUWqUEpj4sH2AB7423U8B8ErDEmVBoJWAZZUKh1c7Agk+qIKzmewYWfhYMYTVTAws8TFLA5k0gtEoLbPiA2/EICQGbN4HQaktg7YLeiurD5kwgtNoZ2DbbKWA+gY8f3ehhP0TA5kwgtEoMbFCQBaP8GKv5foAtmMlYzQK2iynHTCi0ErCMMqHQSsAyyoRCqzCw27y0lKiaMwN2q5eW8lRzZsCSvbTUNZOxmjMDphZGBkw3krIB83gFLI0JhFYCllEmFFq5gemkY/tMKLRSC8soEwqtlgDDQ2ZjSydg40kQWqmFZZQJhVYCllEmFFoFgc28DU7AEmdCoVW4hV1eaDrvFbA0JhBaOQ6J799eQl4BS2MCoZUD2LRX2kLbAVuXZdEny+RcGw7oW/qZyz0yCRibScDITOrDyEwCRmYSMDJTGmDSbhIwMgkYmQSMTAJGJgEjk4CRScDIJGBkEjAyJQD2/hW9IHNG1RDBp3hfPY4mznR6Lh5eYk3n4lU3HMWYjp87g9tXm6LrYj2wqh4PHyMM1T0ix59fon2H85ZFml6fqmFScaaqeIc401tFuDG4fbUpvi7WA6tuuq93Fq/eqpK9PsX6jn/59Snyw5rxAHGm6t2RzZtanabXh3+dV2wMXt/FFF8X64HVGzh1Y9WUzoZI3+n7b+fdMM50/PLP6pAYZ2paWJSpquzG4PddCcXVxXpg1dDMWGDVvY6RvsNjddyIMx0/1YgjP+nSn0SZqrpvDH5fAyyyLm7Swt6/Psb6zmufFrSwyJ2+Nv38Ur2ida8WFlsXN+jD6h0/1neob9l7jOzD/lbXQ5yp2dujTMf4Pqw9S4ysixRniY9xZ3uXMsb7qhYWaXp9ujTNGFPTwqJMVWU3Br+vbpbRdXGD67BLY3na4zrsvHbsJVV1vh198bb4Oiy+LvRNB5kEjEwCRiYBI5OAkUnAyCRgZBIwMgkYmQSMTAJGJgEjk4CRScDIJGBkEjAyCRiZBIxMAkYmASPTvQIb3jR2fR/JHUjAyHTHwKo765uhPH/+9am55e018mbI7HTPwD5dblo/M3o7g7ugev/l3+NXXDDpnoF9bu+dPh8Sq1vXq8lDMf8kwdx1/8Dq3uz1qX6c+0M9fO7WZVul+wd2bWHNkfD179Rd2B8AWK8Pe6vG6f33O3UT+wMAOz1fzxIfXq6DAnl1r8DuVgJGJgEjk4CRScDIJGBkEjAyCRiZBIxMAkYmASOTgJFJwMgkYGQSMDIJGJkEjEwCRqb/A5qDVuPXwyRrAAAAAElFTkSuQmCC)
beta_C <- matrix(beta[1:(p*k_opt)], byrow = TRUE, nrow = p)
## satisfies zero-sum constraints
cat("colSums:", colSums(beta_C))
## colSums: 4.167479e-09 -5.502439e-09 -9.87453e-09 -5.177147e-09
Nonzero <- (1:p)[apply(beta_C, 1, function(x) max(abs(x)) >0)]
cat("selected groups:", Nonzero)
## selected groups: 1 2 3 4 6 7 8 9 10 12 14 15 21 25 27 28 30
sseq <- Data$basis.info[, 1]
beta_curve_true <- Data$basis.info[, -1] %*% t(beta_C_true)
Nonzero_true <- (1:p)[apply(beta_C_true, 1, function(x) max(abs(x)) >0)]
matplot(sseq, beta_curve_true, type = "l", ylim = range(beta_curve_true),
ylab = "True coeffcients curves", xlab = "TIME")
abline(a = 0, b = 0, col = "grey", lwd = 2)
text(0, beta_curve_true[1, Nonzero_true], labels = Nonzero_true)
![plot of chunk cv_cgl](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAqFBMVEUAAAAAADoAAGYAAP8AOpAAUgAAZrYAcwAAzQA6AAA6ADo6AGY6OmY6OpA6ZpA6ZrY6kNtmAABmADpmAGZmOgBmOjpmOpBmZjpmZmZmZrZmtv+QAACQOgCQOjqQOmaQZgCQZpCQkGaQtpCQ29uQ2/+2ZgC2Zjq2tma225C2/7a2/9u2//++vr7bkDrbkGbb25Db/7bb////AAD/tmb/25D//7b//9v///9oE1NEAAAACXBIWXMAAAsSAAALEgHS3X78AAAKl0lEQVR4nO2dDXfTNhSGU9KxlsJoB/vqWrbRdoN12Shxkv//z2Ylcew0dixdffi+8vscoD2HyjZ+uNK1ciVPVgSKydAXQNygMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDo19YcT4xvHhMcDWkl15hy5vr9df56VP0iyH99ApbvHvc+9psS2LgKexIhHH8i4GvsNXiau29ZQyjsBh4C4vTlnRBYWCEEtZIOiyHR9LLWflr+3sHI0wXZw1JZ20/QGE6OIikLihsYGxFVVDYcHR2e8fwFbZ9DGt7EKOwLlyjqol3hC1vLsVtx4gkqpr4d4mL7+/FbceHlywDx7B0iLvBJhSWhiCyDBQWn2CyDHbCZqdPs8nkOuihR0JIWQYrYWVeUf4qLtzKAHaHNlUEjrKzIbQvS2HvHssYkwozaWTxqiOVzJmgXWGFZZc4ObmfS7vE+cvyj4exhVgMWYZEScfXq1+/fjVfDUe/epxQE7F8JRLWPR1yQAZCo3SFFXbCljeT0y9dMxr9h15cWfvqpV2kGqLKMlgJKwOkePPkWnrYyBKjD2A6BMaWZbDNEkthLaWHVodO4OuAQQSm8OUSYTNhhM3WH78MmyVGFxi9K6ywHsMmE9dibM0zHWEFppJl4FziGq8sNKUvCnuOo7hkXWGFXdJxJVmbgimsQuvzn2WEmfnbl2EPDUaLwOTRZXDoEueOi/byEvacoSKQESbjILpSTalxDEtEKIG2Mx1dP2D+ytQmthjNVliQsasvIrueE+1mOm47531LYWudxduOtiQsll1i1yrL7TRj+xrnIBeojEEywyb+pdondx9MhL3hGuck+M90LG/K/LEt5d/7ADOHCoHBo8vg2yXaHXo22qqp4DhE2MzxY+O6bfH6PbowFdFlcBB29APMY2ucl7d3k+lqNZ2Kr5LsSDE1NbvcjmHTKaQ3NdFlcBnDhHWJZQL5POmYTqd43nTgnSV27+a2XyLQNgACSFMVXQZfYXZ7TR1J6xlsbtgJM2Xas5O2+alju7nV3/Y/h6mUpi68VvarV8o/WxdDhNzNTaU0bViWuRkr7YWkQXdzU9Q7aowug12XuLbiukusdC5RjzSNKK2aGtaZ1ugyKBXGMOtCq7DBlGmOLoNeYQbG2QG6haVGeXQZ9G/7kDDKAHwl2fbBF0VPZ8OTYNsHDLQnGxUJtn0IA6NsA1TSEc0ZSHQZHCp/pWucAxIpzIB82Qjb7RIr3UUgJKPvGX1r6+WHVgJKslEBNYZtGXWU2X7i7FdIqhWw6DJYVk1J6kCjCgsSZYC+gMewkXaMdl3ig2RzrwRdoocztGSjIs1iiGiML8oQs8R9JM5Ao8uQoPI3Nu6DGbAvS2HdG1yqeMvsmPIPy7rEzg0u7Sp/FYGabFTYpvVdG1yqiDDDWKLMIcLaN7iEeo8zeHQZrMcwhA0u+4IsA19R0/p/SAQCCWtZ4zz0vyxT+oUt3n1WO9PRSmvPmEN3aMCf6Wgh54Qxg5mOVvacoT97NfFdMqvmOew5uT6X+S6ZVTzTsVWWUXQZfJfMqo2wDdPcfPkvmVU+05Fdx2hTlygqENAgbJNs5OXMSthnjTUd9mSVf9h0iQ/bqQuUB+cNzVw+I2VjiLCsyFRYS26YSZTl2SV25PI5KMMtJJWQQZTlN/l7fOIQXplv1ZT80MMBLc23akp+6CjkNhF1iG/VlPzQMbD3BRtl3lVT4kMPDejuHxlVTQm6Q0Bn+WSJ+Q9fa/IRJgUsyvLoEv2KNqBGM9ukY7Xe0i3koYkIrJ1wWgkzeKFEGX6EjSTZqPAdw7J6yyzCYOabJQ77ltnwFaLqlQUQlt1bZlVL8638HfAts1EHL7XSfCt/7d4yG4PYyYbSAc238ld+aC9SpYb6pHlX/lod2nV/5z6SpvK6pIWaSzyadMydZR8j/YOXpt4xxeTvw8nHgBE21IOyEmkRhB28x7nsEifWHD/2wNMaCl7YmWYFpsMYdtTmWb/SFAwpzVeY3fqwIElHFV1uMRmLoYLNt8zNbgVmCGGdvaFjrxqW5O9+9y1zSxRhTpOG6e1N03nzLnOzWoEZ+jnMgbRht+0oI7oDKHMLMyWfvMNsuAs53uVR0+HCcLlKw53rt43LP34Oy3/WrG3X7RC3JOKWKIPmmHKUC4sPmjW7pKNvc7AowhJuOAQUbA4R1mpFfujjDDEHBWHNQVjKMrfh5gy1B5uDsLZPlT0O3YmGvdf0anMZw1K8tFSBrBqN1pQthlDla4MyaQ6l2oEPfYCGrrADRdLspqZuHdej2xx6D8WyKpRIC/McJjn0DgBZFQqkpRjD5t2ugWRVDCzNQpjvfonms5XZy8O/B5RVMaC0BMIMuw/EzjaigGVVDCQtkbBvv9mKEh1JLQM8qNkIq15tL046ivPWhRSZkHZWJFGELVw3qoIjlbZEwlYPkld34xHfWgJhpqAq/whrEFVaisrf2WSS8xjWSjRpaSp/x0kUab7CFL97RQXBBzVGWAJCWvMew5S/e0UPgfJ+ZR9gZo93NSSFDYK8/DiUMCYdItwr/hlhOrAWF0FY8nUiOTFppfkDPe19zu3RlnRBYWB4z3R0f1hGYTHwjrDNbqWytsQd/y6x85MTCotBzDGMxCCesJDH8Ww/ntNb/mjvWr/x3LGBT09hYKenMLDTcwwDOz2FgZ2ewsBOz4dfMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAxfYYuram/n+jtZe7Nw0HnZbfOk1TobWfvljWBR4t7Vi17wVG2XYX33PIWZm7TZU6X+TtbeVI4Urxxv2d5JZ+6+G+0frp1fePfs6l13kTdU7/Gyv3uewkzF/eY/Sf2drP3cXK7ryvbmSYvX790DdO/ynambmzeECg6xe4+X/d3zFLa+0M0LNN9IVqbvt3Ju32i+vL1z7xKbl/+be5dYN5dGWGXJ/u55Cqs3EJBtJbDXqrsC0qL57FIwhtXti/PrthfpWp9eNoLvhNnfPUURtrhy9bV/eoEwz8tvNC9HX9ctkzeHSBxh4caw9f9x+eln65I+V+ONy/9RIKwxAku3KikSj2GmF6uyxEtRlli1kvjaP6kgwhrtHwRdYuPqPSPM/u6FeQ4zp/V5Divbb0LE9ZbXp/d5DttevvsNr5vPhXvLmLZOd48zHWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBkauw5c26qOf0v4vH4nxdBf7icfOeJvD3MOUqbLWtICv/KF5997Ra/FB+41w2qZAxCLv46X5V/Exh2qmF/XW9+vcjhWmnFvbp7fLDp4vtGCZas6CHUQj7+/cvvxSMMO3Uwh7//OOSwtTTEDY/uacw9TSEbX+dixZc6CJjYXlCYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYPwP23wkq0qTIN4AAAAASUVORK5CYII=)
beta_curve <- splines::bs(sseq, df = k_opt, intercept = TRUE) %*% t(beta_C)
matplot(sseq, beta_curve, type = "l", ylim = range(beta_curve_true),
ylab = "Estimated coefficient curves", xlab = "TIME")
abline(a = 0, b = 0, col = "grey", lwd = 2)
text(0, beta_curve[1, Nonzero], labels = Nonzero)
![plot of chunk cv_cgl](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAABd1BMVEUAAAAAAAMAAAQAAAYAABEAABQAABwAACEAACgAACoAADoAAGYAALYAANsAAP8ABAAAEwgAHwAAJJAAOpAAQQAAUgAAWGYAZrYAcwAAzQADAAADAAoDABcEAAAGAAAGABcIAAAOZ74QAAATCwAUADIXAAAXADoYNgAbAAAdAAAdgf8fYL4gAAAkAGYnAP8oAAAoABwogf8rACssbzo0WAA6AAA6ADo6AGY6OmY6OpA6ZpA6ZrY6kNs6nP8+OjpJAABJAP9Mh75mAABmAA1mADpmAElmAGZmOgBmOjpmOpBmVQhmZjpmZmZmZrZmcy1mgWZmtv9m2/9ro76B//+HQgCQAACQOgCQOjqQOmaQZgCQZpCQkGaQnWaQtpCQvD6QvHuQ29uQ2/+2ZgC2Zjq2tma225C2/7a2/9u2//++YAu+h0y+o2u+vr7bWA3bkDrbkGbb25Db/7bb////AAD/AP//Ugb/gR3/nSr/tmb/25D//7b//9v///86snbEAAAACXBIWXMAAAsSAAALEgHS3X78AAANJUlEQVR4nO2dCZ/bRhmHFa7SLRjKAmHLXe77dLah5lpCWijHcm3KEZo1FAwEQ9cBLNsfHo0OS7I10jun5pX+z2/XVhKPJOvJOzMazRHtACuivk8AqAFhzIAwZkAYMyCMGRDGDAhjBoQxA8KYAWHMgDBmQBgzIIwZEMYMCGMGhDEDwpgBYcyAMGZAGDMgjBkQxgwIYwaEMQPCmAFhzIAwZkAYMyCMGRDGDAhjBoQxA8KYAWHMgDBmQBgzIIwZEMYMCGMGhDEDwpgBYcyAMGZAGDMgjBkQxgwIYwaEMQPCmAFhzIAwZkAYMyCMGRDGDAhjBoQxo1tYPI0Ezyw8nA3opFPY9uoyfV+f3jg/GdBNp7DN/UXtvZoWuMBQWEuEofxzgamw3eYi9d5QhkGYC2jClqc3yyi6tLproAVJ2ObedfITn6vVBCHMBTRh9xdJjLUKq1Q6iMUj0IKYJUYn12tkiSFgXOlwkxbIgDBm0IRtr6LTN+5dW9010IIkbHs1i+/eNLY+5bdhTTdiEOYCai0xEdbQ+rRLZertGmihEGHL5vbdjSyrhDAXkMuwKFJtj4cwF6CWyAwIC5iz5Cf/3UOrdFzoPJ6EME0KSWdN/0iMMNER4Lbicfe7FokVm7VGyVE0NaGQJa4V+20UaUU1Mr6jeNc9NjpFFXiIsLVIOEeINUOJqiqeyrDVxW9WK/G+Ut/PYJGWU21QWzo0zqey60pzyCoXtxqxO9WoqkJr6XigUwKVu95czCSfGZ84naiqQswSZR1tSLuOp50FWCFu2PZMZQk83DgTfFWoZ5nDEWhDlsCDsGUanrq1xCGUebZkCXxkiRZhJ06rJtiGQoQtZ1Z3bUSZZQYsULsm2IaCMNXKvc+2xOAEOpEl8NA01Qe9CnQmS6BShjHul+hNoFNZgpE+D3Mi0LkswUiFHWIu0IcsAU2Y6Ka9PBlRv8RDgZ0Wffkij15JXhsHQ4iqoyjiGprzOQs7pkWgl6ywgNjNTdQ3mjuS3l+k1f34RdVdMycXd3Z2IFL2Tt5rY/oSWpaYVhMba/VZH9PmMc60c2SMSmR1CaWWnaaVjs3FyaNXRYTdHd8YZ5854R7zWuL2KrrdeFM9aGFey60qqNZr0JcsAYSp06cvlT4drY2/Yxnj3FtWWEAQth8Cpth1anjCere189VragCEIEtgXIbJZ3OrdXNDP1JLUNsSZV0EaHNNLXn3rQ8lugTE52HS6902m9t+K37hJbbCQpIlMC3DKBG2ffAomuidXs8EJktAyxLnM+knCLO5LWep1gk7aQH68tHNLb57k8fhhJGz0LLCAm8dSWfFXzOSFiJ+mqYOq/WBh1qo0SUgT/tgNHVRw32YUBaktYBlCRQmVlGdOJsQnQEaC9yX+dRF2rvOCSlzDDkrLFCIMMnURbq73hN4gRYY5DLM8dRFvUvjEF2CkB5g9i6NAyEJy3NH/9q4RJeA8gDz/l/8DujzK42RLEFYEVbiLXtk5itYYV7urDllhQXOB0PkM4MlbzpLjKEacojpYIiutPnMYOJN9T4ux1HmyDG6BKaDIbrS5jODiZ7cel15di7urJnKEpgOhqDsOgmvdyT56YnRIn42pTH25aPSkc4Mlig/nZjWJGxI45oVFrgXls4MJia43A+YKMRN1BWg3dH5jXM205Qo/2QT3CvL05bGPboEJGF6lYV01/nMYLUIa4EsTl0aS1tilF99pB8twlqUdfT8LWYGW0eRyn0cIeqGlzmuSkHHogooZdhcuh6Ol1VmJ/IyjyotyOg6lGNtyGxrhJF6/tqlJpBRpHVHT1MajUHprcIoEabZLEUiNSbqjlJtvUVXl6DV3sdKsjnJqaYyzRI7e/6KMu7zu+W7k8+8/VfEr6rBJH+W5vHRzFGWVv/t9NHwJfLvcKypxPX4sM29P2S2n/OwOkT+Jcss80znXk/QIWO3kgVFm4+Dc63LkSuq4/LG+XHC3599SxR95vt/iz7++B9f/+VjYAxJWHdH0oYxztn+fx2VfKHfrzoMKMIMOpJufpCqunXr+eTly08aST5We1+1vZ+dZWvKZGR/TrdW+639v67KrdXxVvG5LEM62t+Tw89Jjkv8XIWGzzVS+dxRCLULM+hI+r/o6aTW8ek/fuQbUfSurzXcHGZZfv7bWmeSbx4dWnzFSpnQeln2CGnFb2CoCmvrSNrR0rF+//THUfSe94rPfOxzq/zw+fuq+POu9d3C9+06Tr3wr2+11dn8Y9qRtOs+LJ7+XLh6U/TcV59/8/eeiEwig7y1K3NOe1uUCnSNSWuLi0dMa4ldLR1ZW+K3/yNeP/QLrVNs+99P2Wp6KGOuvS+BpsI6Wzri6VdufTi6lXp7p9ZFV6CtUcPyNT2WOJnQ/k6+RTmucd/6jpaOePqpn33wm5f/fet3XvviR/UijEp3E1S/JdFhRB79tsR9CbXSkbxqjV75Z3YH9rZPfPK1zz7tdFlFUpNhQLWHY1pklig0TWlV6/PnYZsffjd66kuuWoAVu9WELK0T1xG2y/ozij0sVRfRpKPeIM9WmofxYWkHVNFryqSXWwuaj0+Czh3lhNu3norJ4y6GzpgLs/Bwkpk054MhnGLnYTKrzNH1YAh3WH70z8WZ68EQfGBizMdgCPs46ljDIW9kXumwTrAzKhX4EaZa/LXhodtayL0dvQhbu+uX6I5A80cfwuYnDy1FmPdOoeFJ85Qlcl4nIixplK7axcoQ+hOriDLMdHmPHgc0hJQ70iJsPkteVFvbGysdmtoCGH4SiDPXz8MyDmuJLFfTCaPC7+F52E5aractgxRAdO3pv8ObwvOwxhyRtspsx30YTVw49CnNQjc3e6vMHosLKbrq9CXNgjD7q8xGLEKun9yRnCVKRq84XGU2ia7wxXmXptC3XvJ4xccqs2GHnFdprkevaCIruwIV5zF3NB69or1rI0IU50earW5uFisdKjXD0MS5D7WBPMAMqJBz/AjUQdOU0cUzve8KRJw7Yyqt9TojMBWxdp/cf8g5yhwVIqwRm3NNOWnV6FGckwLN9QhMBdy2QkVVnB6pjm1pNGFL6QNMWxHmt83Qszyb0qg9f9e3JQ8wCavMEuixjdeLuIm1uiO1DCta5e3tuiSMFnkPeaYNabSWjgfXyY+rvvVh+Krjzp7pI1BaGZa4WldW9j1i2fRPpC8aoq1DHIjTl2anpUNbGAdfeyxHnd4EO9TxYbJaYoaeMFa2DrCXZypKUxgfpkrX9+Dsq46xOAVppi0dmrsOt6+GCQZhR673E2+cZ2qHp+x66OjIIwybMW3p0Nj1MKOrFSVx7cp6LMNGCDnPlM8Q57kMG2F0SYkOOfpEU8Hmowxbq4+PHiUSe3VjtAgzWsdZNGmJhmNElwqysPPUp8PmIOdxcWjNk7APvM9gR6ACpU+H8dL28VR12iMgw1OE6d0YgGMcdHNr3PXc+Uo5I8G4m1v3rtsWLAWq+LhxXqqtfwnaGEhX7fFAbOk4vVlKFmyz2/MXdEFt/E1+Gm9+PawyC6pQy7AkxhqF9bDK7LihPg87uV43ZomIMM8YVzrs9PwFVFBLZAapLXGRzurgc1A6kGFLGCodnkCEMcOBMJvdmcEhiDBmqLTWW+qXCEwwHuMstwlhLjC+D8tmK9VLC9Qxv3GWPpuEMBe4bOkALnAnzOZ+DNOP5/AmQ2Y1D+ki/XgOD2HMDg9hzA6PMozZ4SGM2eEhjNnhcfPLDAhjBoQxA8KYAWHMgDBmQBgzIIwZEMYMCGOGqbDNRTEyutzSSy8GDioPbK8etBhno5d+e6Ux7Ld29lqzNxVDvMhXz1CYuEjZdPflll560XMkvqN4yWoHXar7rqSfXzaNqKImF2evuvyaYJ1rpl89Q2FpH9TzRW1LL/1anK7q3BHVg8YvvKQeoLXTV6ZMrtMPN2F+8jA7e/rVMxSWnmjaq6rc0ksvUE5fSb598Eg9S6ye/m/Vs8QyuW6EFZboV89QWDlFh95kHbVU8h6QhOTLmUYZVqaPp5dNC+mSD69Xgu+F0a9eQBG2uZgZHV5DmOHpV5InpW/TQrzdu/AcYfbKsPT/uP7hs0mLZ9rpNz/SEFYpgXUnA4o9l2EiFytqiTOtWmKRSsdX/aAaEVZJP9fIEitnbxhh9Ktn5z5MHNbkPixJn4WI6iUvD29yH5afvvoFL5OvNWdvEmmVrh5aOpgBYcyAMGZAGDMgjBkQxgwIYwaEMQPCmAFhzIAwZkAYMyCMGRDGDAhjBoQxA8KYAWHMgDBmQBgzhipse5V26jn99/kinqa9wJ9ZZOs0MV/pbKjCdnkPsuQlvvOtm93m5WRjCEsXj0HY+SvXu/gnEBY6pbA/X+7+9RDCQqcU9vqL21dfP8/LMK0xC+EwCmF//d0bP40RYaFTClv86fczCAueirD1yTWEBU9FWP4z1RpwERYDFjZMIIwZEMYMCGMGhDEDwpgBYcyAMGZAGDMgjBkQxgwIYwaEMQPCmAFhzIAwZkAYMyCMGf8HbMPwaqqvhIcAAAAASUVORK5CYII=)
## set a thresholding for variable selection via cross-validation model
## example: cut by average L2-norm for estiamted coefficient curves
Curve_L2 <- colSums(beta_curve^2)
Curve_L2 <- Curve_L2 - colSums(beta_curve[c(1, nrow(beta_curve)), ]^2) / 2
Curve_L2 <- Curve_L2 * (Data$basis.info[2,1] - Data$basis.info[1,1])
Curve_L2 <- sqrt(Curve_L2)
plot(Curve_L2, xlab = "Component index", ylab = "L2-norm for coefficient curves")
![plot of chunk cv_cgl](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAeFBMVEUAAAAAADoAAGYAOjoAOmYAOpAAZrY6AAA6ADo6AGY6OpA6kNtmAABmADpmAGZmOgBmZjpmZmZmtv+QOgCQOjqQOmaQZgCQkDqQkGaQtpCQ2/+2ZgC2/7a2///bkDrb25Db/7bb/9vb////tmb/25D//7b//9v///8cX6YOAAAACXBIWXMAAAsSAAALEgHS3X78AAAJK0lEQVR4nO2di3biNhRFnWkT2oRJW5K+YIa2gcD//2Et2xASbKx7JRufrL1XOxPGlmTY0cPylSj2IEVx7QsAGwgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxGiErW9f1kWxuO61QAS1sNevy/K/7c/fr3w10Esj7PF7WceswgoYghhh+3Vxs9xYm0T6vyGIEzZE1uBidGG9lRouEids91Tc/vt1mSHror9MuESUsN3Tw/b+ZXP7kp41whKJElaOEkth5Z/pWSMsEUMNW+eoYfRhiUT3YUVh9EVFGgSG9WIkC9vOFqH+tVQ/hA1B3KBj3tkc7p4W+9Wi9HZ/dgbChiCyhm1nRXHXdkI5dNw9L6u/jVmDC0OTuPnSNqwvq9fmoTx4rhNhQ5Baw0pj1Rxyy0GEDUFqH+bPGlzEznTkzxpcxM10PPfO+544jXzUBi4im8RKQOugw501uGCmQ4wMMx1dtQ9hQ5DaJIaZjkDLwzKEDYGhhq0fWk44DDYuz3QwAsmGQVjr4D6qhvHUMh/JU1PdzSXChsDSh7njEhGWj3GG9fRh2bjmfRgaHcQJC2Ha65sccYkfjmPMSlwfVsWQmhdDRBxHmJXIMLcw3sgSSPrhOMKsxDWJ1TDROPdLHzYITP6KgTAxECYGwsQwxHRkWb0CiUQIa2YSzashEDYERE2JMUAfRtTUkMTOJRI1NREin4d5di1C2BDQh4kR1ySu2uJvErMGF0T+isFMhxgIEyNOWL6tiyCRyCfO2bYugkRih/W5ti6CRAw1LM/WRZBGdB/G1kXTgFGiGAgTI+YB5uO3vpmO1hhThA1Bag07Po4+t4mwIUhuEl/npSpq2GhkWAzxOr/9B2FjkWUxxHbW1r0hbAiuuBgCPFxzMQQ4yHUfxl5TIzGRG2f0xpLnxtmT9cdTMRZFlLBLorLsNYWwaOJqWLeyPHtNISyamD5s1Tn7FLvXVP9l0IdFMo0aBtGkCovaawrykdok+rMGF8TWizGRG2eIhUBSMQgkFYNAUjEIJBWDQFIxGCWKgTAxaBLFiB10lH8y6JgChqkphvVTgBomxgB9GFFTQ8IoUQyEiZFhMYQza3CRZTGEK2twwWIIMVgMIQaDDjEQJgbCxECYGMQlihE3rH823jPHZA0u4moYe/5OBvowMZKFbWdFmGRMXB8GsaQ+DwuzVuH5JsJGIvWJcy1qdYewkUiN6WhWYK5/OJ/KR9gQJMd0vM6rr/lYs/3eOBCXKEbyPh3urMFF8qL0BvaaGgkWpYvB5K8YOWY62KdjRFKFsRPOyKQKy7TXFMQSG0jaNeigho2MIZC0HfaaGhdGiWLENYlrvhZ4KkQK48Z5KiT3Ye6swQV9mBj0YWIQ5iYGYW5i0IeJQai2GPRhYtCHiYEwMQhzEyM5kLQlDVFTA5Iaqu3PGlwMUMMis4ZzIpol+rAJUez7PzVGiRMCYWLkE9YdNeXPGs5p78Pe/Wvk1NTCU7ojDZzzvt4xWz95HML2q+k/cf68N+qeGjap2fpWNTH9tSr2PsxZTkLanmzPs/7Mwt6BMDGihdkDp0YVZunDsnZ3Y/edesKSP6KslTFDZrb3kyxMbwXmxIQZc0jtwwTXh31+Yd0btl17BaavdZxWH5Zf2IUN265cw/rfa9/HOYXb7ex92IUN2/bDrsA8eSu+u+W+E0a+Gcjx6xFZw/btG7alZN2Z7PiuTj5O081Xaw6tn1bHccuPreW2/9h3NV2ZnV3ypWThj+4N2/rTmmm3ZLn5MuTQe6rlatJzaM/s/EBPuiiy7TXlfIPtOXT8ep+c21OdpYWNEzXV30pZcuj719bjQwnr/fWJudzJCcvQM0c8t718fKA+zHI1SX1YQ7swvZkObVKFCc50aJMqjL2mRiZCWHNrbN9rCoYgyef+YviAd9QxZjKFa7QkS2rWJvuurlUYwtKTKVwjwq5VGMLSkylcI8KuVRjC0pMpXONowmB8ECYGwsRAmBgIEwNhYiBMDISJgTAxECaGX9jr3LzXUWDd8TD0MtsQlGwusUplLTAEHS3shTXJrKVt6tPjS3MLC6ED6zt7upVnR5BNeFfmEqtU1gLDl2Rsf1paC2uSWUsLv1JlMYbS3MJCUM7WGIy/d24Hvbr5uyzJWmKdylrgJnxuq4W1sCaZ5+2VxRhKcwvb3r94vrKlihCxV7LwbuwlhlSeAstSPG+vPN9TWlm1DKW5hYUoKoew0Gw4fg3DR28vsdJsLzCslXO8vZDMXtp2drO0lDZ2Dauw92P+GmYvsFrLYy+sWQJkLc1an8fuwyp8wuwluoRtZ+Fkc2F1MnNp9flj9GGhAXCMEkPl3/3uGtbbSzw0pJYCmw/eWliTzFpa0xYaSrvGfdiNaySVcB9mKbD+rsKFtbBDMuvba84f4T4MrgPCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTQ1VY+JZVTwxQoA4oex9WdthlcPKICqt2bF85jbVFACJsUOoo2fAhv86rhS3PfxXFQ1nrwveR/HkIHPsSooX/OCweOr6qfw7aTo7++Mshr5Ur3nI0NIW9bY+6qsLgd093++3srloIMr992dSfe3lgO6vD5GsLzatjk3hydFOKq096ffzm+s7dkRAVdqgCIcY5RM4+12sQqvUji/Bjc6CKyS8tljWy/Jfm1Zuw5uhjtfisOWm/LjxfuTsWmsLC4oH6h3oF2Kmw0FquFs2Bg5KwCuhm2SGser1aNCdVK1Cmi6awYx/WUsO+LltqWN3I9dSwpiVc/TbhLkxU2NsosenDTpvEu2O/VHdX4f93rz4KO+nDwkn3/3lWiY6FqLDjfdhhlHhSw349HSU2SspXobGrX+2eDqPE5uju6TBKvKkWN7ds+T4ZVIV14l1kqALCxPh0wj47CBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLE+B8Icpl2wQlBewAAAABJRU5ErkJggg==)
cutoff <- sum(Curve_L2) / p
Nonzero_cut <- (1:p)[which(Curve_L2 >= cutoff)]
cat("selected groups after thresholding cut-off:", Nonzero_cut)
## selected groups after thresholding cut-off: 1 2 3 4
y_hat <- predict(cv_cgl, Data$data$Comp, Data$data$Zc, s = "lam.min")
MSE <- sum((drop(Data$data$y) - y_hat)^2) / n_train
y_hat <- predict(cv_cgl, Test$data$Comp, Test$data$Zc, s = "lam.min")
PRE <- sum((drop(Test$data$y) - y_hat)^2) / n_test
cgl_result <- list(cv.result = cv_cgl, beta = beta,
Nonzero = c("Original" = Nonzero, "Cut" = Nonzero_cut),
MSE = MSE, PRE = PRE)
ignoring the fact of one-sum for compositional data: naive
## set mu_raio = 0 to identifying without linear constraints,
## no outer_loop for Lagrange augmented multiplier
cv_naive <- cv.FuncompCGL(y = Data$data$y, X = Data$data$Comp,
Zc = Data$data$Zc, intercept = Data$data$intercept,
k = k_list, foldid = foldid, keep = TRUE,
mu_ratio = 0)
plot(cv_naive, k = k_list)
![plot of chunk cv_naive](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAyVBMVEUAAAAAADoAAGYAAP8AOjoAOmYAOpAAOv8AZpAAZrYAZv8yzTI6AAA6ADo6AGY6AP86OmY6OpA6Ov86ZrY6kNs6kP9mAABmADpmAGZmAP9mOgBmOjpmZmZmtv+QOgCQOjqQOmaQOv+QkGaQtv+Q29uQ2/+2ZgC2Zjq2Zv+2///T09PbkDrbkP/b25Db/7bb////AAD/ADr/AGb/OgD/Ojr/Omb/OpD/ZgD/Zrb/kDr/kNv/tmb/tpD/tv//25D/2////7b//9v////VZo7RAAAACXBIWXMAAAsSAAALEgHS3X78AAARZ0lEQVR4nO2dDX/bthGHlXrL5ilNls1zm662m63WYqpb1Gxr12lWJHz/DzUCfMHbgbwjQUon3fNzIhskCIp/3uHwRi6UwIrFsU9AoCGCMUMEY4YIxgwRjBkiGDNEMGaIYMwQwZghgjFDBGOGCMYMEYwZIhgzRDBmiGDMEMGYIYIxQwRjhgjGDBGMGSIYM0QwZohgzBDBmJFDsN0npX7ZaJ6BLWbTJyiP/bDJ8a5xElQWkJbaDThaWCZut22cEUiC0nbAtUKSQbBtcz7bn8AtPwOnVufZBl9l/+Nn9ctPfUlgWXBamAScC3R6qN3MffqfviQoTX+nXflvCOMF+/nH/31qTwPYcvhXfGZ1njZrw668vvveJKgsOC1MAs4FOj3kbir2EHBSmLbV+kE3CoJcLlHFt1a9ZQ/5yoRLRFsYUBaUFiYB5wKdHnI3XJmJNOiWw5BPMOgM9BZt+/H9majDyusSigMkDTUw6Fyg00PuVlZFUZlAEph2+Od/o2+AIp9gUK3S6hHZPyyYvizb3qTBNVjiXEbUY5i7BErbbwbqlVGwn4FToAqmr3Dw5YAkuCwgDdotr2ADD7YbGiNmFAyslfUWfcUP/wbzDLWwwYECcC7Q6eF2Q95dQNoIvfIJBlaiTTsM8OtwHbaN9wWShlZh4LlAp4fbbfDB4EYrDunpYIYIxgwRjBkiGDNyCPb5zZP64UpzB2ypNj1BW9oPJznaNU6CygLSkrsBRwvKxO32HGcEkqA0/Z3Ca4Ukg2DPzfk8/+YjtOUDcGp1nufwq7x7VD8ERwGSoLLgtDAJOBfo9FC7mfv0d31JYFr5nT6/eYzLRTBesA9Xf63N5N0juOV9fGb1Fpu15vMfyuv7ri8JKiuRFiYB5wIkYXeLPQScFKY9a/2gGwVBLpdYEt5GzZZ3gK9MuUSshQFlQWlREnAu0Okhd8OVmUiDbjkEGQUDzkBvMbYf3Z8JwcrrEokDJA01MOBcoNND7qarovh7RUlgmnr/p/gbYMgnGFSrtHpE9g8Lpi/L81Vf0uAaLHEuI+oxzF0Cpr0bqFdGwT4Ap0AVzFxh/8sBSXBZQBq0W2bBhh3s85uBMWJOwaBaWW8xV/w9LqzHWtjQQAE4F+j0cLsh7y4gbYReGQWDPEHbDgP8OliHPcf7AklDqzDwXKDTw+02+GBwoxWH9HQwQwRjhgjGDBGMGSIYM0QwZohgzBDBmCGCMUMEY4YIxowxgi2EKZhQsBF5hRQzCbYZcZwhXM1c3nyIYMwQl8iMsYLtv35S+9vF4mU8PUIEG01haZIyCKY1U7s/duYVl0inEmpT06RmEGz39mNlaR15RTAKrVmFaqkMgt2++P47bWFvI58oLnEIRaiW8m/38UHH4X5xrbZfxBOQRbAhhGrlFwyVV1wiChthFOaaiWCnT1H9m1wwJ+hA9noJMYA/PJKFCf3YBpdrXuISTxmtVGBe2QXbvTLurydKFMH6sRa2mdDCDvfVHPBt3DclLpGAGx5OKlgTbPT0dAj9wP7wSBYmLrGH1sJqgcYJBphPi+6qlzpsHEFzebxgh4chK6TFJRIonOZyBgtLWlF/XgGFrcDUtEEHMq+4xC7cCkwEO3GA/t4Mgh3uF+AsAExeoRfXH+YJOu5vyv/XRMVEMCxA73yOsL4ruE/nrRCXCOJNssko2HgLE8FSmK6oeDhF6rBTJTGccuwoUUjRJ9iy/Fk6+4/umurLWyEuMQauwFrBSpmWy+VmY/7ZXDN1TYlgIEXKwpaGzcZ23zdI19QxgUdUQrXogg1DBOsFGLBcLgO9tKZOFqnDjoO3ysERLJBrquGV3WtAURGsmyKysKVvXpMMr9SboI3iEruJKrBAronaYfvbUqpeCxNiwgFL2LwmaDjvb1/+XVwihaAB5vpDbVmeeQ0QTIccejVROvTYvYK8pQjWQTglwLEuNXbWVL9g6bxCgmBKgOsOR/climAT4I8wLx1/qOYRrGf1irjEBmhKgOcP1UlYmAjmEHb5LnMLlm5q9eYVAIIxS98fTj4ehlu9Ijj4HfSBeU0tmMytpwCs2VuqeQVDrl4RwRqCCizyhydiYUKDDg9tBRb7w8nrMNzqFaHB7aCvZgGo3ILliBLFJQJDYEvIH+axsJWZl3hNO0ERLMIbAgubXxkFGz/zVzC4Q2CNP4R758cJJnPrM+F00NfmZTQBBBrpEs3MX6JHFJfoEg6BOfUXKNB0USIu78ULpsIhsJ7hlGMLJvhDYKYCU/2CFcWmqP41R0G7xJf/+Io4/VcEq4nmZAfxPCiYkQmYSooNOnZvPwKdGYi8FRfvEgtvzLJrOKWyqI37RFIXbFhfCibjYSPw/KFSXv8G5AJhtRTJwiSsH0I8wmwD+nj8qxYspZaSBX1zkKzAgiCjMi3V7Y8krJ8eb4S5Y/zLizNSB5PFEFMSD1jW7jAWrDGvvis104K+y6UzoE+F8R3HG72gb/dqcbMCt4lgmri/FxiwxJqXyjPivPriqe/NEJfoEt3wsDUvE9MrWLBe69LkmNOxvZE5HTDegKUXHkYNZew1wgm2TrrEZk6HvHvFx9Zc1h+m+nurIRaMeSlsHfbV4/Y6MeK8v9WKraUOi/BXFPnhYdyzgT0qNqxv3hNG4NJdojec0tXfizcvhQ/ryx9wmSXy0BclWBH6w8T0+aaXF29eCluHlVptF4ubjh3l3Ss+4HRRsL+XYl5KBjCz43X2JpvLvkukIH2J+Wm6omx4qOLuw3bUi3jwHD0d8tx6H7g3CuzvpblDDcHC1lAdJnPrXZyGctg7n1rvQIUgGBjWy7tXAuLVKap+jp7tjSqG+kNFEgx6MamsD7N4fYeFFxzWBjewc8ODUofdQXvIu1ccCq933q+9Nq15UTs3PCSsz0PhorzO3mg4xdZhAxDBslGEg19w7zy9N8qD4hIXtBVil+YSwbnzUGcvtTfKAzm8cq1kfVgSoHMDrL0ClzgMWR+WB88fxtGGP89msD9UtPVhIyzsfCmCcMNvK4edvSPCwxpZHzYe+zb6pvWlQneYybyUdP7mwOvshYP5TOalJKwfR+wPXblSz0oZBUaw/e217pQnPvXhEgRT7ssdCrepnOg7HOcONRjBVjemy1DCeojaugBvCD8rZSyYB6t8XT02Wybh+ASdvaF5BX2HIzo3PJCC6bVhYybhnCn2QSlFUHl5fYdjOntDUC7xzjTEVtIOswTRRlR5he2uPOalsEHH4uVHHXlAlNHIi8eLnKptm16JgZTAJeYhx2IIbX+XLBhQeeXrigrJsRhCe8tLmiLguUNjXhvfHyp3VlRG69KgBQOn4LRTBNa/iiOSsxVMUwBNr3DeIXaNHomxgpUVnNkArIY4T5foj6UEetWxYRsV5ml6eYwWDHfo8xFM1WG8xpNLuW1l0ho9EtKXSMHO2fDNSxm9Kn9oG8qT3KeYhjNmgsDFLIYoWtkg87JTACYxL0WxsDU8zQ136PNxiVFXlFqm5h1OUj5WMPMivuGHPgfBvGjeta5GLyeSn8i8FFqwdXJx2EW9ygNqK5fRoXLiQ8TTh0aBnOaWNK9LWQyRDubreYduJD+ZeSmcYNuOtZcX9CqPAgrml6EznKzyqhkbJV7CYoiOtjLQ2TuleakM7bBzf5WHo5Y1r+AVvsVc5qWk4YwhmrOx8c3LmQIwrXVpKO2wy3v3itOzYSuvyLyoTx8aBVKw7cIMUw4+NFPBFNBQjs2rsbDpzUvhBDODyitiPwd/lwhVXmG04a5GmeemRE7C0RM7Mh/6tPECw6Ch3HZFOV2885iXQs/pWNyNE4yjSyyaptcyNK+m63COno0QbNCxuqA6DIzkXWdY+OY1ddPLAx8lHu4v4j3OfivZeRCAW3kpz7xmvRulHRbTtLusbW2Wrlwbrx9qTvNS+Knal/Fqe7DdZa3LCS2yrUahIoK1uL7QdYatebXDlEVgXrnPvouZBDt5vJrLCzT8hnL0+Ly5EcHCOAMI4x29sq9GoXLxLjGKCgFfuGyng+Zf3EAl16ypnkOfqGD+tDV/iZdjW8pxhtYlHocLDutTQYanllN3HS2S97hUwTr6CpXjFYtGrNNQS2VYvZL2lyfrEotQLS/IAGuu6ebKU0EIZp6qkqzDqqfk9B36dASLqq2EWq0rnHlEuQ+Mha07H1m/T70u+NRcYlEUHdUWbFubzbHD+BD0RNLu9wwMOPSsRF4wrLa8qLDw1Zp06jUVQh22YhnWF4FhFcoNKiq1NmHnrq/WUdtdIUjBVqlH/jb0rF45hmCQCwy9YLLamm0mLxVkHUYebu4/9EQUHgkXWHnBcL2rY1onqpbKECUOPnRWAJnarqbQBS7rFxs6wEHGCTlCy+iGM2L1SnXDEugpMp2t/blSsP9bAgK2XU9hkHF6aqk8z+nQ9M+t9x4DWaT/RmsaZ7aO7irWCfSCpx9khOR5Tkfv6pVlXY3UP97fCQH71HU0gOqoThcYeEEdxxdHmlRDZR4L67qYoRGEFmEFSXq59MGgLMaulG9akVrEyzgfs6xe8bu/p0arC9dhgQt0TMs1rBNWS83TW18UMwjmm1NQh5VKqYRUp15nhcw0vFJ01GGpv7O4yNioKrFUWGGxUEvNJZh7zUOjSPxNrbNipTahUTVSzb3gJCuzuMQZ67BqVu7mKvKAhTfEz82uLDNZWH7BNtVPdc3rFwhZFa5io2IvVcUcgrXXp779mx/3DzwKEqhjt1AmrkrVzGJhFD1G0yh1JhYVMksdNotAnkVd+bucEfPUYbEL7P0bubM7VdDh6tx0apmpHRbOSEf83b1zytzOnhmDjuyMODe+zGdhs3LUV9tPigjGjEudqs0WEYwZ4hKZMaVgDpvFvFzNXN6MTCiYZJo7kwjGLJMIxiyTCMYskwjGLJMIxiyTCMYsk/RWsEMEY4YIxgwRjBkiGDNEMGaIYMwQwZghgjFDBGPGMMH0k6mqNbT2N0Km/e0iXjGdollkTSjJZiKUtG0PTyjJZiKUpB+kMeA7VQwTzL6nhfDGlnZXfTXX19hszeulSe+GqTMRStq9fmp2xZdkMxFK0k/H2335SCupZZBgh4fH6DdCJvNQ6NfIG2v3+2/uiCXZTKSSVL0rqaQmE6GkrdbVSEUsyTBIMPNggTv/N0Im8wTv1EMYAw4P399TS7KZKCWVVCZCKanNRCyp2pVYkmGQYNqgq7vD/kbIpB/5gf1y65u6OiKUZDNRSiqrlur9TpSS2kykkpqnvJJKqiELtlosQk/f74n9TMi7UWcqd23iB3xJNhOhJP3p7Ir/TlUmUkn7W+dZosR6bHhYPzzwIPj7tZmp1349XEk2E7EOm+c7lVbpHnwOwbT9H7578n4jZNIOAR8lOt4NWZLNRCjJ+jTid6oyEUqyetG+U8Xgdljpuk1Qi387rc1EbocRS7KZCCXN9p0qB3BH/k4V0tPBDBGMGSIYM0QwZohgzBDBmCGCMUMEY4YIxgwRjBkiGDNEMGaIYMwQwZghgjFDBGOGCMYMEYwZIhgzRDBmsBcsMbdsfZOedGYnltk5j3peIQvOVLDy+tMEKxXOfWbTcDaC7W/Nwp3y49d/ftSaNBuqxT2Hh78sFjdbM7909a2Z0a53/eauWfwDvOPzJDkbwVY3Zl6m/njx6MzDNYt7Xj8d7q9Laa5N6uqlMT+96+Ku3j5o6c8xOBfBtETltdcfh4dHXSM5LrFM1WsO9D+tj5ZmdWd2vW9XJXHxiecimP6olSo/tK21gq0WpbN0BfvboxasNSuzXQSbiz4L29+a6iywsNK0agurt4tgc9FXh+mP3ZePnmB1VWbqsHq71GFzoYM8vejKRom/faiixHrDeqGDQU+wb80ShMO9iRKr7RIlHg9jM+R2MBOPeG6ClVZTLeChXn/p6RCmQQRjhgjGDBGMGSIYM0QwZohgzBDBmCGCMUMEY4YIxgwRjBkiGDNEMGaIYMwQwZjxf7WmoaIuJuV/AAAAAElFTkSuQmCC)
beta <- coef(cv_naive, trim = FALSE, s = "lam.min")
k_opt <- cv_naive$Ftrim$lam.min['df']
beta_C <- matrix(beta[1:(p*k_opt)], byrow = TRUE, nrow = p)
## does NOT satisfy zero-sum constraints
cat("colSums:", colSums(beta_C))
## colSums: -0.2672194 -0.1827702 0.1566777 -0.2521267
Nonzero <- (1:p)[apply(beta_C, 1, function(x) max(abs(x)) >0)]
beta_curve <- splines::bs(sseq, df = k_opt, intercept = TRUE) %*% t(beta_C)
Curve_L2 <- colSums(beta_curve^2) - colSums(beta_curve[c(1, nrow(beta_curve)), ]^2) / 2
Curve_L2 <- sqrt(Curve_L2 * (Data$basis.info[2,1] - Data$basis.info[1,1]))
cutoff <- sum(Curve_L2) / p
Nonzero_cut <- (1:p)[which(Curve_L2 >= cutoff)]
MSE <- sum((drop(Data$data$y) - predict(cv_naive, Data$data$Comp, Data$data$Zc, s = "lam.min"))^2) / n_train
PRE <- sum((drop(Test$data$y) - predict(cv_naive, Test$data$Comp, Test$data$Zc, s = "lam.min"))^2) / n_test
naive_result <- list(cv.result = cv_naive, beta = beta,
Nonzero = c("Original" = Nonzero, "Cut" = Nonzero_cut),
MSE = MSE, PRE = PRE)
random select a component of the compostion as reference level
## mu_ratio is set to 0 automatically once ref is set to a integer
ref = sample(1:p, 1)
cv_base <- cv.FuncompCGL(y = Data$data$y, X = Data$data$Comp,
Zc = Data$data$Zc, intercept = Data$data$intercept,
k = k_list, foldid = foldid, keep = TRUE,
ref = ref)
plot(cv_base, k = k_list)
![plot of chunk cv_base](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAA1VBMVEUAAAAAADoAAGYAAP8AOjoAOmYAOpAAOv8AZpAAZrYAZv8yzTI6AAA6ADo6AGY6AP86OmY6OpA6Ov86ZrY6kNs6kP9mAABmADpmAGZmAP9mOgBmOjpmOv9mZmZmtv+QOgCQOjqQOmaQOv+QZv+QkGaQtv+Q29uQ2/+2ZgC2Zjq2Zv+2///T09PbkDrbkP/b25Db/7bb////AAD/ADr/AGb/OgD/Ojr/Omb/OpD/ZgD/Zrb/kDr/kLb/kNv/tmb/tpD/tv//25D/29v/2////7b//9v///9ja24oAAAACXBIWXMAAAsSAAALEgHS3X78AAATRUlEQVR4nO2dDX/bthGH1XrL5rJNl85xm66V02yrtYbuFvY1y6pZi4Tv/5FGgAQBHA4kwDfxqPv/3DoGCZLgozscjoC4ESxS2pz7AlhpYmDExMCIiYEREwMjJgZGTAyMmBgYMTEwYmJgxMTAiImBERMDIyYGRkwMjJgYGDExMGJiYMTEwIiJgRETAyMmBkZMDIyYGBgxMTBiGgPY4cfyf6efih/eo1veFUXxI7pF1nqM3aC3HMCxDv7RkSJ5EUXhHhO55LgirAxtJdIKv+hYFL/CvcIaAdheXejbR7GH5222BOrIZj5GbtBb9uC2HMsd37knPpY38x1yD+D1IZccVxTYzZffCr9IEnz3b6QyruHA3v7w3/IWHr1Pl95y+tkzvHpLaQq/gM9bcIPeondodPi1IgSL/OsBe6G7RBVhZVgrkVYgRfJg0Gu0aCyXePj1bcAlHn1fJLQb/fl/uEvENgRcIgIsYGHwY4xcclwRVoa1EmkFUoQ0oE2jASuv9uDdJLWlvBj/81f1br8F+jBsQwCYcon+3fN5eTcFueS4InQ3pJVIK5Ci/ZmAoadtbq7n4es6ODB0Q0vQ8R/3Tsmbt/ecjNcRhZxpZ1HQKNxWIq3AGnYmCzv2AFZFbr9FbggBE16fgn9m38LjIZccVxRoLGwl0gqsYWfqw+TF4i5R3r7TL2idMcJ61WG5vRNqYUhUgFxyXBFShrcyJqw//fTbvFFifQuPyDCkGYcFRmijjMP2foe1R86IWARyyXFFWBnaymWOw1hzioEREwMjJgZGTGMAe//p6/L/r66u7tEt31+Veo3WKWttYzfoLXoHU+wdHSlSF3EFjolcclwRUoa2EmkFUvTi6g9vvN1CGgHYo7rQ77biEZ632RKoI5u5jdygtzzC2/JiK74HJ35x7xWpI4Ay5JLjigK7+fJbgRSVBL//CKmMaziw767+Lj/zL7xPl97yyvt86i3i/WdfbeM26C16h0bvP3ujCHlF/vW8ABeC7RJVhJUhrURagRWVB4Neo0VjucT3n/015BJf+L6ocaP3AZeIbAi4RARYwMLgxxi55LgitAxrJdIKvwhrQItGA/bptjo1suUe+fxVvdvzQB+GbQj1YdIl+ncP8YjwpiCXHFeE7+a3EmkFUqR86xmAoZ+T5uZ6Hl7XQYGhG1qCjufunZI379ELALyOKOhMu4qCRuG2EmkF1rAzWZj4Kh1YFbk9j9wQAibl9in4Z/Y7eDzskuOK8MbCViKtQBt2nj5MXizuEtXtezVdWC87rI/AjpiF+VEBcslxRUgZ3sqosP7V83mjRH0LXyDDkGYchocjYgxgZZyPDCf8MyIWgVxyXBFShrZymeMw1pxiYMTEwIiJgRETAyMmBkZMDIyYGBgxMTBiYmDExMCIaQiwDWsKTQhsQF1WSAsFVvSrdjXglETEwIhpocBYIQ0FdvzitTjebjZP/GdsDGywciNdNAIwyUwc/pxet0UX7xIbUEUtvWEEYIdnbypLS63bossGFqQlRgB2+8G330gLe+b5RHaJfWRcoIE1JjAhTneba7H/0J+IxcCSlTtSpMYHNkndy3WJrjNkYAuWCQdtWlMBs4KOyKwXyxd0hnQs7OIEQw1qwC7cJU5pYYePlfsbOUq8LGBt5jU2sNNdNRd87+em2CUmKNc/NaDJgOlgY+RMx8XJjjYmBTaRhV2WSxRjWxhiPo1kqp77sP4KpKKGATu9jF6yGX9oViMkFTXQwoJW1F2X1S3bH07dh01U93JcohQDoyKkAxsH2Olug84CiKnLahWSPhwh6Li7Kf//kEiMgcUoF64/HDGsbwvuw3X76UJcopgI2PwWdgnATLJ3dGDch00k5PEy0SjxQpQ3LhECKkSWyR+RyQjSqjI4NdVVt59W7xIbf4gDy7SKInPnuS00NbV6YMIaLkNgFi1/YiKnps6mOtQAHVgmIC1pZaYW92Fnkz8lIPNpFcrirFrch51BdgdmAct8XzjV45XDU4QoA2tRbvdhhRNnZNBLjvp4pd6EbWSX2CIPGOQ12TjseFuiGt3CVi+wOsUL46ccOB9vn/yTXWKsctCBKSCOdYmh09xkyCFXE4VDj8PHmLdkYAGZ9GENzHWHg1NT3cDCdVmIIDALlxhjEg4DG03O8i+k/6o6r1mAjbx6Zb0uMbctzA8PZwMWqNtXlwMs0y7RRBuDgYWHWp11WUYwQZ/Xw+XKyqxgftrnYdOsXlmpcvjMsnaGYj5gPLc+ReCJikn2Cnv0NSmwiVavXASwmlTdgcEnzrQsbHUC3+cggWnjcv3h5H3YNKtX1iiT77W7LxV0jAvsHFHiGl2im6A3wy+3AxvHwnZqXuJ12gUysFo58IcSWD32UmE9luwdBmz+mb+rk8kdAn9YCPRxyjBgPLd+sGx/aD9dziYBVs38TfSI7BKVEH9Y1LZVW9gELrGfGFgt53FKbmayibr/wgBRA7YSeesd6v5LE0O+RUUDy/P6B0wljXaJT/71eeL0XwYmlQMLK7RxNfNtALAconIm/kYHHYdnb5BkRkTdflqPSzT5Xssf1llE/wmzImXR8o8XG9aXwPh5WJq89Q6NPxSWP7SAadNSv0JHTbAwDuuTlQOXWPnDIsOBGQsLmZfgBX3TCvjDrPaHwk1ItQYZUAuNEsm7RG/6fJ3eMH2Y/fxLe8Giu+W8GGIy5W5YrzMbwB+6FtZhXmKxC/pIC07uVcCayQDCe6LiDsI6NHhB3+Hjzc2On4dBObnDJttb1H2YsCZkN1Fht3VJjfHEeffh67HfDEHXJVrmZY++rPShl+81fViExpjTsb/hOR22YO4wt3KH4IlKHcbH+kMRC+wh6BL1nA5+94qRAqbNq3GHQkeJhZug78hsQMX1YZ/f768DT5yPt5LYA/dhUjDZq80Lrl62A0ORYF4iPqzX7wlL0GW6RJjsze3eC0/Qjw+sDOvLH3SZZf9Dt4oiMBhtKH+oJ0aZx8tOvrc7FQUV14eVrPabzU3LjvzuFSmYO8xdXnb60EpuJH08F5qaIircvAqHl7+iaA3AyLlEmDt0Z/ZauNxnlFHZKEdjZDr4e+uV3FSUsM1LDMpGOUqwsAesD+O59Y2cZG9jXjrZ22SjzDPKpg9LUQIwNKznd6+gyV7Ye4Wmz6e5Q6kEYNiLSXl9WK3cG3oZ89K80PUOqUrpw7bYHpf+7hU/2dukNIydWeGGM21jKmD9dDkuUaDmZewsMBs73R8KBjaGbGAGkmVeQjjmZRL0PZTiEhNXiF2AS7SSvZVMgteNNvzJvb3MS0Q/XrkWvD4Ml52Kaoyq6cdg79UrueGI14cNlJfZKMy3EfnvcLOfqPRTyvqwGS2MgkB4CM3LBIdOVGiijQn7MF4fFpA3K0rAVFRhBRkDgnmjhUaJiweGmldmzQcI5g6HmJdYLDAK8ibZOKmoKndo+qxRzEvEATveXsukfOK3PqwXGMwd1pDMYgcRzB3OA2x3o1KGHNZbyk12PugOC21UiXNFWxXzxSpfVF+bzZNwLHWnogo47hrBvEQ0MLk2bM5JOMsV8IcaUu6komx/aKLE4eYlIl3iVg3EdjwOq6QzG7Y3dFNReVG4kfxI5iVig47Nkzcy8hj10K1aqEsEE0VNtNGkogR0hkNzh1CDw/oyfPzg/oLm1uuBMpLptSYCOHmoQakoqDFWr0iHeQFTBOzOy4y9zDc5NOYFw/ihY2VH0cDQKTjasnbXFwBMYMG8M2o2EwGSVsEmaSiwek7Hw+/8EHJVLhEfK9vu0J93OL55ieHAyohEbUCWr6wKmMDWe9lPmHNnns14mQ0oziV2yx16Wc7QjJUd87LW6I1rXVIxA2d+0QA+h61JReWueQ2aZNOleAt7wKe51Vrtu1dQ83KXpiATAaZwhpVigakX8Y166FYtB5jw5mwoF2imH0pGrnmNOlCGigT20Lo4rNehKcjNbFidl3DHygLOs5nIvET0NLewea393StwoJw5YX3de004UIaKAbZvMa8Vz613MxuWeWnr0s4Qm8Y2mYZGiet+9woYKNtDL+s55aQDZagxcolSa1sfZpuX4wwR85p0oAw1eOC8wnevoM8oMzh33p/GNmFsaJQyDrukd68AWk6skeWueY06BaBLkcD2G/XUa8xDt+p8wFzzsmlpeoX11RrOQHl68xJxwNQzyl1bnqPPoZco9ImyO1DOrIHxHANlqMhJOHJix8iHXqi8cVfmhPX68b89S3TSgTJU5JyOzXZeYOdwiei4CwSKElA9K2q+SN5RbNCxW3kfhj6gBJ1X9abKBs9skbyj+CjxdLfyxyserSyHgWIzRbQ4j3kJfoApZQeG5osNnbSGa155np/HvET8VO3VvtrepZU5tEyo4Tz2ssZd81qX1GUDQx+fVCNjENY7mQ1xLvMSiwU2h/B5hlXPZf9dTwzVvdfEUwC6dKnAAlFh5uagmoGyCS3mHihDLRTYpC4RT+4iYXz1lQCNdRXu8qHzaKGzpqYDFsg+IT2XoZW75jV3JO/ossL6YLfl06oeUJpx14yPvFp1IcByo6qPsofHmRfWF9q8QNd1bvMSiwU2qksEfZbIfEfo0ioKp+s667gLKgKY+laV4JyOcAe3AGC2XQkvk2G6LefV5U4Yb79E7/y0RJyFPbR9ZX31tUath3b80XzyIncXlUXLGjIDxweeKC9A0RNJg9COofc724fOHWxTwoPnwLyg020ZWhpRYV6hfP5xF1RCH7YbENY3gZm9Jq4FYOqdqQ9x5UaBvgtsaAkY0ze0Co/WueMMR5HAdqGv/I08dAaBdQAsOgwRA56rPizkAps+S+2YiVZaYqG0RHQf1gmrffVK3T80P6L6aQEogobo+ztgUdD/tfZZqG2JBUWFUEOjxKhDI/fPAMwMwCY11GGO7t9ZUxk/TZ5ZgDOXVuP7gCNcwhAZ1zzjsACxsArLHG0gyN+2rnxYVuTuouoIMhZpXmKc7+mImPnbYQQIsF7SwEqvl1sdVvkH8IKZhwoJMpZHS8w2tx56tbCR9FVzMKuPqywM7lm4aSc8yEi6iXNqKLCY1SvBgCEEEKMZ/Lvq+cDBs9zjVMFSyj3TIkJLzGNh7bEeGvjJe9YadKDKrjCL0naVgVjQNi0itMRsq1dco5F9Sp7hN71S0bYRYqqPlyFBR+MCa8Nq84IEaInzRolIWJ/XgR1Cs9nYAMrbOsDCCgQ9F0jLCzqaBVjwrraof0QiLQqNAt1YUNhFdDSnhSXc/fSwXlnUFW5R2qoKgi4QaqYHmNZTeWG7wJARdQMrXK9XcbryDcp1gYRJ1ZoDGH4XoaQjaxjYf5i/fQ6Iy3Mtag1G5WgWCwvc0xz8hP8orPfndkmgFrUCVJXOCWw06S+jKa7Wislonj7MuqfOTy808GAWI9OHDbjuZWumoAM+tNC3WeAMELr4QHfNtoRrOUFHoqldqi5hXuKqxMCIaaHAWCExMGJaKDB2iSFNCWyAin7Vroack4imA9avMlcaUomBEavEwIhVYmDEKjEwYpUYGLFKDIxYJc5WUBMDIyYGRkwMjJgYGDExMGJiYMTEwIiJgRETAyOmnsDkV1NVazLNvxIqHW83/grckPSi3YQzmUoJZ9o3h084k6mUcCb5xQw92qTUE5h5UUvCK1uaXeXdfLiOrabfL530cpi6UsKZDk9f613jz2QqJZxJftva4ZP7tDNp9QN2ennv/SuhkvpW6KeRH6zDn77cJp7JVEo6k6h3TTqTrpRwpr3kqlAlnkmqHzC1UH3r/iuhkvoK79CX+gGdXn57l3omUynlTKUqE0k5U1Mp8UzVrolnkuoHTBp09ekw/0qoJL9CIrZxDzd1d5RwJlMp5Uxl11K94CnlTE2lpDPpbw1NOlOldGC7zQZ6+m5P7FaK/DTKSuWuOn6IP5OplHAm+dvaNb5NVaWkMx1vrS8TTevHBoT1/QOPBH//oKbqNc2LO5OplNiHzdOm0irtg88BTNr/6ZvXzr8SKkmHEB8lWt4t8kymUsKZjE9LbFNVKeFMhldam5T6j8NK162C2vjX05pKyeOwxDOZSglnmq1NlQPYJrdJiTMdxMTAiImBERMDIyYGRkwMjJgYGDExMGJiYMTEwIiJgRETAyMmBkZMDIyYGBgxMTBiYmDExMCIiYEREwMjJvLAAlPLHm7Cc87MvDIz5VFOKyShlQIr738asJLw2Fc2jVYD7Hir1u2Uv37/l3vJRG+o1vacXv5ts7nZq+mlu6/VhHa565dbvfYHeWXkIrUaYLsbNS1T/vrg3pqGq9b2PH19ursu0Vyr0t0TZX5y18223t5n5c9ZtBZgElF57+Wv08t72SNZLrEslUsO5H+Sj0Sz26pd75pFSVR84lqAyV81qfKXtLUG2E4ucrSB/eNeAmvMalctgmRg86jLwo63qjsDFlaaVm1h9XYGNpe6+jD56/DJvQOs7spUH1Zv5z5sLskgT665MlHiH19WUWK94WEjg0EH2NdqBcLpTkWJ1XaOEs8nZTPJ42AiHnFtwEqrqdbvpN5/znSwphEDIyYGRkwMjJgYGDExMGJiYMTEwIiJgRETAyMmBkZMDIyYGBgxMTBiYmDExMCI6f9CG07gaAzk/AAAAABJRU5ErkJggg==)
beta <- coef(cv_base, trim = FALSE, s = "lam.min")
k_opt <- cv_base$Ftrim$lam.min['df']
beta_C <- matrix(beta[1:(p*k_opt)], byrow = TRUE, nrow = p)
## satisfies zero-sum constraints
cat("colSums:", colSums(beta_C))
## colSums: -3.469447e-18 2.168404e-18 2.602085e-18 6.938894e-18
Nonzero <- (1:p)[apply(beta_C, 1, function(x) max(abs(x)) >0)]
beta_curve <- splines::bs(sseq, df = k_opt, intercept = TRUE) %*% t(beta_C)
Curve_L2 <- colSums(beta_curve^2) - colSums(beta_curve[c(1, nrow(beta_curve)), ]^2) / 2
Curve_L2 <- sqrt(Curve_L2 * (Data$basis.info[2,1] - Data$basis.info[1,1]))
cutoff <- sum(Curve_L2) / p
Nonzero_cut <- (1:p)[which(Curve_L2 >= cutoff)]
MSE <- sum((drop(Data$data$y) - predict(cv_base, Data$data$Comp, Data$data$Zc, s = "lam.min"))^2) / n_train
PRE <- sum((drop(Test$data$y) - predict(cv_base, Test$data$Comp, Test$data$Zc, s = "lam.min"))^2) / n_test
base_result <- list(cv.result = cv_base, beta = beta,
Nonzero = c("Original" = Nonzero, "Cut" = Nonzero_cut),
MSE = MSE, PRE = PRE)
GIC tuned model
linearly contrained group lasso creterion: CGL
GIC_cgl <- GIC.FuncompCGL(y = Data$data$y, X = Data$data$Comp,
Zc = Data$data$Zc, intercept = Data$data$intercept,
k = k_list)
beta <- coef(GIC_cgl)
plot(GIC_cgl)
![plot of chunk GIC_cgl](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAulBMVEUAAAAAADoAAGYAOjoAOmYAOpAAZpAAZrY6AAA6ADo6AGY6OmY6OpA6ZrY6kNtmAABmADpmAGZmOgBmOjpmOpBmZmZmZrZmkJBmtv+QOgCQOjqQOmaQkGaQ29uQ2/+2ZgC2Zjq2Zma2//++vr7bkDrbtmbb25Db/7bb////AAD/ADr/AGb/OgD/Ojr/Omb/OpD/ZgD/Zrb/kDr/kLb/kNv/tmb/tpD/tv//25D/29v/2////7b//9v///9MZenrAAAACXBIWXMAAAsSAAALEgHS3X78AAAL9UlEQVR4nO2dDZujthVG2TbT2U23djZtEs+0tdMJbZdMmky2S9ct5v//rSKB+RRwL8iY1/ueJxt7MJKAg66ELHCQEiiCa28A0UFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWB4UPY8dn8L4qenZ98iAyfetK0PnDkIlvkKMVV8Omn6Mf/pWOLbNLnKYtMdu1djbtrORZJ8SAsNoUn2U5/+MX1iX3j/MSkOTYOVpId31YujoydZblK6S56+SRZZJa1kS0yFltL7Un769giMfOFvfz4n6z8Y7bPybPrkzQ/wo5PYrPNjb22uTRXdmTsKstVSneRK5XjVD/9q5OTbFG2bT93alhaxJPRRSJ8hUT3WV9sleNsOm9v85A6hMlrmKOU9qLjLy/t+OdYlFX0TiyVLcok/tclTLBpUry1Ydnmd49hrsVx6p+FnX76WF9qQ2K7VnQzdpYlqGBZ4/fJnhTDi2ycbtUf2aL0w8duG5YV0dk0xyIhvoSZrY87tTzX4mhbziqjj63FUfTv7mFpZewuS9CC9YZc19Gb0o5lmTmEyU4mIb6Exc79zrW8fHSnsWd3h1aj4sjYXZajlM6ixBFyfQrLO6aO3ZV2WcZZoIa5G+dnly/bOrV6VMIaJusVvHTjn2OROSNOP09YlDq69eIzToa3Niy7tOhuQt66uTat5wot7rZOjoxdZcmiTtK9AHIsMlvWTitb5LoOk6aUwZEOMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCYIywgl+CCwmakJX1QGABR7f1iwqIo6lmP/GAQvBqWEZbbojEPLCIsMsoozAsUhkIRE5cSxjZsgB/GV1lWGG0Nsz5hZBAkYax7KZAw9vZzRMIKrimMnccmVbAJd7XF4aa+0tWERbZ6UVhF7VDUhcXBBYQdv3xSpa1iIduwkoaw8D40Q/Ob9Pj277kwT21Ysi1G/X/TUdaXtqpatJVTuKiFxO/+ULxL3j3FXoVlxjJVjRo2+M3N2RZjYY1OpyMM7tO8hhX1LPXZS0y2d/+QhcSGLdaukq6w3ff78g/PNcxwfN0NiI60tNWDQ9jx7fvzHxcQJkxbNV0zcr1F1nodxl6hBxa9DqOt+Vx7LJFIWcNYIsmBGfzNYaDEEsYBewpDQzTHjcJwuMpE0n7YhmlYgTCigcIAWFlIJGNQGDAUBgaFAcCQCAaFAUNhYKxDGAc7BlldSORw4jAUBgyFgTFb2PF1sDm4ZmqzDfOGz5B4etilBzNZu5r4KE5LhOiFxXau/M6xQvLuyU5NNa/nNLKnIpFJSITFwb15OT04lJkaZmANWwiBsORPpYu/dVuqZGuMhfPaMDLI6rr1ZBi9sPDufehuw8bTEq+IhCVf7bP/nPfFTs+aTEIm7N1TVsco7FpMCInBq3182ZDIa+d+1tjp4OiUEEm3/tvyT0e3fnLWTShMyNwL5+lZN6GwAXwOTY2nlW4U27Be1tiGESEUBoZMWNZ83f3zq/3wusqsiRh1SDw9bI5v38d3nQF5QVoyH7Ww5N1TJqz2lZcICrsEihoWsoatAHEbFgRKXxTmjZV363lF1mbNwvjc5hFknY6tNhwKsnbAxwCPI6xhx9dBPqDoL+sufFZpDxNDYux6iqUwrQg+T7GHCcKWqmG0Ncb62jAyiHSko3eNrO692jtX4XWYL/RjiY+9475m5u/pYUNhl0Q/lrjt+S2Bc8063HNu/UL4uHslI/xtdw4chV2C2SMdydY+Vt0xuZ7CfOEzJI6mJfOZeOHc/B2rcWYKYx/fiULYgl9gRhwD7mNNQ1MlUe0fmd6GLXa7EYU1WfP3YZZiUJG+HKxSGDsc/Uinau/S8BXnJV4JfRtm55Dyhr5rMWEiqelvcCLpGpCFRNtNVPbqKewirLPTQRqsvlvf4TPvNYIJ4zBVnbUL41TFFqsWxp9Vz1GHxMP9+SFgGmYL4w91F2iF5fcZFV8ty/Ek7HO31ULynI7iXlnHIxHnZD0OJ5Y6kAgrvrhc/g5M2spRhsTzQ0c5NHUttG3Y8Y0d+33N5yWuAPk3ztqhRAq7CKu+DiM5YENTRCmsmIKzgomk7DXOr2EOmxe6GYJjwJbZIdHcazQ1rYbPeQxY3YaF92nYGxGTvoeGUZgvpowlhhujTYV3YWzDUvlYYibsyrOmaMuy6rFEksOxxEWZHxq0bVj+EwMxZ/5OwPvFiHAiKccSp+G/b8uhqQviawYRxxKXoLI11IaJGjh8YZGdnLPqbr7clm43AIVVu7leY6JgOCleoglrHIl1CosUwVAmDDgkNieWSsL/4oGzrP4yWzfehtV2U7yrQmN+1AritbAG9oAnTLGb8kbdz/Vt3URPXmVBU3v7YMKEtsoaKDkskaI1Gc5orN40G+CRPal9DhwSJTRPYtFh8SFsTHvjzg7xZuV/1T65NWH1wyZwIDnjJZV6rFlq2ZLvQ4ebEta+PUl0Ho+d8aONm6ATIbcVtfehze0Ia5/EGlsja41d/o5W50nx+ZbbsEgVcsQ9a0FvpNaNH8hJV907K9+YMFUno7P+8GrSXp+iyzeY03hm0MKiKBo4L/uSyNYfbklqq8y/4laFB2Rh5W6KTvJUcx6rqtfEzS8zGrXlNSSaH/lwfx+9nDBVz1vScqmC4bSNr+U1mpNPYQMzdBYSpupkKHqFQ9dl/myNFdZh9tz67hy45X5o4NyGja6Vyuw2G7i+lWqZTdnovg0UglzDZMg7JZIGTtaN94zXNqz/mfYrEzZ+gD124/1yY9dhI8hrg6ynv2TVcnD7whQHeO6wxxJ8BsJ84LljqC299p7CJFw5GFKYllUEwxwKk7BwN34IChNxXVsMiWBQGDAUBgaFAcCQCAaFAUNhYFAYAAyJIBTfqDeu2ilsvTjnK1DYeqEwMEphEYVhwDbsBqAwMCgMgKWuw4gnovofFxTGREsnojCwRBQGlojCwBJRGFgiCgNLRGFgiXjxCweFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYExTVgYnB+NU71TJEq2gfzHAM/PTlKUVCVSlBSX2StKqhIpSjLPx5uwTznThB123XfyROZoyn/RNix2TlFSmUhRkvnF1mJVeUlVIkVJ5jdij2/2upJKJgk7Pe477xSJzIPhxL9oe/z91ztlSVUiVUlpsaqqpHMiRUmx8WpVKUuyTBJmnxe2a75TJDq+fd//a94tTo9/fdCWVCXSlJSmRRXRlFQmUpaUr6osyTJJmKnQ+dlRvVMkMk/yk+5cuCmaI0VJVSJNSVnTkv9qq6akMpGqpCyAbtQlFaiFHYKgHenHI3EzkfBsNImyVc/9B3lJVSJFSea1tqp8n/JEqpKS7aa2SFfFpnfrp3c8FPE+tDP1yt2TlVQlUrZhy+xTVivrmS8hzNT/05+fGu8UiUxAkPcSa9FNWFKVSFFSFdOU+5QnUpRU+dLtU87k67AsdNtOrX2nTKS+DlOWVCVSlLTYPuUBYKfepxyOdIBBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBrywnrll4aZ/0lk1saya82jmFUJwo8Ky468Tlhn2vWWX4WaEJVt740728sUf98bJ+YP85p7T41+CYBPb+aWHb+yMdrPq17vzzT9mJigCNyPssLHzMs3Lq31tHq69uefLp9PDfabm3i493NnqZ1YNdsXnk279uQa3Iswoyo69eTk97k2LVAuJ2VJzz4H5Z/wYNYedXfWhvCsJJSbeijDzUpjKXkxdK4UdzF2OdWHf742wslod8rsgKWwZxmpYsrXNWauGZVWrqGHF5xS2FGNtmHk5vtk3hBVNmW3Dis/Zhi2F6eSZm66qXuLvHvNeYvFBGJjOYEPYN/YWhNOD7SXmn7OXeD1snVFfB4NExFsTltWa/AYe7fHnSAe5DBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFg/B9sjojWQHOdNQAAAABJRU5ErkJggg==)
y_hat <- predict(GIC_cgl, Znew = Test$data$Comp, Zcnew = Test$data$Zc)
plot(Test$data$y, y_hat, xlab = "Observed response", ylab = "Predicted response")
![plot of chunk GIC_cgl](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAe1BMVEUAAAAAADoAAGYAOjoAOmYAOpAAZrY6AAA6ADo6AGY6OmY6OpA6kNtmAABmADpmAGZmOgBmZjpmZmZmtv+QOgCQOjqQOmaQZgCQkDqQtpCQ2/+2ZgC2Zjq2/7a2///bkDrb25Db/7bb/9vb////tmb/25D//7b//9v///+5QTkNAAAACXBIWXMAAAsSAAALEgHS3X78AAAIUklEQVR4nO2dDVviRhRG47ZKt4JuF2270G7aCjH//xc2MxkISCLzkYx523OedTWQO2NynA/CnVDUIEXx0b8AhIEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMVKEFTAFEwpLiIUhECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBgIEwNhYowjbP/5e3QsBJEqrFq5t2k+XShD2BQkt7Bq1aiihWVjhC6xWt3+hbBcjDKG7ReXHSLCpoFZohgIE2MsYdXjsVP0zO+BKGhhYiBMDD9hr0/F7d8Pm1GLhii8hL0+Lff3L7vblzGLhii8hDUzikbYybyig0tTmQloYWVvC2ueiysaovAew4pioEeshsY2hE0Bs0QxECaGn7Bm+CqLYj1q0RCF3yzxYdP8633TK75oiMJ3Wt+0MYTNAc8usbjZ7OgS5wCTDjEQJoafsN3A1aeUoiEKv0nHKnD48ikaovCdJY5fNETh1yVuhy7wJhQNUXh2iYxhc4FZohgIEyP9/bDYoiEK33eca3vJfsyiIYqAaX3o5B5hU5DcwvaL4mbTaxNhU5A6hr0+ra1PhGUieQWmFbW9Q1gmUoWZFtZQ/nD57ibCpiCgS7zr3aNa2ctWJYmkeQiZdPQbiy0aomBaL4ZnTsddfaWFsaAvEyFX6wOv1yNsCrj4KwbCxCBVW4zUVG0W9GUmOVWbBX15SU/VZkFfVph0iIEwMbhPhxjcp0OM1Pt0xBcNUSTfpyO6aIiCvEQxmCWKwXIjMfzGsOfAGb1P0RAFy43EYAwTI4+wNsHDK83D7HTcscsNGYw93bfdOFZ2HjOwdVHwRTrKeaGDJb53RGd7HjYuwwcKPI++UpnPb3Q11h6r+98j6LijPVHtt6HY033bL1dZa+/t7zAYdvJMURf1mzNU9B2B1/G0v0tRd+e8ePPVcywDh3jcere2eBB2KBBhl/sjrHs6LsntrGjGsIuT0Rc+3hi2JVV7LngJI1V7PngJu7ICc6C7RNgU+HWJw8uN3PqwuufdTYRNQeos8dBPsgIzE6nCaGGZ8e4Sh5Jwhi8MI2wKfCcdJOHMBN9pPUk4MyGghb2bhMMKzEwETOtJwpkDWS7+wnhwaUqMkKv1zBJnQEAL64UVmJlJHsNYgZkXP2Fm9WV505+cyArMrPh1idYJH0c1BzxfOJsrvFyamgN+XaKdWgSmdCBsEnjhLAbCxPB54fz4J4shZgMtTAyEiZEn8xdGg8xfMXh7RYzkzN/ooiGK1Mzf+KIhCmaJYiBMjOTM3+iiIQoyf8Ug81eMsTJ/w4uGKJIzf/eL4q5k9Uo2kteHPW/sNaumBQbHQgTJeYnNU+WSFZjZ8BvDhm+/RwvLjF8Le+cd52YMWzKG5YMrHWIgTAwfYebeKdc+O4wVmJnwEGbyfkPTtD2Khii80ty+R92mGWFT4CvsaahP5F5TeUkVxp1wMpOa5sa9pjKTfC2RFpaX5Ndh3GsqL7xwFgNhYiBMDISJgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBhTCoMpmE5YDNH1EZhYeCTzPAtCgQgTC0SYWCDCxAIRJhaIMLFAXvyKgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxcgmrVofPMeh+Cgwsi8BP6TzctDOwxi4wsMbuxq6BNXaB12vMJMzcYqz9iM3up8DAenvtLrdv2LlDD6yxCwyssXrY1PufNuE1doEeNWYSZm7e1/7ddj8FBobeM3V7862tJrDGLjCwxp1RZM94YI1doEeNmYSZjymo7OeTdT8FBtq7ygX9ybtzFlhjFxheY9wxdoEeNWYSZu622P5O3U+BgabTCPubd+c9sMYuMLzG9sNfI2p0gR416rQwS9CoktrCgmusVss6qsZDoEeNOmOYJUZY6BgWLWy/cPuG1ngM9Kgx2yxxeZwlLsNmiYfdTT/z+mvEeQ+ssT7rSwNq7E57YI1doEeNeV+HmTMR9TrMBDavUW6CBnITFFFjFxhWY2nzdtfhNZ4EXq+RKx1iIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDH0hJUuOdY3/+zwyeH/EeSElSar1uRdIkwCk6PZ5qFtv9qWtmvX51Qr823/8xeTX/b6vGm3zcM/fjHCmmc+fXcPtiHV4+82pcxF3v92WdwcURO2s/l+ptVsb18ab0agydnc2tTN/WJtUtP39y/ttnl4ZztQm6zZPuhCqtXty66xYnZpilrYb+fFzRE5Yce1YqZL3K5dCrv51pxsm4vZnOul2zbnv+0SbXqne9CFrNa2KT7aJ1zu6HlxH3WQ76EpzEj4Y2O97RemY7PrdG7a037/j+0RD9vtYOfycW1irQt5sAWYJxpvh2Tfs+I++Fh7URN2MoatD/MJ15W5x1+fv90ftt+0sK7NmJCHzWULOy9ujqgJq0szX7CzxDtzhs2YZr7cSGROeVksa7d9MobZtmYfdCHV6u5N5EVxH32sfcgJM/O49nXYV9tpbY/TOtPZGS12uXC73TSvwyzx83EnF1I9/HI6S3Qt7Ky4OaInbCwC17TOBYSJ8f8VJgrCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDH+BYtGjRAeIvR3AAAAAElFTkSuQmCC)
ignoring the fact of one-sum for compositional data: naive
GIC_naive <- GIC.FuncompCGL(y = Data$data$y, X = Data$data$Comp,
Zc = Data$data$Zc, intercept = Data$data$intercept,
k = k_list, mu_ratio = 0)
beta <- coef(GIC_naive)
plot(GIC_naive)
![plot of chunk GIC_naive](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAtFBMVEUAAAAAADoAAGYAOjoAOmYAOpAAZpAAZrY6AAA6ADo6AGY6OmY6OpA6ZrY6kNtmAABmADpmAGZmOgBmOjpmOpBmZmZmZrZmkJBmtv+QOgCQOjqQOmaQkGaQ29uQ2/+2ZgC2Zjq2Zma2//++vr7bkDrbtmbb25Db/7bb////AAD/ADr/AGb/OgD/Ojr/Omb/OpD/ZgD/Zrb/kDr/kNv/tmb/tpD/tv//25D/2////7b//9v///8gEBskAAAACXBIWXMAAAsSAAALEgHS3X78AAALtklEQVR4nO2dC5viSBWGGbXtmXWkd9ZdpVuFtTfqsLOXcUUR8v//l6kkkNsJOadSCfXB9z4zD3SlLiEvdUmlEhYpgWJx7R0gNigMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwQgjbf0rTX7aO/wpb8k2fpDTVSxXcjdoNksoSwvqiCbm1y9RF23UTCkFS2F44VkoCCNud9mf3s7jls7BrZZpd66Mcvv9f+svPQ0FiWXJYO0jYF2n3VNHy7+m/hoKkMPeZ9tl/H8YL+/z9fz6dd0PYcvyxu2dlmnPSE/vs+B4Gg6Sy5LB2kLAv0u4po6XdFkIOaoftnD/pi6IgVJOYdr9a5ZaD1Fb2NInqGiaUJYW1g4R9kXZPGU1XZk+Y9JXTEE6YtAdui6v73e9nTx+WHZe2HCHIt4JJ+yLtnjJa1hV1yhSCxLDjD//ufAIV4YRJvcrZR6f+y8LcYdkNBnn3YD37MqIf03xLpLDD1tNXQGGfhV2wCnNHuPXhhCC5LCFMihZWmGdme98xYkBhYq/strgjfvxJTONbw7wHCsK+SLuni6b8dglhI3yFEyZ2oqfzMKFdl/uwXTeuEOTbhYn7Iu2eLpp3ZvJJqw7OdIBBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGxhhhCzIFEwobkZb0QWEAbGvv5xa23W6l4PvmO4fi1TGzsG3a/L4QK7MKK6oXhY1hNmHb3Fb+f0Sed0zZJs4jrKpa7MMEvhuOMq+wWtWiLoFohbF6ycQojG3hBaITRluXUQkr4UxHPFTf6mRVC06W9UhXFcaKV6d2xlMXtlvEIWzLs+gWDWHJY+Km5pfp/v3fC2FB+7D9F6+mtLXzMl0Bt03potYk/un35bvDh9ddUGGHp/Iyza86yvrT8rysSWfQkSwe06KGlfUsDVfDDk+ZqkYNG7rUVs1RsQ8r6Apbfbs+/xG2hqVO2cM/DE1iZUuR930gCNu//3j6I7iwrAt7220QLwmjrSZxn4dd7rmocoArDOsvOeE4ZIjIZjoorJdZ5xLV3Kew+CZ/9dxlHwYm7C4dNUASxhnF1C1jSxVr3Yq4VxbGdTkarrmQtLsrFGbi+sLYh5mIow8jF4moSSQaKAwYCgODwgBgkwhGtMI4ZBwmGmFc9qYjDmFc9naR6JpELnu7TLTC2IcNE4swDjiURCGMti4TtEncv10sN9JKbZ6HBcMubJcvvV4JEY7Pq3TjFmtXK1WVWRMvNMJ2i0f3cnwWlB0+vOZrid3rKY3uMVbEC4Www9fnyvO3TsPnapiDNWxCgvZhhydnLGEfNiF2YcnDx0Tuw4bTkqCohB2+XGf/xNss/bMmXuiEfXjN6hiFXQuPJnHxZr1jk3gtoptLJHo0w/o/nv/sDuv9syZejD1x9s+aqAk5NTWcloyGfRgwFAaGTljWfT3888v15bjGrIkac5N4fF7u33/cPXTmdxVpyXjMwg4fXjNhtSsoKihsCgw1LGENiwB1H7ZYGH1RWDA4rAeDwoDRDTqerM2hImvihbKG7d8uignFcFkTNZ5N4k56KKIyLRmDhzDWsFhgHwaGdqYjfNZEjX0u8aV/3jdrLN+sRacUFgr7XOJTz6Ppi5W/x+clhc3F2BPnQtTmkWvrZ2KssHJtffLr7qJFCgtFyCYx25Y/B19YXE9hofA8cW7+jtUwI4XxrkwRg7AZL2DyoR29RDg1xYd2tPHtw+ZZl8iHdnSI+3oYH9pxkSiFccDRj3ap9ipN3sy0LpG22tj7sHwNKW/ouxYeC0ndeIMLSWNA1yTmw0TjqJ7CJiG+QQfpEPewXsKNQ+54LAImjDMfdWIXVv3oc0phjqiFnWzduzBzk7h5PD1TysJoYTVb7MNOaIQV9xmVVyr1BBJ2z6YENM/pKO+VFZ6wNybrYTirKKARVl64nP8OTNoqMDaJp2dYcmrqWlj7sP27fO73LR+sEgH6K87WqUQKm4Soz8NIAdjUFBrhh0pGYeUSHHOjeKfCJp6SGb22vmvzrtfWT76gcnST6G5d8U17azRmqkPmW3uvEpY8pklvi3joe2jY3QmbbCrNZy4xWTptJu5M2FTVq4V2LjETdpOrpgJVh8rW1LNpUc8lTk2wEcLEtm52LtF0wII1YdvJG0NrH1b8xMBurpW/nliqy7Z9LTtAqTNdWlAuJI15LnFboP6Sn4/x6Fas3nONtqXMAn1qqm5KI6y5qGfcYQ5USc+Z9WZzM3OJnUU6gwIay0TGHeahgaG9Q71xYTVbWgHhBt86W8pCTIMWTGHNQcNWVVtCnioN5WTqUG27hSbsbMf0Mc2HRZlZbxRDh6qIDNsk2vuscyrDl16X2cXCNTtnGP6gCjP3WWnYUaFhGK+1ZR/+QArTHnn/wzJU/qhsWrtlSw0mTDvCSDu2BiLq8wtwrm09a0dtEvWHyvIlVo2/A9ryaJ9hhSmwjyIHhQUaYpra5/7tNybMYxQ5JCzIENPUbV08675dYZb2c2jjqCGmaZAhnuTdcJMY7vQ4z009jL+chXa3+s4Yb1iYYRSpyGmsfeMgQxX51oQFoj2n4peJsX1WVWcKkxhxYtvJxXAykvZEDtokuh/5kK9HowoLNlNsHqr0RQ4p7MIKHUhhAUaF7czCMnptfXcNHO7a+gCjwslhDTsRZJwxDUH7sP5n2iMJGzWBPjm3fB7mR8y2WlCYI+wESUg6u0RhjkCjwqDUrjuwSewQm63mdQcKA6BvyEphkdLXrVJYrNSuO7BJBIPCgKEwMCgMADaJYFAYMBQGBoUBwCYRi8ZsB4VFT3PdNoVFD4WBsW1cYKGw+GEfhgyFgUFhAMx1HkYCsa3/MaEwJpo7EYWBJaIwsEQUBpaIwsASURhYIgoDS8STXzgoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMP2HJ4vRonOqdIdHhaaH/McDTs5MMJVWJDCXtztkbSqoSGUpyz8fz+EwFfsI2q+47fSJ3NPW/aJuUH85Q0jmRoST3i61lVH1JVSJDSe43Yvfv1raSzngJO76sO+8MidyD4dS/aLv/3VcrY0lVIlNJaRnVVNIpkaGknfOaqzKWlOMlLH9e2Kr5zpBo//5j/695tzi+/PXZWlKVyFJSmpZVxFLSOZGxpCKqsaQcL2GuQhffjuqdIZF7kp/2wyXLsjsylFQlspSUdS3Fr7ZaSjonMpWUNaBLc0klZmGbxaLd0g+3xM1Eym+jS5RFPY0f9CVViQwluddaVP1nKhKZSjo8LWtBtirmP6z3H3gY2vskX6l3/ni6kqpExj5sns+U1cp65nMIc/X/+OfXxjtDItcg6EeJtdZNWVKVyFBS1aYZP1ORyFBS5cv2mQq8z8Oypjsf1ObvjInM52HGkqpEhpJm+0xFA7Ayf6YCznSAQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBga8sJ61Zcmyf9FZtbCsWvPo1hVCcKPCsuNvE5YZDr1n03Azwg5P+Y072ctv/rB2Tk4bipt7ji9/WSyWu3x96eabfEW7i/rV6nTzj1sJisDNCNss83WZ7uXNurYON7+554vX4/NjpuYxD9085NXPRV2syu1et/5cg1sR5hRlx969HF/WrkeqNYlZqLvnwP13fpyazSqP+ny+KwmlTbwVYe6lNJW9uLp2FrZxdznWhX27dsLO1WpT3AVJYfMwVMMOT3l31qphWdUqa1i5ncLmYqgPcy/7d+uGsLIry/uwcjv7sLlwgzx301U1SvztSzFKLDckCzcYbAj7Jr8F4ficjxKL7RwlXo+8zpjPg0FaxFsTltWa4gYe6/HnTAeZBgoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAw/g+Jt9zYBG7dJwAAAABJRU5ErkJggg==)
y_hat <- predict(GIC_naive, Znew = Test$data$Comp, Zcnew = Test$data$Zc)
plot(Test$data$y, y_hat, xlab = "Observed response", ylab = "Predicted response")
![plot of chunk GIC_naive](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAe1BMVEUAAAAAADoAAGYAOjoAOmYAOpAAZrY6AAA6ADo6AGY6OmY6OpA6kNtmAABmADpmAGZmOgBmZjpmZmZmtv+QOgCQOjqQOmaQZgCQkDqQtpCQ2/+2ZgC2Zjq2/7a2///bkDrb25Db/7bb/9vb////tmb/25D//7b//9v///+5QTkNAAAACXBIWXMAAAsSAAALEgHS3X78AAAIUklEQVR4nO2dDVviRhRG47ZKt4JuF2270G7aCjH//xc2MxkISCLzkYx523OedTWQO2NynA/CnVDUIEXx0b8AhIEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMVKEFTAFEwpLiIUhECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBgIEwNhYowjbP/5e3QsBJEqrFq5t2k+XShD2BQkt7Bq1aiihWVjhC6xWt3+hbBcjDKG7ReXHSLCpoFZohgIE2MsYdXjsVP0zO+BKGhhYiBMDD9hr0/F7d8Pm1GLhii8hL0+Lff3L7vblzGLhii8hDUzikbYybyig0tTmQloYWVvC2ueiysaovAew4pioEeshsY2hE0Bs0QxECaGn7Bm+CqLYj1q0RCF3yzxYdP8633TK75oiMJ3Wt+0MYTNAc8usbjZ7OgS5wCTDjEQJoafsN3A1aeUoiEKv0nHKnD48ikaovCdJY5fNETh1yVuhy7wJhQNUXh2iYxhc4FZohgIEyP9/bDYoiEK33eca3vJfsyiIYqAaX3o5B5hU5DcwvaL4mbTaxNhU5A6hr0+ra1PhGUieQWmFbW9Q1gmUoWZFtZQ/nD57ibCpiCgS7zr3aNa2ctWJYmkeQiZdPQbiy0aomBaL4ZnTsddfaWFsaAvEyFX6wOv1yNsCrj4KwbCxCBVW4zUVG0W9GUmOVWbBX15SU/VZkFfVph0iIEwMbhPhxjcp0OM1Pt0xBcNUSTfpyO6aIiCvEQxmCWKwXIjMfzGsOfAGb1P0RAFy43EYAwTI4+wNsHDK83D7HTcscsNGYw93bfdOFZ2HjOwdVHwRTrKeaGDJb53RGd7HjYuwwcKPI++UpnPb3Q11h6r+98j6LijPVHtt6HY033bL1dZa+/t7zAYdvJMURf1mzNU9B2B1/G0v0tRd+e8ePPVcywDh3jcere2eBB2KBBhl/sjrHs6LsntrGjGsIuT0Rc+3hi2JVV7LngJI1V7PngJu7ICc6C7RNgU+HWJw8uN3PqwuufdTYRNQeos8dBPsgIzE6nCaGGZ8e4Sh5Jwhi8MI2wKfCcdJOHMBN9pPUk4MyGghb2bhMMKzEwETOtJwpkDWS7+wnhwaUqMkKv1zBJnQEAL64UVmJlJHsNYgZkXP2Fm9WV505+cyArMrPh1idYJH0c1BzxfOJsrvFyamgN+XaKdWgSmdCBsEnjhLAbCxPB54fz4J4shZgMtTAyEiZEn8xdGg8xfMXh7RYzkzN/ooiGK1Mzf+KIhCmaJYiBMjOTM3+iiIQoyf8Ug81eMsTJ/w4uGKJIzf/eL4q5k9Uo2kteHPW/sNaumBQbHQgTJeYnNU+WSFZjZ8BvDhm+/RwvLjF8Le+cd52YMWzKG5YMrHWIgTAwfYebeKdc+O4wVmJnwEGbyfkPTtD2Khii80ty+R92mGWFT4CvsaahP5F5TeUkVxp1wMpOa5sa9pjKTfC2RFpaX5Ndh3GsqL7xwFgNhYiBMDISJgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBhTCoMpmE5YDNH1EZhYeCTzPAtCgQgTC0SYWCDCxAIRJhaIMLFAXvyKgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxcgmrVofPMeh+Cgwsi8BP6TzctDOwxi4wsMbuxq6BNXaB12vMJMzcYqz9iM3up8DAenvtLrdv2LlDD6yxCwyssXrY1PufNuE1doEeNWYSZm7e1/7ddj8FBobeM3V7862tJrDGLjCwxp1RZM94YI1doEeNmYSZjymo7OeTdT8FBtq7ygX9ybtzFlhjFxheY9wxdoEeNWYSZu622P5O3U+BgabTCPubd+c9sMYuMLzG9sNfI2p0gR416rQwS9CoktrCgmusVss6qsZDoEeNOmOYJUZY6BgWLWy/cPuG1ngM9Kgx2yxxeZwlLsNmiYfdTT/z+mvEeQ+ssT7rSwNq7E57YI1doEeNeV+HmTMR9TrMBDavUW6CBnITFFFjFxhWY2nzdtfhNZ4EXq+RKx1iIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDH0hJUuOdY3/+zwyeH/EeSElSar1uRdIkwCk6PZ5qFtv9qWtmvX51Qr823/8xeTX/b6vGm3zcM/fjHCmmc+fXcPtiHV4+82pcxF3v92WdwcURO2s/l+ptVsb18ab0agydnc2tTN/WJtUtP39y/ttnl4ZztQm6zZPuhCqtXty66xYnZpilrYb+fFzRE5Yce1YqZL3K5dCrv51pxsm4vZnOul2zbnv+0SbXqne9CFrNa2KT7aJ1zu6HlxH3WQ76EpzEj4Y2O97RemY7PrdG7a037/j+0RD9vtYOfycW1irQt5sAWYJxpvh2Tfs+I++Fh7URN2MoatD/MJ15W5x1+fv90ftt+0sK7NmJCHzWULOy9ujqgJq0szX7CzxDtzhs2YZr7cSGROeVksa7d9MobZtmYfdCHV6u5N5EVxH32sfcgJM/O49nXYV9tpbY/TOtPZGS12uXC73TSvwyzx83EnF1I9/HI6S3Qt7Ky4OaInbCwC17TOBYSJ8f8VJgrCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDH+BYtGjRAeIvR3AAAAAElFTkSuQmCC)
random select a component of the compostion as reference level
GIC_base <- GIC.FuncompCGL(y = Data$data$y, X = Data$data$Comp,
Zc = Data$data$Zc, intercept = Data$data$intercept,
k = k_list, ref = ref)
beta <- coef(GIC_base)
plot(GIC_base)
![plot of chunk GIC_base](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAt1BMVEUAAAAAADoAAGYAOjoAOmYAOpAAZpAAZrY6AAA6ADo6AGY6OmY6OpA6ZrY6kNtmAABmADpmAGZmOgBmOjpmOpBmZmZmZrZmkJBmtv+QOgCQOjqQOmaQkGaQ29uQ2/+2ZgC2Zjq2Zma2//++vr7bkDrbtmbb/7bb////AAD/ADr/AGb/OgD/Ojr/Omb/OpD/ZgD/Zrb/kDr/kLb/kNv/tmb/tpD/tv//25D/29v/2////7b//9v////wNpY7AAAACXBIWXMAAAsSAAALEgHS3X78AAAMFUlEQVR4nO2dDbubthmGybrsJF3m03br5nO22W3GGvqVLAuLN5v//7uG+P4Q5n2xDDziua+m9gGEgNt6JWQJBwmBIlj6AIgOCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwXwk4/pf+7/Bz9+D/rmo9RFP1kXWNSfZauKNecOvs69fduWWQOIora+7QcsmyRbZn1LC1n0V90jqL33a2GcSAszg70w+ck7uZbrRlIY07zs3BFuSbuXJZzuuHHdsbn9GJ+tFyD7vFZDlm2aGCzPv2z6C8yBj/+y5LYzu3CPvz4n/QSnnufrnLN5ZdewSvWpEXh187nbXBFuabcoOL0PjfUXdQ/ns5W1k1Ei2zLbGdpOQvLIrOzbtS4gquQeHr/YSAknvuxKCnD6C//tYdE24qBkGgRNlDCuh9jyyHLFtmW2c7SchaWRZYTuIYzYenRnnoXKVuTHkz/85fXbp8G6jDbigFhWUjsX72+r95FsRyybJF1M8tZWs7CsiheSJg12+ri9iJ8kcYuzLriSqPj3+0rZS5e3AsyvYpoKJiOLhosFO2ztJyF7cQWKmHnCcLyltsn4YohYUmvTrF/Zj9092c5ZNmigZPtnqXlLGwntlAdZg7WHhLN5bv8ak3jolmfVVjt2slawiytAsshyxZZltnPUtKsv/z8ad5WYnEJz5bbkOo+bOAOzcl9WNyvsGJLjpYSYTlk2SLbMutZrvM+jMwJhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGLcIC8g9uKOwG9KSISgMgKjxHk1YFEXjG6HxD4Pg1bCoMP3Vj5L25217LClswtWnsMWERXnxojApRUxcRlgtSx0TvazDFBstICy75vk/H6/+BNYtLLeVCbth516xZmHN4nXDvv1itcJoy45IWMFswiLaGqO+LuG+sTjcNTeaS1jRImTVNUzj0jSFxcECwhrNeBavIVrCwofQdM3vktObf+bC5qzD2Iy/TuGiERL/8ofi3fnrt/EiwkRFa6tCe42OMHhI8hJWlLNkAWEjW0/rqvKDvrD9d4fqjwVK2GjRmd5V5QUWYac378o/lhA2QqP5KBPmWexc5X3YMNq+xQ3HzsSBsNOrYHcMgt+8nZA2mdL5sfG7uVuFXZ72yTG11Yi54rRJs8Ev2Liwu1Fhjuqw9D4hqxXNa5lGNiArQ3H1y0jo4e33nJ2/poQZppcw4dWv1XpmK5m5t/78aIyFN9Vhkq18joSr/XplGv53RIrGuIEI43iCBGggaWnL21g4gfUKa9nabunqsk5hEW01WXlIrBoZtFWwZmGt9jtt9VmZsFYgpC0LqxTGotVmxSGRgdDGioXR1hhrE0ZGoDAA1hwSiQUKA4bCwKAwABgSwdALi7NBNfvr2w6kJU6RCIuDB/NyedIpo7B7IBB2/lM1IOr7/lCbybsmYliHgaEXFr58F7IOWwUiYeevDul/py81AXE+YdvqL5YJ+/ptWsZWKmwDU1kmhMTgxSFeZUjcxNwjbxodfg/gtiNp1v+5+nNNzfqNPqgF9sZ5S2OCnXZNnV6l9Vtrfphw1zeypeLlsg4z88MuT7tFhG3DVgcXMzCT5PgwcQbmDWzRViIVllZfL3/46mDZoJiBGX7Rv0ljT4cr1CExjXmnN+/il71ZsYmZgZk998MyBXNhYR6VQLWwNN6lwizVlCDtInj8HCRFCQutJWzyru+G389BEtdhQaD0tZQw9XOQ1o83XVMWvJwT7a2wLfRWyRodj9pwKNi1a7ZgKxGXsNOrIO9QdLdrl7TnRM+X70xMDImx5Wk30rT3xfc50ROErbuE+f6knCZe1GGFLP9tJfKeDve7donnsvR9ic+2ft8x2PnrCn1f4mP2dclqGx1bwqsb5y3guzAvarfNhERfvmSZeOPc/h2rcRYX5uWgRYUwoC8wc7YubL1dUwN40/kxtQ5b49j6q3hhK/H4+7AtQGFgSIdq75PwhbJ/isJcoa/DsjGkK53QtwEmDCQ17Q37QNKxtMQtspCYNROVrXoKu0srlY2O++GsX4zN+jlw2JFJYTNwr34xCrsLRfFiHYaA+1+nU4fE40P5S3w9im7GZhtyphmYYmbuUZSOuVMcllZYPs+omLnXxUxwvpJ2eeb9DlM812nqYUme01HMlbX8LmnSWK3d9VzMKkw6dWZ6E1IirPjiEu4LzJwZhcmaGvo6ThkSy1/+Re2amqkOk2mY9Aus2jrs9Drr+30F9wXmnMgmY9w+PVT+jbO2K3FTwhTFS1HH2eB9mAskxUsz4bBTAtk1ZXBWtUWCtoZqemivCakUZrk3FrFiYe6e4yG471JNDx1tQm6yhAnvbWV7GquTotZm1/cmaEJuUJi7B68ImhqaZvxgE1Jdh4UPSejJN84un+MhKDeKzK4c2ZS+xHBntKlYozBp1S+9wIL7Lnk7Q3Rk0r7EVBj+qCnx3dJovJT2LokCr6oJ6X9fYon0skgqOMUdlbRj8Zp9ZUhE70tMWo1leUNNsLPbDyyR2NfWYflPDMS4I3+bl+XqdpIWibzVN35cE3qChQNJsfsShZdFUgirpoaLm+4p7dVN3IeJrrH0psrdYA1x0VLfh01jPcKEzXj5TZWjo0qE6rcnbAxxMHRXvKpdKqEwYavPdfGSEEX9QL51YcKG2v0Ghl7PMv8YMSQWyDoW3Q8MFR1b43aawnLGuxiKrZZ4HsHQwVGYKBYm8xavpPqMsA5rIKu5HH0Xo6aRJUNigYOej1mgMAmLNDXGobABVlO8OtwszDxv2941DC9sNbZchsQrX5Z5IGzpgyhwKczydfTaJvRNZC3FqwNLGBg312HDj5elMFewWQ8GhQFDYWBQGAAMiWBQGDAUBgaFAcCQCAaFAUNhYFAYAAyJYFAYMBQGBoUBwJAIBoUBQ2FgUBgADIlYtAbcUdjqaQ9ppbDVQ2FgRK1JsxS2fliHIUNhYFAYAHPdhxFHRM0/7idsWmImuiURhYElojCwRBQGlojCwBJRGFgiCgNLxJtfNCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDIyJwsKgfCBf/U6R6PwYyH8NsHxioyKnOpEip7javSKnOpEiJ/NU3gnnlDFR2HHffydPZK6m/Cdtw+LkFDlViRQ5mZ9sLTaV51QnUuRkfiT29Pqgy6lkmrDL86H3TpHIPI5W/JO2p99/s1fmVCdS5ZQUm6pyKhMpcoqN10yVMifDNGHZU0r37XeKRKc37/LfIhZwef77kzanOpEmpyQpiogmpyqRMqd8U2VOhmnCTIHOPx31O0Ui8/xg6cmFu6I6UuRUJ9LklFYt+c+2anKqEqlySgPoTp1Tjl7YMQi6kX48ErcTCT+NJlG6adl+kOdUJ1LkZF4bm8rPKU+kyun8uGssUhWxG5r10xseingfZkP1qtOT5VQnUtZh85xTWiqbO59DmCn/l7++bb1TJDIBQd5KbEQ3YU51IkVOdUxTnlOeSJFT7Ut3ThnT78PS0J01arN3ykTq+zBlTnUiRU6znVMeAPbqc8pgTwcYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYMALGxhaFu6Gx5zV48rqIY9mWCEEngpLr79OWGrY9ZHdB2+EnR+zeTvpy2//eDBOyhX53J7L89+CYBdnw0uP32YD2s2m3+zLuT9mICgC3gg77rJhmeblxaExDDeb2/Pl28vTQ6rmIVt6fJkVP7NpsC/WT5n5swi+CDOK0mtvXi7PB1MjNUJiutRMOTD/jB+j5rjPNn2qJiWhxERfhJmXwlT6YspaJexoJjk2hX13MMKqYnXMJ0FS2DyMlbDzY1addUpYWrSKElasp7C5GKvDzMvp9aElrKjKsjqsWM86bC5MI8/Muapbib97zluJxYowMI3BlrBvsxkIl6eslZivZytxObIyo74PBomIvglLS00+f0d7/dnTQe4DhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBj/B/rwr8OHThAfAAAAAElFTkSuQmCC)
y_hat <- predict(GIC_base, Znew = Test$data$Comp, Zcnew = Test$data$Zc)
plot(Test$data$y, y_hat, xlab = "Observed response", ylab = "Predicted response")
![plot of chunk GIC_base](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAe1BMVEUAAAAAADoAAGYAOjoAOmYAOpAAZrY6AAA6ADo6AGY6OmY6OpA6kNtmAABmADpmAGZmOgBmZjpmZmZmtv+QOgCQOjqQOmaQZgCQkDqQtpCQ2/+2ZgC2Zjq2/7a2///bkDrb25Db/7bb/9vb////tmb/25D//7b//9v///+5QTkNAAAACXBIWXMAAAsSAAALEgHS3X78AAAIUklEQVR4nO2dDVviRhRG47ZKt4JuF2270G7aCjH//xc2MxkISCLzkYx523OedTWQO2NynA/CnVDUIEXx0b8AhIEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMRAmBsLEQJgYCBMDYWIgTAyEiYEwMVKEFTAFEwpLiIUhECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBgIEwNhYowjbP/5e3QsBJEqrFq5t2k+XShD2BQkt7Bq1aiihWVjhC6xWt3+hbBcjDKG7ReXHSLCpoFZohgIE2MsYdXjsVP0zO+BKGhhYiBMDD9hr0/F7d8Pm1GLhii8hL0+Lff3L7vblzGLhii8hDUzikbYybyig0tTmQloYWVvC2ueiysaovAew4pioEeshsY2hE0Bs0QxECaGn7Bm+CqLYj1q0RCF3yzxYdP8633TK75oiMJ3Wt+0MYTNAc8usbjZ7OgS5wCTDjEQJoafsN3A1aeUoiEKv0nHKnD48ikaovCdJY5fNETh1yVuhy7wJhQNUXh2iYxhc4FZohgIEyP9/bDYoiEK33eca3vJfsyiIYqAaX3o5B5hU5DcwvaL4mbTaxNhU5A6hr0+ra1PhGUieQWmFbW9Q1gmUoWZFtZQ/nD57ibCpiCgS7zr3aNa2ctWJYmkeQiZdPQbiy0aomBaL4ZnTsddfaWFsaAvEyFX6wOv1yNsCrj4KwbCxCBVW4zUVG0W9GUmOVWbBX15SU/VZkFfVph0iIEwMbhPhxjcp0OM1Pt0xBcNUSTfpyO6aIiCvEQxmCWKwXIjMfzGsOfAGb1P0RAFy43EYAwTI4+wNsHDK83D7HTcscsNGYw93bfdOFZ2HjOwdVHwRTrKeaGDJb53RGd7HjYuwwcKPI++UpnPb3Q11h6r+98j6LijPVHtt6HY033bL1dZa+/t7zAYdvJMURf1mzNU9B2B1/G0v0tRd+e8ePPVcywDh3jcere2eBB2KBBhl/sjrHs6LsntrGjGsIuT0Rc+3hi2JVV7LngJI1V7PngJu7ICc6C7RNgU+HWJw8uN3PqwuufdTYRNQeos8dBPsgIzE6nCaGGZ8e4Sh5Jwhi8MI2wKfCcdJOHMBN9pPUk4MyGghb2bhMMKzEwETOtJwpkDWS7+wnhwaUqMkKv1zBJnQEAL64UVmJlJHsNYgZkXP2Fm9WV505+cyArMrPh1idYJH0c1BzxfOJsrvFyamgN+XaKdWgSmdCBsEnjhLAbCxPB54fz4J4shZgMtTAyEiZEn8xdGg8xfMXh7RYzkzN/ooiGK1Mzf+KIhCmaJYiBMjOTM3+iiIQoyf8Ug81eMsTJ/w4uGKJIzf/eL4q5k9Uo2kteHPW/sNaumBQbHQgTJeYnNU+WSFZjZ8BvDhm+/RwvLjF8Le+cd52YMWzKG5YMrHWIgTAwfYebeKdc+O4wVmJnwEGbyfkPTtD2Khii80ty+R92mGWFT4CvsaahP5F5TeUkVxp1wMpOa5sa9pjKTfC2RFpaX5Ndh3GsqL7xwFgNhYiBMDISJgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxECYGwsRAmBhTCoMpmE5YDNH1EZhYeCTzPAtCgQgTC0SYWCDCxAIRJhaIMLFAXvyKgTAxECYGwsRAmBgIEwNhYiBMDISJgTAxcgmrVofPMeh+Cgwsi8BP6TzctDOwxi4wsMbuxq6BNXaB12vMJMzcYqz9iM3up8DAenvtLrdv2LlDD6yxCwyssXrY1PufNuE1doEeNWYSZm7e1/7ddj8FBobeM3V7862tJrDGLjCwxp1RZM94YI1doEeNmYSZjymo7OeTdT8FBtq7ygX9ybtzFlhjFxheY9wxdoEeNWYSZu622P5O3U+BgabTCPubd+c9sMYuMLzG9sNfI2p0gR416rQwS9CoktrCgmusVss6qsZDoEeNOmOYJUZY6BgWLWy/cPuG1ngM9Kgx2yxxeZwlLsNmiYfdTT/z+mvEeQ+ssT7rSwNq7E57YI1doEeNeV+HmTMR9TrMBDavUW6CBnITFFFjFxhWY2nzdtfhNZ4EXq+RKx1iIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDH0hJUuOdY3/+zwyeH/EeSElSar1uRdIkwCk6PZ5qFtv9qWtmvX51Qr823/8xeTX/b6vGm3zcM/fjHCmmc+fXcPtiHV4+82pcxF3v92WdwcURO2s/l+ptVsb18ab0agydnc2tTN/WJtUtP39y/ttnl4ZztQm6zZPuhCqtXty66xYnZpilrYb+fFzRE5Yce1YqZL3K5dCrv51pxsm4vZnOul2zbnv+0SbXqne9CFrNa2KT7aJ1zu6HlxH3WQ76EpzEj4Y2O97RemY7PrdG7a037/j+0RD9vtYOfycW1irQt5sAWYJxpvh2Tfs+I++Fh7URN2MoatD/MJ15W5x1+fv90ftt+0sK7NmJCHzWULOy9ujqgJq0szX7CzxDtzhs2YZr7cSGROeVksa7d9MobZtmYfdCHV6u5N5EVxH32sfcgJM/O49nXYV9tpbY/TOtPZGS12uXC73TSvwyzx83EnF1I9/HI6S3Qt7Ky4OaInbCwC17TOBYSJ8f8VJgrCxECYGAgTA2FiIEwMhImBMDEQJgbCxECYGAgTA2FiIEwMhImBMDH+BYtGjRAeIvR3AAAAAElFTkSuQmCC)
Regression with compositional predictors
Generate data
library(Compack)
p = 30
n = 50
beta = c(1, -0.8, 0.6, 0, 0, -1.5, -0.5, 1.2)
beta = c(beta, rep(0, times = p - length(beta)))
Comp_data = comp_Model(n = n, p = p, beta = beta, intercept = FALSE)
Comp_data2 = comp_Model(n = n, p = p, beta = Comp_data$beta, intercept = FALSE)
sparse log-contrast regression with compositonal predictors
m1 <- compCL(y = Comp_data$y, Z = Comp_data$X.comp,
Zc = Comp_data$Zc, intercept = Comp_data$intercept)
plot(m1, label = TRUE)
![plot of chunk compCL](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAABC1BMVEUAAAAAAAMAAAgAAAsAAA0AABcAACoAADoAAGYAAP8ABAAAFAAAGKEAHCAAIaEAKrwALdsAOmYAOpAAZpAAZrYAzQAA//8DOrwEOrwGZtsUZtsdkNsge9sikNsyAAA6AAA6ADo6AGY6OmY6OpA6ZrY6kNs6nP9J2/9Ytv9mAABmADpmAGZmOgBmOjpmOpBmZmZmtv982/+QOgCQOjqQOmaQZpCQkGaQnGaQtpCQ29uQ2/+d//+2ZgC2Zjq2/7a2///bkDrb/7bb////AAD/ADr/AGb/AP//OgD/Ojr/Omb/OpD/ZgD/Zrb/kDr/kLb/kNv/tmb/tpD/tv//25D/29v/2////7b//9v////SZhhMAAAACXBIWXMAAAsSAAALEgHS3X78AAANJUlEQVR4nO2di5/jNhGADaUsy5tle7y5HmVp91rosserS8qFzfX6Yrn0QhP9/38Jlq0kdizbM5JG1sjz/aDJJbJWzueRZFuSCyWwopi6AAIOEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMEcYMMmGb1Wr1Wqn1avUJLH2dbPfpa2DusGwdUj+sMMlxqber1RfwktigErZ98bXavPhaa3j4DyD92uz1wwoibFumegDvOS61+hJ0yDil1ocj6Ofoh0rYWher3pkN4AD88sVXVarN56AI23xRHxIwcKl3n0FT4lNvP4H9HANQtmHmR4IdUtV+7D773+TCylprBQpzh9S4klghFLb79JXS7QesgJWwh1ewNqyq5BAKMKnLihwRN7jU65SFbVevzBtQCbWw8gCEdzr+C6+McKk1RO1YyhG2OVYUoP2puie6ftl7HmOLagtwqamEJdyGGV/wOgDVrddZwntbuNS6yLvPaVLrViLRXmIdLK/1K6INA5+HrVHnM7jU4CI7pE73PEwgQoQxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4QxQ4Qxw0dYIVDgK2zzuMrmO88dthXw+Arb3d5Ur+uL+862xdmdY6mEXnyFbd9/3nptbntvsSh4Qhlh9+tLx1IJvXi3Ydvr3jasOEoENpnCKN7ChrbdPOnWk9huj9CGVJilYbOlE4EIQglruGn89jdORRJ7A5BGWADE3gmUwmA1IiZHMUcqzNZ1DMOMnVEK2/3F1q51gsTt159rrHlf6bg2v1w3mprCLDk31Vk6+/1YshkrZUZ4R9ju9qpv20Fh7gA0Zox/lbh92nOJl0qY7U/NR12kNiwOc9CW+nmYy989NIsTFYCUDIXpP12Yl/y05SnsoKx+m5OzXIW1PWXkLF9hquMsi/oxa2GqG1vsneUuTLWaM/PvVErmwgyEdZVxdjYLYRZlbM/TZiLMHlQclc1GmLL64adsTsKsReKmjPSOs9sgHEpseng1ZiGGauubmJZB2dBhbnGx2mHkLICwSsvm95ZtUxTWZ4eLsgDCNu/c2ydDLIrDzei3pqJnt6zKxvY1CfzHdJx9pO9TVtY62y77hg9Eo8+hvTGLXDgX/Dsdu9viUq3tE/qmF2bBaGN6YkbbS0x4upGWVhSdKjN5ZfM6D+tQKjtt6RJXRjkZwr1UEalK2e6hJF3wmUeYakXUUI8kFURYu6CHHslkpRnBW1gOyz509PSfwU2Or7ChSemORZoCS0Ql6izEtcTmK2bbpLDc4myfdU9RKAsSYQcsN6Wb/0pEm3cbNrDsg2uZpsI+jqD94eTapJfYxFpky1nlhNJgwpYX90v07UiGwoa686faJoo1kLDt07vyf5ZlUnyyTpSxH6SrjbY8nQJAvi67gGWMzUMYqNitWjKuM2CVWJzdredQJWrg5W46iyVNOh1dMAU/hlokZ9AqUVnPjX2yThj0lUQjLYYzgLDDyg7I5Sr5CnO6Xh/JGSLCAmedOC6lj+FstkMERnHbdXJnMGHrvstPI9sukhv6i8D1YN07C1mWRu6Qr7fXLj888wjzufFcdUJojNG2Ycu96FUDh6ymwq9FIAkyWJW4cBlfWB5lb/7a8vmqD4e/QY1Xx8kyis4fYJXo1IYVxRt/wGyQoETPnu4+yM7PzwMUps7S6+vhbYvie96djqkF+hrTLVllK5QxWmHf/eXxfZjhitHNeRvbB1cgYzBhu9vi4t99y+z1blsU3/i5ef9oT1k5FH2gco/mzdNY6emtw7sAgITtbq8276CfpKId/Kj5QSnr/FGX8z0uBum1eQkrw+vYwQ9iDNqtL4WhL/6W/Ni8f+nEwdx5l9afIm3j3I2ZUhb7s+gQxhARtnSIsF+Y971d+fo3PkZTK/Ys8g65WxxqKLS5NbuNg6rYKwvQWQS3YfCL9c1K7WfmsyFhcJznYXIHLcwJLeyHFoktOgU6oVXusZpxkDJKLZ8ecyoj+qQkJlLLL+B/pSfzqqSP3HM6Qizs++b9yM9+ymni034KvjBWW6pZBXcOgJX5HvgXRo8mf1kayA3M9//leqXjW38275u7YekonlYBY4IsWQymt+hq/eXjx+cnJQUcXVY7XcL4oo6w31g+H66jPRnz2DpGxrIIV6xgvoiFffNd897NTm/E9XPIfzBix7Sq4OcHweRDb2DeqCX2mbE6wn5l3g93gkZ/VvRu1dj7GeOE1xUuWKEjf8v/ogeS/rb49j8P/wIe2sF46WrL9DZCEnLPgCfO+qo7/tLUmz94z7wPbGOQ+kTbffukfQGrxOqGGPZpYEVx9vfoC6t4qapI2xfxyN9lvEE4L30Dy5C4L+JeYqwAC2GqIvzV49BtAeWJs3up4ASKK0P4i/3B226QMLdVDwmFHS/ih8yV4OZM+L4WLMLchrlR1InhNe2huJdG0DeGtGELc2UdXSXu/oo7EWjSc0/TOb9hSG59kpzLUEaY+nj82kjv3WaHP+gMzRgDmlNPSmHr40jtJLz0QKSL6FIBZZUICLDpIRrBQ3Zlx3ds/dAy6H9yb8FiQaWL7kqc74nz0DLo3Y8Sg2p8HOWFU9+BpEPLoAcoHh2EwxlJL3T7DiQdXgY9WQjHnhLfmPAfSDq0DHqi8NVFO5DUtUzEcNZFMpAUmPU0UA7Ej3KXNtTtFRadDtppE3Fuqs9m6SLqSS6xBkFQjq33K1lIqCeShR5BNAC006GqVS4tKdJfBp182l/MEUbei4MlvUhzlEmaUXX5R1iyy6DHmU8bW5d/G5ZghEWJq4qITdcB715iSsugx1OlprGlsujWtybfxiH4wHI4IGGLS5cFwhwujkA5nSMdlelkaSDC6t7G9rp/CJT1WZdFoOX3LJOdA+TqyLS2FGxMh7kRZrmDsqdPGHSo9uBUdFgWMZjclkINJB0YKNAj7DgKZ2R5AKeiRyYFWwokbKDnPrztgvf6lg0m7GScAmnDNm9X8/keYx80wHn92CbJyNLA10vETg+jGao9ASnpyuI8jJi0fImwMRLzJcJGSM2XCBsmOV8ibJD0fImwIRL0JcIGSNGXCOsnSV8irI+kLm80EGF2EtUlwuykGl5KhFlJV5cIs5GyLxHWIeHqUCPCTkhbF4mwBCdDgEk8vJREWIv0dYmwBhx0ibADPHSJMAMXXSKsgo8uEaZ46RJh6YwQBTJvYdxsqXkLY6hrxsLY1YWGeQrjakvNUhhjW2p+wnjbUjMTxt6WmoswiifKTcQMhGViypC1sHzi6kiewjKqAk/JTVi+pgy+woaeDOFcKBcyDqoWAYT1PhnCuVAIaJ4pnDABhE3wZAiiZz9zwFvYwJMhiiLgyip0j+jmhX+nY+jJEP6Lg4mjEyh7iagIi/G4+xyg7dYfVnPr0SFm0IQSZu90rHNYuigtIkWYEArKyRB5rA2WGLldmsoeEcYM/xPn3mcGFwIFnsLMYz4cwQUhKjVh1gmldtl8a336bLDsXVOn86OmJixe9iKMfPOw2YuwYJtb160Pl71L6nR+VBE2ddYJpXbdXIRNlDr05kJsRBgzRBgzRBgzRBgzRBgzRBgzRBgzRBgzRBgzqIXtn68JYY15SuDmcYF6pOPmCTTr7TXu8dTwnPGl7kItbAkvn97xJXQosb6rWj+6E8YafCzoQwxcDFTO+FJbIBa2+d0HqAMKfLRWDySGj+xfnP0DmrMeMYsIGkTO6FLboBW2e/YRokoswRzauLELYAd6mg5NzjU+Iy4UtbDlFaYNK2v4M8TO4EYHwWP3glSY15gmRShsURSX5cEKFLYwE2FAv1SdeHsN23OTdSIRBi11L6QRtqzG2aFKCK7gsY8JB/+syDYM20v0nYyQULceVReh9xz8s+o6C9WUIoT5+0pJmA5IcBtWBy+mfUzhPAxd6i5ypYMZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZIowZeQizjYPpHxtzHEqHGuaaBiKMGfkI215X00jKl5/88U4dhdUzfHbP/lYUV+tqWOviw2qomU75wU2AGUBRyUfY4qoaiqpfzhrCqhk+T57vbi9LNZd1yov7/QbFjfl+4h2Ak42w6rFZT+/0y+5ZM8JUNfpaf6b/r/3oKnFxU6Wsq0TLmvzJko0w7ac0oqcynAhb6ImdTWEf32lh1feVO8zEz8nJRlhfhG2vb6oqsR1hZWiZCDPfT7wDcLIR1teG6ZfN23ctYZfq2IaZ7yfeATj5CDv2En9aR9hjPfHgUk9AKDuDLWEfVpMudrdVL7H+fuo9AJOHsDacajg0uQkrwwY+Z4kjuQnLHhHGDBHGDBHGDBHGDBHGDBHGDBHGDBHGDBHGDBHGDBHGDBHGDBHGDBHGDBHGjP8DgdshRl9OWVEAAAAASUVORK5CYII=)
coef(m1)[1:10, 90:100]
## L90 L91 L92 L93 L94 L95
## Z1 0.92406820 0.92547909 0.92393318 0.92468211 0.92482516 0.92499664
## Z2 -0.66133490 -0.66328703 -0.66439024 -0.66621092 -0.66758205 -0.66890518
## Z3 0.58299125 0.58434121 0.58253225 0.58270957 0.58266568 0.58257745
## Z4 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## Z5 0.04906269 0.04925977 0.05387758 0.05507937 0.05690799 0.05856053
## Z6 -1.44762276 -1.44781669 -1.45108388 -1.45185109 -1.45315509 -1.45438880
## Z7 -0.29034027 -0.29063399 -0.29528538 -0.29726352 -0.29908001 -0.30087369
## Z8 1.02544144 1.02813418 1.03215631 1.03528559 1.03798521 1.04059747
## Z9 -0.01258296 -0.01473562 -0.01643933 -0.01842473 -0.01998444 -0.02156272
## Z10 -0.11472875 -0.11498203 -0.11470839 -0.11519857 -0.11511888 -0.11536527
## L96 L97 L98 L99 L100
## Z1 0.92528297 0.92544175 0.92555603 0.92573786 0.92614382
## Z2 -0.67023889 -0.67159139 -0.67288982 -0.67415268 -0.67543773
## Z3 0.58278803 0.58302517 0.58317946 0.58337527 0.58363210
## Z4 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## Z5 0.05989629 0.06128825 0.06265640 0.06387050 0.06490137
## Z6 -1.45558553 -1.45673691 -1.45775937 -1.45876054 -1.46061136
## Z7 -0.30234935 -0.30400968 -0.30576860 -0.30729023 -0.30818420
## Z8 1.04290265 1.04505940 1.04710022 1.04896656 1.05050565
## Z9 -0.02305226 -0.02429200 -0.02541595 -0.02649356 -0.02785752
## Z10 -0.11561523 -0.11568306 -0.11564942 -0.11550258 -0.11480866
Cross-Validation tuned model
cvm1 <- cv.compCL(y = Comp_data$y, Z = Comp_data$X.comp,
Zc = Comp_data$Zc, intercept = Comp_data$intercept)
plot(cvm1, xlab = "-log")
![plot of chunk cv.compCL](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAt1BMVEUAAAAAADoAAGYAAP8AOjoAOmYAOpAAZpAAZrYyzTI6AAA6ADo6AGY6OmY6OpA6kNtmAABmADpmAGZmOgBmOjpmZmZmtv+QOgCQOjqQkGaQtpCQ29uQ2/+pqam2ZgC2Zjq2///bkDrb25Db/7bb////AAD/pQD/pTr/pWb/tAD/tDr/tGb/tJD/tmb/xAD/xLb/0zr/07b/09v/25D/4Wb/4ZD/4f//8ZD/8dv/8f///7b//9v///9s1d3dAAAACXBIWXMAAAsSAAALEgHS3X78AAANjElEQVR4nO2dDXvbthGAlabzZDernbW1vU1Src1R3DTKsiiLVon//3eNAL8A4igCBEDckfc+6aMaJCGar+9wpEBqkTGkWKTeAcYNFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkYMFkaMIMJOu93H/iawLTv/9m1Q0+fdbvfSXmn3/g+96WisBDWBbdnxRX0xlnyF+nnpfM/dN3CTcgmw6x2EECaO5ud/9zWBbeLAG3asmj4Z64imr/pfxCnf7HNvE9iWH90X5cVYIty0fxu5AHrPXMbx/R/AJuUSYNe7CCHs9GL8HQJNYFt2/N2IHaum8wfjD/IEhMLH4phcbgLbPr3/70vzAizJjN+mWAB09lVYKv7CWpuUS4Bd7yKEMI/Dcv7wP1OFTVOeXnetqDt+/NTOKx571psSMzPC5AK4s6rJTDFiCbDrXYQQ9hX4mwL2Gmr7/MUYneyaRB5pRdkxF3gE8p/+lkAT2NYv7Gge4+NLV2fn377Am8glwK53kTbC8qa2CrumAn0cA99zt/vPh94msM0iwoxAqooOo7PT7kvHJnJJR1RCpB3DRKm3q34Xl6YCXdipIxeZBxwaMsy2fmFG5VObbG1zbLJ3a5NiSceuQ4SpEr+YVaLRBLYNLutFfj3/rrd9MlNifhCAuqy9F1BbnzAowRcp0eys9AVsUpk0d70LyudhZnYxSvCv5lsCTWBbX4QBO1CfoLU6K1LEN2CTagmw6x3wlQ5isDBisDBisDBisDBisDBisDBisDBisDBisDBisDBisDBisDBisDBisDBisDBisDBi+AhbMDGIKMxjW6YLFubDa7PpXeS3ZGE+sDCmjzGErYfh8dYTZqQIW8t/5cva4mca0qabEp2FFT9gh4URE5YAb2GHG3k6993TxW0HCSORFcfGV9j54V6+7q+eL207MMKw2yKYEk9/fdJeO7ZlYaHgCCOG9xh2uos2hrEwANxVIva6g2BKhLYBPggYHmF4bWVTEQZtO1FhCWBhxPAu6+/K/GdWHTMQhjUlAidZNeeHW4uuWVgorISdH1fdq5zediycgbAE2EVY57mWbdcewpBX9mNDoejAawtrSvTvmoWFwk7Y+SHPiObVQvuuJyosAXZFh6wEt47GWFgMHMr6S8V9X9cTFYY1JXpH2PW1+Fe+1P8u/iz/0wpFhN6wCvMdw64HQiHIxmaUKpGFhcP70pRV1wNTopIVB7z/CGBNiRcvTVl1Pajo0IMMI1iFJbo0hV9YAlBf6UCfEhMwzhjmfR6GtLLHmhITjWF6kGGMMqzC0n28ogxjGIUlAPUYxsJMkAvjlNjGQpgoOQ4/Pie9+MvCKryFHW4Wy23EqdrKzygLxbHxFSYKyO0ykytc6NpHWJkVq+aZ4ytMNG5vtVO1wFO1q7oDozCCKTF+hLEwDRth3ZN7MzmG3cYdwzglqmAv642f5w4NYSLIMArDmRJDdO0nTA5j8v+QVfYsrE8YviAbGxLC1JTIwnoXX64S7boOV3SgEoY2JW7kvMTl8K5ZWCiIzPxdqxNLZ80YM38HPi+xqAebYaz8ZAxZoTg2DjN/HTOiV8FSAgnDFGRYU2KUrm2AUyILG7w42raCVlbUB7jZYp0Sr/7VdfP5wK4tYWE6tkXH4cdn4IFtPl1bogpDVyhiTYl5QZ8L8ynrh6MIU6dQ4QCrsCLCfG6ZHYg6jCnC5lzaj3JDnydASpynLAHiKrGm47IIArCmRP+bIbxgYQoj3QzhhSoMXaE4NiPdDDGcpsBYoywUx4bCGJY1EYZMGNaUGKVrJ4yUiKOyRytse2le4r0o+oGaP4ow/VLVDLEbw96u9kv4E2fxRQOb+96Zv76owpSJpTPEtqwv/pkr5I2yhuyZW++F/mmmOnU7NVhTYq4k/3d4A6XEPLz2t1m2N8MvZIQJWJjAbgzLXe0XC/hpzBsZTUC6jCSsSYko6o6xIVIlCuAZO3ODhfmANSWmvNJR0JpFheXmCKzCCrZd3ygwrOshqHUHBmEJcBCW5hNnDUPY/OoOB2H7hCmxBEiJSW1hTYnlGHYftGtnOqe9pQOrsChdD6TWNNNhjKywuc61d0mJjoX9SMJSBhnalCgv1PvcHxYOICWysPZi//vDwgEUHfPJh5nb/WHpI0yfWIrs4+dxwH1/GEwlLP0whjUlRul6OCxs8OJo217ESIlzyoo2wk53S/EMMMcrU7GKDm2eIvJHzEbARtjmVs61QVB01EDzFBMEGc6UKJ+X+OYJSVlfAAobP8jwChP3hsGTcAZ37YWaEtXnKU5/JLNKiffyRGyDMCVmxhNLJ20rsy06FlfPovII2vVwuucpjh1kOFNipK69WetZMcWlKrLCwOFtLGHG8xQnfVJmLaxjCs6FT17SCZuyM19hubHvnvQICz63HgCY9tZ2NsZohjkldk9yO91d/TNFSsw0Qea30o5QgWAWdonDDXTZanxh+meaCarGEbA5cQ7wCNlIqBFVf/mz+QHnlKTZR9g29TS3Nsa0N/3ruvUiJM4uIE6JsrYI2nUg1BSoBdkYVSNeYduOm8M8ug5Ea8xqnBkRNqOUOCC8+rsOQWuOR31Sdt1cE56jsP2A8OrvOhytMetaIfJ5Gc6UiLhKLDCLDBlgrVCLUXvgFBap60C0pgyYZT7mL4wbAHlhBVAZryXHbCpn0S7nYWhTYuuMTK8aMz0zzmAME+wXi1ersF0HBxSmFo7hHyuGVNjhJre1cbzOMXJKNIMMLELoD2OWk3DkE2/Cdh0c41IVVISQH8Ys53Qs7vELy+ATaf1aY9isiDQlCjb4x7CCzjI/0+bEhQGxMHEHC9Yqsc2FqrE+KSPLRM7D2nQVIeowRtOb7VTtrq+2H951dNbrtTmmhc2KOFMiVWEFxpgWMiuysCjo+ZF4VpyDsEx3pswUJshMhJWoQRZCGNKUiP3zMEuaIJu4sMvIC43gQzywCSuosmL+H8lhzFeYuJdW3D1GQ5g2ktE8hfYVVl4ZXsZ8bn1oFGGeQYYzJcqnqnSNYSLCcrbfm7PrcQsLcgqNU5j45pXuaVOnO7lsO/7tRsNQTqRJZkXriaSXpA3qOim1MHKlh8MYtiFe1tfUZ9HKlKphIE2JWfF9Ha6fYKIVVrJuBdkAb0iFbd1lWXSdllaQ0RnKfKvE4V1joL7sQWcom+gHmJaoU6qGBBnSlBin6+S0LnsMCTIWNjK1M+96cTTmLayA1EnZ7IW1TsrcnHFKTELlTL3VxQ4WlhLl9qS1DLTUOwTCwgraibFMjql3S+6YticsrKZxZh1q0VJi8c7VPzVJszCVVpnfhFrHuBZSWC1ofW3SrMbCWqyNEkQNt7DvVV4TK6Ooi0x7WxYGUzrLzOO3Lp/152qv2MTO0bW2snaewcI6aUKtnSJrXlfmmgP8rvVztdyCZuV1C2WvWFgPSoo0g+K1edTfWblRHK2roavTkQYLs6bRZpPSANZqXaE6ctkLFjYMxV4j4OLPga5SsjAf+NIUMSgKO9x0TR+YgbAEhJhbL9hfPTtvywwgzNx67e4V9HPrw0EwJV6KsOnz2mx6F/s9PYVVt/sNeMasUwS6hes81o7TQ6C+8RwmPGvH6SFQ33gOE5614/QQqG88hwnP2nF6CNQ3nsOEZ+04PQTqG89hwrN2nB4C9Y3nMOFZO04PzKiwMGKwMGKwMGKwMGKwMGKwMGKwMGKwMGKwMGJEFCZm7zg8r+UAfRk7yOluYX78HaRn153eO36wW31870E8Yae3q+zww8p29b31ry5+6+3Sfkfse3bdafGX4LIn7l+VbRJP2F78JtbfwLN59XfbOJDPlLYPGoeeXXda4BS+f/kJsTCB+IO1xfo3lw8Ad+jY6ZhmbjuduUTY+fEfmFNilslnBVtjfVjF/K2Iwtx2+uaV/Y5sb9GOYZvFYlk/y9RuZSwRZrnTzfrWe5LvN1phksON095ZH1bHMcy1SnQ9pNYj3lZOFxz0FecK8YS5/urWh1XkLKfazEGY2047J2fUEVb8QdnvIIrzMMed3jp+Ny9qYUwUWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxWBgxJiCsY47N9rZ78k0zNa2ZFiOmOxJgssLy4+8mLDcces9iMCFhpzt5m0r+8qefV8JJtaC4hej8+LfF4nYvp3JufpFz2cSqP91Xtxgpz1xFzISEbW7lzE7x8mqlTA+WtxC9eTo/LHM1S9m6uZLhJ1Zd3JfLHe9aScV0hAlF+bEXL+fHlRiRlJSYt+aNYoH0J9Rs7uWqRUqU0UUiJxIXJm6kKL2Il9JU/iJirRa2ETdKqsJ+XQlhdVhtihspWdg49EXY6U4OZ60Iy0OrjLByOQsbi74xTLwcflhpwsqhTI5h5XIew8ai+OqKpVIl/vmxqBLLBduFKAY1Yb/ImxjOD7JKLJZzlZgOGTPO58EkMuL0hOVRU9wC5Hr8+UoHEwMWRgwWRgwWRgwWRgwWRgwWRgwWRgwWRgwWRgwWRgwWRgwWRgwWRgwWRgwWRoz/A17LCG0lZOlhAAAAAElFTkSuQmCC)
beta_est <- coef(cvm1, s = "lam.min")
head(beta_est, 10)
## 1
## Z1 0.90048789
## Z2 -0.62012673
## Z3 0.58372632
## Z4 0.00000000
## Z5 0.01756609
## Z6 -1.38814631
## Z7 -0.23387376
## Z8 0.94222496
## Z9 0.00000000
## Z10 -0.09116647
sum(beta_est[1:p]) # satisfies zero-sum constraint
## [1] 5.516289e-09
y_hat <- predict(cvm1, Comp_data2$X.comp, Comp_data2$Zc, s = "lam.min")
GIC tuned model
GICm1 <- GIC.compCL(y = Comp_data$y, Z = Comp_data$X.comp,
Zc = Comp_data$Zc, intercept = Comp_data$intercept)
plot(GICm1, xlab = "-log")
![plot of chunk GIC.compCL](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbAAAAEgCAMAAADi9tZbAAAAqFBMVEUAAAAAADoAAGYAOjoAOmYAOpAAZpAAZrY6AAA6ADo6AGY6OmY6OpA6kNtmAABmADpmAGZmOgBmOjpmZmZmtv+QOgCQOjqQkGaQtpCQ29uQ2/+2ZgC2Zjq2///bkDrb////AAD/pQD/pTr/pWb/tAD/tDr/tGb/tJD/tmb/xAD/xLb/0zr/07b/09v/25D/4Wb/4ZD/4f//8ZD/8dv/8f///7b//9v///8cH0dIAAAACXBIWXMAAAsSAAALEgHS3X78AAALWklEQVR4nO2dDXubthqGyb7S9izZ2VnSndmrGzfrmq6rs3qz+f//bIgP8yFhJJDQ++DnvrY6kZCQua1XQgaSpASKJHYDiBsUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgeBF22G4/DCcZ09Ljb19GJX3abreP3Y227/5pJ+21jUxJxrR0/9h80XKeTfU89u5z+8VYpMwxNL0HH8LU0fz0x1CSMU0deM2OVdJHbRuV9Nz+RByyYp8Gk4xp2dF9bLxoOcpN993kGaZ9ZjL27/4xFClzDE3vw4eww6P2OTQkGdPS/e9a37FKOr7XPpAHQ1f4UByT80nGtI/v/nqsXww5qfZuigxDZc/KUvEJ6xQpcwxN78OHsAmH5fj+b12FTVIWXredXrf/8LEbVya0bDAkpnoPyzPMlVVJeohROYam9+FD2LPhM2VotSnt02dtdLJLUnGk08v2mcC9If61d2lIMqYNC9vrx3j/2FfZ8bfP5iJ5jqHpfcTtYVlSV4VdUkF7HDPuc7v98/1gkjHNoodpHamadGiVHbafe4rkOT290kTcMUxN9bbVe3FJKmgLO/TEIv2Am4YMPW1YmDbzOZnslNnX0btTpMjpaboJP7PEz/osUUsypo2e1qv4evy9nfZRD4nZQTDMy7qtMKUNCTMF+CIk6pWVvgxFKpN60/tAPg/To4s2BX/Wd2lIMqYN9TBDA04naJ3KihDxxVCkyjE0vQeudIBBYWBQGBgUBgaFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYFBYWBQGBgUBgaFgTFGWEJC4l+Yj7KkDwpz5yHmzqULG4wBEaCw83XIMxYTCgODwtxhSDxXid1sdlYobKgapUuSsZhQGBgU5g5DYm8N+dBVjmFyxjEKO1NBYvzlgqEwMCjMnQsIifUw1Pqi4NzvVZlmHRMa45PlC6snermI8r+zvxt7lBxn8UASJigqxoPC3Fl+SBw7hhkqpLBJ2Rm7F/nx/erNiLKuiFsGjsBUYcf7u/z16du3zmVHIKWTxWOqsMN/37ReXcqOQIQw6JB4iT0MWlh6uJ1vDEt5KiZ6aYoYoDB3sENiSWPSIfEyDK8sQpjnsqQPCgPDx0rH1Wq287C83tjBFjokqvOw4/3NjMLin4pBCytEra8vSVhUPK10bL5+RWGz4GGl40a9bPSljoBjWNxhDDokBio7XHNEYxQ2puZLjYq4wqLP7uMAKayyFckYQ+LY2inMNTtYWbvaLzEmwgqLP7uPA66weJ2MIXHCDijMKTtYWesdXFxURBaWRp3cRwJaWKROxpA4aR/z37FOYVN2cmlREV3YxU098IXN38kYEqfSXr0P3t0obCpJ4/+lh8hFCDt93VLfeNve+4IEBhAW51Ltk6fqh0YDPI9xDIleOHWyehJSOaMwH1X7pw6GhbSyoy3ru7MlCWtN8AthVXAU9TS4SSxKWKsvVZ2s6nT+WsSQ6BH98VQU5qnqeThFyaWMY0sXlp463UKuAbkAYSd8rYEwJM6Ery/PKGwmfAmLyiUJW8S3nRclLOUYFqqsaOQLW1+rOy3v/FYdDeylKhthm/xJbeW9sd6qjkVj6iG0hWexEHb4YZX/vPtee8LelKpjMV2Y9JB45hmWU6qOxfKFnXmG5ZSqo7H8MWz3UsXE3QvjrGPOhzQTy1li/tRRg5J09kfI+mVkJ5MeEs8y80OavTJ2GIMWhtzDIGf2NrPE2/KyNWNQnPkhzT5ZqrBAVQvgEsewQGVFI1/Y5jrd9E0TS1Af0gzSzBrbtcTNjdLms2oR4A1jtmuJmbCd/gzLKVWLYJQw6SFxYWuJLRYpbGlriS0WOYY9JcrY09XKsMGZkzSwIwGC1SzxeN+7lpg/A31U1WJw7mTSQ+IA1febY8pKwH0YAxcWpOyMgM3sKUy2MO2rVgoTPYYlafcTRWE5TtIoLDpiwyKFmRErjGOYGTdhM4TE/m88KCxH2BiWpI3b6bWsgZJT9oqEpFVFChtG1DBGYcM4COMYJgFZwvqhsApJY9gZKKwBgjMKq7GNigyJQqAwMOLO7C0v56SwBjHHsMap19lmUFgbG2d+QmJnT+1z5f5GBBAGdal2B6uo6EVYd08RhfkoG4v5hrHGnpqPl248Yrq/3EC1U5qERgxhpa1GFscwe2yCuecxrAiGtoeLwroMdzIPwpofCgqbxhxRsb0Pp6cCUliXGYXV8dd+Tk1hGjOMYYUwt1jYKDo+O1jZ2JyV5msMozBvzDO7pzBvzHQ6NmZBiMJMnBfm5zxsJBRmJPQYNh4Km5fJa+IUNivauuGoGsZnBysrg55DOyEkjpoYdmuYkB2srAj6ph4UJpIwf/djcq0U1sP0vtBbM8ewEJyZHmBP65f6VO16QV3LghaG/Mzf8yTDF8TEYKow5KdqW7A8YcvtYTlGYdAhEfmp2jaYFtRHC/NxrSZniUP4i4peaqKwIRYqDPWp2sPoh3lsSBQlzHNZSWgPoxkhzN+HmMJsmNo3PK5zcaXDhgUJW/h5WEl7YdE9JAoStvCVjorWtyKOwsox0NMsjD3Mkp6oOKzB8+oWVzosMX/dYmFDmrAwZQVSR8VGSGxrNHY3CouGWVg9m9DUnO6F9dyI8dnByoqkIaR5W3L5gx41/dtKKcyJ1j3jSevf8qfyORvl5v4m8402TMoOVlYquZmHtqqk3d3qZ6NQWHwKYY2h6TR6FWoajiwe4TC2BeOzg5WVSqsHafPDpPVD43e/LZiQHaysWGohafsdNrpbHRlDNGBSdrCygknSB5Mw07NRKEwE2RiWvwwEu0Df4FIYGBQGBoW58zC8SThCClsqD1H3Hk7YoNBgG1/I1mFq8FS3nMMkZ+swNXiqW85hkrN1mBo81S3nMMnZOkwNnuqWc5jkbB2mBk91yzlMcrYOU4OnuuUcJjlbh6mBzAqFgUFhYFAYGBQGBoWBQWFgUBgYFAYGhYERUJi6OfrOYfNX+h1oZg63iX53oZeaXRv9ZLxvrp/q7sgJhBN2+GGV7l6ubDd/sn7r6l1vru0bYl+za6PVJ8GlJenG5RNsJpywJ/VO1rYNXF/937YfqLutHTqNQ82ujVY4dd///ChYmEJ9YG2xfue77986Vex0TFO3RqcuPez4+hfJITFVwevGfmPrw6pujw8ozK3RL67sG7K5ETuGrZPkWk0OrN56sbGUHmbZ6Hp765Zk7RYrLGf3wql11ofVcQxznSW6HlLrEW+TX3Xo9nHQCSfM9a1bH1YVs5zmZg7C3BrtHJxF97DiA2XfQBHnYY6NzjZ3GMOECyNBoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA2MBwnqusdnc9F98U1+aVl8Woy53BGCxwrLj7yYsM+y7ZSFYkLDDbX6bSvbyzf9WykmVUdxCdHz9c5LcPOWXcq5/yq9lU5v+eFfdYmT4k3YCWZCw9U1+Zad6uVo1Lg/ObyF69eZ4f52puc5T19/m3U9tmtyV+Y53rcRiOcKUouzYq5fj65UakRohMUvNElVG7k+pWd/lmxYhMe9dEDERXJi6kaL0ol5KU9mL6msnYWt1o2RT2K8rJezUrdbFjZQUNg9DPexwmw9nnR6Wda2yh5X5FDYXQ2OYetm9XLWElUNZPoaV+RzD5qL4y+DXjVnid6+LWWKZsUnUZLAl7Kf8JobjfT5LLPI5S4xH3mecz4MhIuLyhGW9prgFyPX4c6WDhIDCwKAwMCgMDAoDg8LAoDAwKAwMCgODwsCgMDAoDAwKA4PCwKAwMCgMjH8Bnwnu1o1naM4AAAAASUVORK5CYII=)
beta_est <- coef(GICm1, s = "lam.min")
head(beta_est, 10)
## 1
## Z1 0.85898896
## Z2 -0.56310354
## Z3 0.58218204
## Z4 0.00000000
## Z5 0.00000000
## Z6 -1.32214479
## Z7 -0.14017233
## Z8 0.82199822
## Z9 0.00000000
## Z10 -0.08630774
sum(beta_est[1:p]) # satisfies zero-sum constraint
## [1] 6.958372e-09
y_hat <- predict(GICm1, Comp_data2$X.comp, Comp_data2$Zc, s = "lam.min")