%%%
%%% This is the beginning of the actual thesis.  If you don't know latex
%%% then start with the LaTeX manual by Lamport and another easy
%%% reference, like the paperback by Jane Hahn, LaTeX for Everyone, PTI,
%%% 1991. See also $TEX/latex/sample.tex and $TEX/doc/story.tex, where
%%% $TEX==/usr/local/lib/tex
%%%
% Start this dissertation....
%
\chapter{Simulation Study} \label{simulation}   % level 1

\section{Simulation 1: Main Effects Model} % level 2

\subsection{Example 1}
In this example, there are altogether $D = 10$ covariates: $X_1, \dotsc X_{10}$. They
are taken to be uniformly distributed in $[0,1]$ independently. The sample size $n = 1000$. 
We use the simple random subsampling technique to select $N = 50$ basis functions. The 
perturbation $\epsilon$ is distributed as Bernoulli($0.5$) taking two values $\{+0.25,-0.25\}$. 
The true conditional logit function is
\begin{equation}
f(x) = \frac{4}{3} x_1 + \pi sin(\pi x_3) + 8 x_6^5 + \frac{2}{(e-1)} e^{x_8} -5
\label{eq:truef}
\end{equation}
Four variables $X_1$,$X_3$, $X_6$ and $X_8$ are important, and the others are noise variables.
We fit the main effects LBP model and search the parameters $(\lambda_\pi,\lambda_s)$ globally.
Since the true $f$ is known, both $CKL(\lambda)$ and $ranGACV(\lambda)$ are available for 
choosing the $\lambda$'s. 
\begin{figure}[H]
  \begin{center}
   \includegraphics[width=5.7in,height=6.5in]{sim1/gg.eps}
     \caption{Contours and three-dimensional plots for CKL($\lambda$) and GACV($\lambda$).}
   \label{fig:c11}
  \end{center}
\end{figure}

Figure \ref{fig:c11} depicts the values of $CKL(\lambda)$ and $ranGACV(\lambda)$ as functions 
of $(\lambda_\pi, \lambda_s)$ within the region of interest $[2^{-20},2^{-1}] \times [2^{-20},2^{-1}]$. 
In the top row, are the contours for $CKL(\lambda)$ and $ranGACV(\lambda)$, with the white cross ``x'' 
denoting the location of the optimal regularization parameter, $\hat{\lambda}_{CKL} = (2^{-17},2^{-15})$ 
and $\hat{\lambda}_{ranGACV} = (2^{-8},2^{-15})$. The bottom row shows their three-dimensional 
plots. In general $ranGACV(\lambda)$ approximates $CKL(\lambda)$ quite well globally.

In simulation studies, since the true model generating data is known, it is possible to compute
$CKL(\lambda)$, which is the true optimization criterion for choosing the smoothing parameter. 
However, in real life examples when the true model is not available, the LBP model will relies
on $ranGACV(\lambda)$ for tuning parameters. Thus it is vital for $ranGACV(\lambda)$ to be a 
``good'' approximate for $CKL(\lambda)$. By ``good'' we mean that two functions should have 
similar variation patterns in the parameter space and achieve their optimal values at same locations. 
It is not necessary for the two functions to have same range in terms of their functions values.

After choosing the optimal parameters $\hat{\lambda}$, we fit the main effects LBP model. The
direct solutions obtained are the estimates for all the coefficients. With a little extra work, 
it is not hard to get the estimates for all the functional components in the model. In order to 
select important variables or functional components, we compute their importance measure and make 
a thorough comparison. Here, both $L_1$ and $L_2$ norm scores are calculated for the individual 
component $\hat{f}_1, \dotsc \hat{f}_{10}$.

\begin{figure}[H]
  \begin{center}
   \includegraphics[width=5.5in,height=4.2in]{sim1/LL1.eps}
     \caption{$L_1$ and $L_2$ norm scores for the main effects model (Example 1).}
   \label{fig:c12}
  \end{center}
