function [optRes,pvalue,gppRes] = Ex4gPOCRE(t)
%function [optRes,pvalue,gppRes] = Ex4gPOCRE(t)

% An example to show how to use gPOCRE to analyze simglmdata.txt with the
% tuning parameter chosen by EBIC.
%   Step 1. Screen variables using gPOCREScreen();
%   Step 2. Fit with gPOCREPath();
%   Step 3. Select the best model with SelectModel();
%   Step 4. Calculate p-values using RBootsGPOCRE().
%
% Input:
%   t -- t*n variables will be selected at the screening stage.
%
% Output:
%   optRes: results of the selected model;
%   pvalue: p-values from Boostrap;
%   gppRes: results of gPOCREPath();
%
% Dabao Zhang     August 19, 2016
%
% Revised by:
%   Dabao Zhang, 12/18/2017 -- Input & indices of screened variables. 
%
% Copyright (c) 2013-2017 by Zhang.

if( nargin<1 ) t=1; end

YX = load('simglmdata.txt');
[n,p] = size(YX);
p = p-1;

% Settings:
%option.covidx = [15 16];    % covariates of no interest
option.delta = 0.01;        % 0.01, 0.02, 0.05, or 0.1

%!!!!! You may want to standardize both of your X & Y before proceed!!!!!
%[X,idxX] = Standardize(YX(:,2:end));    % Standardizing all covariates
%Y = ...

% Step 1. Screen variables using gPOCREScreen();
option.maxvar = t*n;
option.maxcmp = 5;
%[tmpX,tmpXIdx] = gPOCREScreen(YX(:,1),YX(:,2:end),'binomial','logit',[15 16],option);
[tmpX,tmpXIdx] = gPOCREScreen(YX(:,1),YX(:,2:end),'binomial','logit',[],option);

% Step 2. Fit with gPOCREPath();
option.maxvar = min(size(tmpX))/2;  % or n/log(n); % Should be smaller than previous one.
option.maxcmp = 10;
gppRes = gPOCREPath(YX(:,1),tmpX,'binomial','logit',option);

% Step 3. Select the best model with SelectModel();
%optRes = SelectModel(gppRes);
%optRes.XIdx = tmpXIdx;

% Or the first sparse model with the largest R^2
bSparse = [gppRes(:).bSparse];
rsq = [gppRes(:).rsq];
idxOpt = find((rsq==max(rsq(bSparse==1)))&(bSparse==1));
idxOpt = idxOpt(end);  optRes = gppRes(idxOpt);
optRes.XIdx = tmpXIdx;

% Step 4. Calculate p-values using RBootsPOCRE().
rbRes = RBootsGPOCRE(YX(:,1),tmpX,optRes,'binomial','logit',100,option);
pvalue = ones(p,1);
pvalue(tmpXIdx) = rbRes.pvalue;

%-------------------------------------------------------------
% Visualize Some Results
%-------------------------------------------------------------

% Plot lambda vs. (R^2, #[beta~=0])
plotrsqnzbeta(gppRes);

% Plot lambda vs. (beta, R^2)
plotbetarsq(gppRes);

% Plot lambda vs. (beta, #[beta~=0])
plotbetanzbeta(gppRes);

% Plot all components:
plotcomponents(optRes);
