\documentclass[11pt]{article}

\author{Jeremy Morris}
\title{Multivariate Analysis \\ \small Project 1}

\usepackage{amssymb,amsmath,amsfonts}
\usepackage[dvips]{graphicx}
\addtolength{\hoffset}{-1.05cm}
\addtolength{\textwidth}{2.5cm}
\addtolength{\textheight}{2.5cm}

\begin{document}
\maketitle

\section{Generating bivariate normal data}\label{gen}
Definition 2.2 from our text states that we can transform a sample from a standard normal random variable ($\boldsymbol{Z}$) into a multivariate random variable with the distribution $\boldsymbol{X} \sim N_{m}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ using the equation
\begin{eqnarray}
\boldsymbol{X} & = & \boldsymbol{AZ} + \boldsymbol{\mu}
\end{eqnarray}
The matrix $\boldsymbol{A}$ is defined as $\boldsymbol{\Sigma = AA'}$ \cite{text1}.  Theorems A.4.3 and A.3.3 in the text suggest that, if $\boldsymbol{\Sigma}$ is positive definite, we use the spectral decomposition $\boldsymbol{\Sigma} = \boldsymbol{T\Lambda T'}$ and set $\boldsymbol{A} = \boldsymbol{T\Lambda}^{1/2}$.  It follows that we get the correct decomposition of $\boldsymbol{\Sigma}$.  The function \texttt{mvrnorm} in the MASS library of \texttt{R} uses this method.