\end{figure}
 
Figure \ref{fig:c12} plots two sets of the $L_1$ and $L_2$ norm scores for the $10$ variables, each 
set obtained respectively using $\hat{\lambda}_{CKL}$ and $\hat{\lambda}_{ranGACV}$, in decreasing 
orders. Simple observation tells us that all the lines drop dramatically at the beginning and then level
off gradually. If a selection is made by using the ``eyeball'' method, we may either pick up the 
two variables corresponding to the highest two scores, or the four variables corresponding to the 
highest four scores. Different decision-makers may give different lists of important variables. To 
make an objective decision, we need a cut-off value decided by data-driven or model-driven method.
In this dissertation the proposed sequential Monte Carlo bootstrap test algorithm is used to compute 
the threshold value, which will help us differentiate important variables out of the rest. 

Let us focus on the left-hand sided plot first. The dashed line indicates the threshold $q = 0.21$, 
chosen by the proposed sequential Monte Carlo bootstrap test algorithm. By using this threshold, 
variables $X_6, X_3, X_1, X_8$ are selected as ``important'' variables correctly either by $CKL$ 
or by $GACV$. In the right-hand sided plot are shown the $L_2$ norm scores for each of the variables. 
Using the same threshold $q$, the same list of variables is selected.

The procedure of the bootstrap tests to determine the threshold $q$ is depicted in Figure \ref{fig:c13}. 
We fit the main effects model using $\hat{\lambda}_{ranGACV}$ and sequentially test the hypotheses 
$H_0: L_{(\eta)} = 0$ vs $H_1: L_{(\eta)} > 0, ~ \eta = 1, \dotsc 10$. In each plot, the variable being 
tested for importance is bracketed by a pair of *. Light color (green color in a colored plot) is used 
for the variables which are in the null model, and dark color (blue color in a colored plot) for those 
not being tested yet. The null hypotheses of the first four tests are all rejected at level $\alpha = 0.05$ 
based on their Monte Carlo $p$-value $1/51 \doteq 0.02$. However, the null hypothesis for next component 
$f_5$ is accepted with the $p$-value $10/51 \doteq 0.20$. Thus $f_3, f_1, f_6$ and $f_8$ are selected as 
``important'' components and $q = L_{(4)} = 0.34$. 

\begin{figure}[H]
  \begin{center}
   \includegraphics[width=5.5in,height=6.4in]{sim1/boot0-4.eps}
     \caption{Monte Carlo bootstrap tests (Example 1).}
   \label{fig:c13}
  \end{center}
\end{figure}

In addition to select important variables, LBP also produces functional estimates for the individual 
components in the model. Figure \ref{fig:c14} plots the true main effects $f_1,f_3,f_6$ and $f_8$ and 
their estimates obtained using $\hat{\lambda}_{ranGACV}$. 

\begin{figure}[H]
  \begin{center}
   \includegraphics[width=5.6in,height=4.5in]{sim1/uni_logit.ps}
     \caption{True and estimated univariate logit component.}
   \label{fig:c14}
  \end{center}
\end{figure}

The mathematical expressions of the main effects are: $f_1 = \frac{4}{3}x_1, ~ f_3 = \pi sin(\pi x_3), 
~ f_6 = 8x_6^5, ~ f_8 = \frac{2}{e-1}e^{x_8}$. In each panel, the solid line is the true curve and the 
dotted line is the corresponding estimate. In general, the fitted main effects model provides a 
reasonably good estimate for each important component.

Altogether we generated $20$ datasets from the true model and fitted the LBP model for each of the
datasets. Figure \ref{fig:c15} depicts the $L_1$ norm scores for the fitted models for all the 
datasets. The dashed line corresponds to the data set used above, which is the first data set we 
generated. The dotted lines are for the other datasets. In all the runs, variables $X_1, X_3, X_6$ 
and $X_8$ are ranked at the top, and $X_6$ always has the largest $L_1$ norm score, $X_3$ the second 
largest score. 
\begin{figure}[H]
  \begin{center}
   \includegraphics[width=5.5in,height=2.5in]{sim1/mult1.eps}
     \caption{$L_1$ norm scores of 20 simulated data sets (Example 1)}
   \label{fig:c15}
  \end{center}
\end{figure}
We observe that in 17 runs, $X_8$ has a higher $L_1$ norm score than $X_1$, while the other 3 runs
give $X_1$ higher score. The theoretical $L_1$ norms for the variables can be computed: 
$L_1(f_6) = 1.55, L_1(f_3) = 0.84, L_1(f_8) = 0.48, L_1(f_1) = 0.33$ and zeros for the rest variables. 
In theory, $X_8$ has a larger $L_1$ score than $X_1$, however their difference is not big. Thus the 
order of their $L_1$ norm estimates could be switched due to randomness in the simulated data. However 
we should note that this randomness in the data does not prevent the procedure from differentiating
important variables out of all the candidate variables. In the following is drawn the $L_1$ norm plot 
for the second data set, which belongs to one of the $17$ runs that rank $X_8$ higher than $X_1$. 
\begin{figure}[H]
  \begin{center}
   \includegraphics[width=5.5in,height=3in]{newsim1/L1L2.eps}
     \caption{$L_1$ and $L_2$ scores of the second simulated dataset (Example 1)}
   \label{fig:c16}
  \end{center}
\end{figure}

\subsection{Example 2}

The only difference between Example 2 and Example 1 is that in Example 2 the covariates $X_1,\dotsc, 
X_{10}$ have the truncated Normal distribution $N(0.5,0.3)$ on $[0,1]$. They are independently 
distributed. The sample size $n = 1000$ and we select $N = 50$ basis functions. The true logit 
function and perturbation $\epsilon$ are both same as in Example 1.

Firstly, we fit the main effects model using the original data and tune parameters by both $CKL(\lambda)$ 
and $GACV(\lambda)$. The optimal regularization parameters are $\hat{\lambda}_{CKL} = (2^{-9},2^{-15})$ 
and $\hat{\lambda}_{ranGACV} = (2^{-13},2^{-20})$. Figure \ref{fig:c17} plots two sets of the $L_1$ and 
$L_2$ norm scores for the $10$ variables, each set obtained respectively using $\hat{\lambda}_{CKL}$ and 
$\hat{\lambda}_{ranGACV}$, in decreasing orders. It shows that variables $X_6, X_3, X_8$ and $X_1$ have 
the largest scores.
\begin{figure}[H]
  \begin{center}
   \includegraphics[width=5.5in,height=3.5in]{sim1_normal/Srun2/L1L2.eps}
     \caption{$L_1$ and $L_2$ norm scores using the original data (Example 2)}
   \label{fig:c17}
  \end{center}
\end{figure}

Secondly, we transform each of the covariates by using its true distribution function, thus all 
the transformed covariates are uniformly distributed. Then we fit the main effects model using 
the transformed data and tune parameters. The optimal regularization parameters are $\hat{\lambda}_{CKL} = 
(2^{-10},2^{-14})$ and $\hat{\lambda}_{ranGACV} = (2^{-12},2^{-20})$. Figure \ref{fig:c18} plots 
two sets of the $L_1$ and $L_2$ norm scores for the $10$ variables. Again $X_6, X_3, X_8$ and 
$X_1$ are ranked at the top. Also we notice that the optimal parameters chosen using the transformed 
data are quite close to those chosen using the original data.  
\begin{figure}[H]
  \begin{center}
   \includegraphics[width=5.5in,height=3.5in]{sim1_normal/Frun/L1L2.eps}
     \caption{$L_1$ and $L_2$ norm scores using the transformed data (Example 2)}
   \label{fig:c18}
  \end{center}