\section{Kernel Density Estimators}
The bivariate kernel density estimator with kernel $K$ and bandwidth $h$ is defined by
\begin{eqnarray}
f(\boldsymbol{x}) & = & \frac{1}{nh^{2}}\sum_{i=1}^{n}K\left\{\frac{\boldsymbol{x} - X_{i}}{h}\right\}
\end{eqnarray}
If we assume the data to be bivariate normal, the kernel function $K(x)$ will be defined as
\begin{eqnarray}
K(\boldsymbol{x}) & = & \frac{1}{2\pi}\exp{\left(-\frac{1}{2}\boldsymbol{x}'\boldsymbol{x}\right)}
\end{eqnarray}
The kernel density estimator can be calculated using the function \texttt{bivden} provided by Everitt \cite{rsplus}.  

The parameter $h$ is known as the bandwidth of the estimator.  There is no closed form solution for an optimal bandwidth, but there are some suggestions.  Everitt's function \texttt{bivden} uses $h = an^{-0.2}$, where $a$ is some constant provided by the user.  This choice for $h$ may not be optimal since it does not use any characteristics of the data set.  H\"ardle suggests that if we can assume the data to be normally distributed, Silverman's rule of thumb can be used to choose the bandwidth \cite{hardle}.  Silverman's rule of thumb suggests taking an estimate from the sample based on the sample variance so that we have the estimate
\begin{eqnarray}
\hat{h}_{rot} & = & \left(\frac{4\hat{\sigma}^{5}}{3n}\right)^{1/5} \approx 1.06\hat{\sigma}n^{-1/5} \label{rot1}
\end{eqnarray}
H\"ardel further notes that this definition of $\hat{h}_{rot}$ is sensitive to outliers, so it is suggested to modify the definition of $\hat{h}_{rot}$ and take into account the 25\% and 75\% quantiles.  Then if we define
\begin{eqnarray}
R & = & X_{[0.75n]} - X_{[0.25n]}
\end{eqnarray}
where the $X_{i}$ are normally distributed, we can use
\begin{eqnarray}
\hat{h}_{rot} & = & 1.06\frac{R}{1.34}n^{-1/5} \approx 0.79\widehat{R} n^{-1/5} \label{rot2}
\end{eqnarray}
The two estimates in \ref{rot1} and \ref{rot2} can be combined so that we have the final version of Silverman's rule of thumb.
\begin{eqnarray}
\hat{h}_{rot} & = & 1.06 \min \left\{ \hat{\sigma},\frac{R}{1.34}\right\} n^{-1/5} \label{rot.final}
\end{eqnarray}
The function \texttt{kde2d} in the \texttt{MASS} library of \texttt{R} uses the $\hat{h}_{rot}$ by default, whereas the \texttt{bivden} function assumes that the user will provide a constant.  This small difference could change the smoothness in the plots if the proper constant is not provided to \texttt{bivden}.  For this reason, the function \texttt{kde2d} will be used so that no unwieldy calculations have to be done.

Figure \ref{indiv} shows contour plots for four different values of $\rho$.  These values are $\rho = (0,0.25,0.75,1)$.  The plot shows that as the correlation between the data samples goes up, the contour plots get more elongated along the line $y=x$.  This makes sense because a bivariate sample that is completely correlated will have the same observations in both elements.
\begin{figure}[ht]
\begin{center}
\includegraphics[scale=0.6]{contour.eps}
\caption{Contour Plots of the Kernel Density Estimator}\label{indiv}
\end{center}
\end{figure}
%\clearpage


\section{Testing for Multivariate Normality}
Johnson and Wichern suggest that a good test for multivariate normality is to look at Q-Q plots of the components of the sample \cite{text2}.  This is done by plotting the sample quantiles $x_{(j)}$ against the quantiles for the standard normal distribution $q_{(j)}$.  Sample quantiles can be thought of as the order statistics from the sample.  Equation \ref{qeq}
\begin{eqnarray}
P\{Z \leq q_{(j)}\} & = & \int_{-\infty}^{q_{(j)}} \frac{1}{\sqrt{2\pi}} e^{-z^{2}/2} dz = p_{(j)} = \frac{j-\frac{1}{2}}{n} \label{qeq}
\end{eqnarray}
can be used to find the standard normal quantiles, $q_{(j)}$, by using $p_{(j)}$ as calculated from the ordered observations.  If the points $(x_{(j)},q_{(j)})$ are approximately linear, we can be fairly certain that the observations do not violate the assumption of normality.
\begin{figure}[ht]
\begin{center}
\includegraphics[scale=0.6]{qqplot.eps}
\caption{Q-Q plots to test for normality}\label{qqplots}
\end{center}
\end{figure}

Figure \ref{qqplots} shows Normal Q-Q plot for two bivariate data sets, with the least squares line plotted in red.  The first sample, $\boldsymbol{X}$, was generated using the method described in Section \ref{gen}.  The Q-Q plots show that this data is normally distributed.  The second bivariate sample $\boldsymbol{Y}$ was generated using \texttt{R}'s exponential random number generator.  The Q-Q plots show this data to be far from normality.

\section{Visualizing Trivariate Data}\label{tri}
Everitt suggests using a scatterplot matrix or conditioning plots to visualize trivariate data.

\subsection{Scatterplot Matrix}
A scatterplot matrix is a square matrix of bivariate plots for each pair of components in the data set.  In the case of trivariate data, six plots will be created.  Everitt also suggests placing least square lines and lowess lines in order to better visualize the relationships between the different variables.

Figure \ref{scatter} shows a scatterplot matrix for data on population, income level and murder rates in all fifty states during the 1970's.  Also plotted are the least squares fit and lowess fit for each plot.  It is apparent that the scatterplot, in this case, does not fully capture the relationships between the different variables in this data set.
\begin{figure}[ht]
\begin{center}
\includegraphics[scale=0.6]{scatter.eps}
\caption{Scatterplot Matrix}\label{scatter}
\end{center}
\end{figure}

\subsection{Conditioning Plots}
Conditioning plots, or coplots, can provide more information than scatterplots.  To construct a coplot for trivariate data, we need to specify a conditioning variable.  We will use the same data set as in the previous section and choose to condition on population.  The \texttt{coplot} function is used to display six bivariate scatterplots, for the variables representing income level and murder rate, where the data are separated by population density.

In Figure \ref{coplot} we can see the top panel which gives six bars representing the density of the population variable.  Then each bivariate plot, from left to right, starting in the bottom row plots income level against murder rate for each of the six bars in the density plot.  The coplot gives a better picture of the relationship between murder rates and income level.  Notice that in all but one plot (bottom left) there appears to be a downward trend, suggesting that as income level rises, murder rates decline.  This relationship is not as clear when looking at the scatterplot in Figure \ref{scatter}.
\begin{figure}[ht]
\begin{center}
\includegraphics[scale=0.6]{coplot.eps}
\caption{Conditioning Plot}\label{coplot}
\end{center}
\end{figure}

\section{Application to a Bivariate Data Set}
From the data in Section \ref{tri} we will look at a bivariate data set consisting of the income levels and murder rates from the 1970's states data.

Figure \ref{mvqqplot} shows Q-Q plots for both components of the data set.  The income level data appears to be approximately normal.  Data for murder rates appears to be normally distributed, but the fit is not great.  This could be because the murder rates are integer values, and are thus not normally distributed.  A transformation could be used to improve the normality of the data.
\begin{figure}[ht]
\begin{center}
\includegraphics[scale=0.6]{mvqqplot.eps}
\caption{Testing for Normality}\label{mvqqplot}
\end{center}
\end{figure}


Figure \ref{kernel} shows a contour plot for the kernel density estimator of this data.  The contour plots shows some interesting clusters of data.  There is one large cluster around the income levels of 4500-5000 suggesting that there is a low murder rate in middle income areas.  There is another smaller cluster in the income range of 3500-4000 suggesting high murder rates in low income levels.  There is another small cluster around the 6000+ income level suggesting a high murder rate in a few upper income areas.  These observations would suggest that there may be at least one other factor that should be included in any analysis attempting to explain murder rates.

\begin{figure}[ht]
\begin{center}
\includegraphics[scale=0.6]{kernel.eps}
\caption{Contour Plots of the Kernel Density Estimator}\label{kernel}
\end{center}
\end{figure}


\clearpage
\input{bib.tex}

\end{document}