\end{figure}

In this example, it shows that there are two ways to apply the LBP model to the data with Normally 
distributed covariates. One way is to handle the original data directly, and the other way is to 
first transform the data to uniform distribution and then apply the LBP model. Both of them select
the important variables successfully in this example.  
\section{Simulation 2: Two-factor Interaction Model}  % level 2

In this example, there are $d = 4$ continuous covariates, independently and uniformly distributed 
in $[0,1]$. The true model is a two-factor interaction model, and the important effects are $X_1$, 
$ X_2$ and $X_1 * X_2$. We generate $n = 1000$ samples and choose $N = 50$ basis functions. The 
distribution of the perturbation $\epsilon$ is the same as in Simulation 1. The true $f$ is 
\begin{equation}
f(x) =  4 x_1 + \pi sin(\pi x_1) + 6 x_2 - 8 x_2^3 + 8 x_1 * x_2 - 6 \nonumber
\end{equation}
There are five tuning parameters $(\lambda_\pi, \lambda_{\pi \pi}$, $\lambda_{s}$, $\lambda_{\pi s},
\lambda_{ss})$ for the two-factor interaction model. In practice, extra constraints may be added on 
the parameters for different needs. In this example we force all the two-factor interaction terms to 
have the same penalty parameter, or equivalently, we set $\lambda_{\pi \pi} = \lambda_{\pi s} = \lambda_{ss}$. 
The optimal parameters obtained are $\hat{\lambda}_{CKL} = (2^{-7},2^{-6},2^{-8},2^{-6},2^{-6})$ and 
$\hat{\lambda}_{GACV} = (2^{-8},2^{-6},2^{-8},2^{-6},2^{-6})$. Two sets of optimal parameters are 
pretty close. The ranked $L_1$ scores are plotted in Figure \ref{fig:c21}. . The dashed line in the
plot denotes the threshold $q$. This experiment shows that the LBP two-factor interaction model, fitted 
using either $\hat{\lambda}_{CKL}$ or $\hat{\lambda}_{GACV}$, selects all the important effects $X_1$, 
$ X_2$ and $X_1 * X_2$ correctly.
\begin{figure}[H]
  \begin{center}
   \includegraphics[height=3in,width=5.7in]{sim2/L1.eps}
     \caption{$L_1$ norm scores for the two-factor interaction model.}
   \label{fig:c21}
  \end{center}
\end{figure}


\section{Simulation 3: Incorporating Categorical Variables}  % level 3
In this example, there are both continuous covariates $X_1, \dotsc X_{10}$ and categorical 
covariates $Z_1, Z_2$. The continuous variables are uniformly distributed in $[0,1]$ and the 
categorical variables are Bernoulli($0.5$) distributed with values $\{0,1\}$. The true logit 
function is
\begin{center}
$f(x) = \frac{4}{3} x_1 + \pi sin(\pi x_3) + 8 x_6^5 + \frac{2}{(e-1)} e^{x_8} 
+ 4 z_1 - 7 $. \nonumber
\end{center}
The important main effects are $X_1, X_3, X_6, X_8, Z_1$. Sample size $n = 1000$ and basis 
size $N = 50$. We use the same perturbation $\epsilon$ as before. The model in $(\ref{eq:obj3})$ 
is fitted. The ranked $L_1$ norm scores are plotted in Figure \ref{fig:c31}. When using the 
threshold chosen by the bootstrap test procedure (denoted by the dashed line), the LBP models 
using $\hat{\lambda}_{CKL}$ and $\hat{\lambda}_{GACV}$ both select the important variables correctly.
\begin{figure}[H]
  \begin{center}
   \includegraphics[height=3in,width=5.7in]{sim3/L1.eps}
     \caption{$L_1$ norm scores for the model incorporating categorical variables.}
   \label{fig:c31}
  \end{center}
\end{figure}










