\section{Gaussian Mixture Models and Glaucoma Diagnosis}\label{S:GlaucomaProject}

\newpar
Glaucoma is a condition caused by elevated pressure inside the eye due to a build-up of excess fluid. Left untreated, this pressure impairs vision by causing irreversible damage to the optic nerve, which in turn may lead to blindness.\footnote{Based upon
techniques of Goldblaum, et al., in Comparing Machine Learning Classifiers
for Diagnosing Glaucoma,
\emph{Investigative Opthamology}, \textbf{43}, 2002, 162-169. This project assumes some background
knowledge in probability.}

\newpar
An early warning sign of glaucoma is the loss of peripheral vision. As a means of assessing one's
peripheral vision, ophthalmologists conduct a \emph{visual field analysis}, also known as an \emph{automated perimetry exam}.
In this test, the patient sits in front of a concave dome, closes one of his or her eyes, and focuses the open eye
on a target in the center of the dome. In a random time pattern, the examiner flashes lights of varying intensities on 52 different locations within the dome. When the patient detects the appearance of a light, he or she clicks a button. A computer then plots a graph containing the 52 locations, along with the corresponding thresholds of detected light intensity.

\newpar
Figure \ref{Fi:visual_field} illustrates the result of a typical visual field examination performed on one eye. At each of the 52 field locations, the recorded value measures the extent to which the detected threshold light intensity differs from that of the dome background. Units are in decibels, with the smaller values corresponding to locations where a greater contrast between light intensity and background intensity is required before detection of light occurs. Consequently,
Figure \ref{Fi:visual_field} example suggests a deteriation of peripheral vision in the upper left quadrant.

\begin{figure}[h]
\centering
\scalebox{1.2}{\includegraphics{c:/images/OR/nonlinear/project_glaucoma_diagnosis_sample_vision_field.eps}  }
\caption{Sample Visual Field Readings}\label{Fi:visual_field}
\end{figure}

\newpar
Data for each eye, when combined with the individual's age, can be viewed as a vector in $\R^{53}$. The rate at which
peripheral vision is lost can then be estimated by comparing exam results obtained over a period of time. However, only a direct
examination of the optic nerve, conducted by an ophthalmologist (hereafter referred to as the ``expert''), can lead to
a definitive glaucoma diagnosis.

\newpar For the sake of simplicity and in order to visualize the underlying ideas better, we use data vectors
in $\R^2$ instead of $\R^{52}$. We assume that data has been collected from 16 individuals, with each individual contributing two visual field readings that constitute an ordered pair. Eight of the individuals have Glaucoma.

\newpar
In the context of glaucoma diagnosis, a \emph{machine learning classifier} can be viewed as a tool for predicting
the presence of glaucoma based solely upon the results of a visual field exam. More specifically, given the results of
a visual field alone, this tool predicts, with a certain degree of probability, that an expert would diagnose a person
as having glaucoma.

\newpar
To begin the construction of a simple type of classifier, we define the \emph{multivariate Gaussian} probability distribution
in $\R^n$ having mean, $\muvec \in \R^n$, and standard deviation, $\sigma$.\footnote{The function in (\ref{E:glaucoma_gaussian}) is a special form
of the multivariate Gaussian distribution. It assumes that the components of $\vecx$ have a common standard deviation. A
more general form of this distribution does not make this assumption and instead utilizes a \emph{covariance matrix}.}
The probability density function for this
random variable is given by
\begin{equation}\label{E:glaucoma_gaussian}
f_{\muvec,\sigma}(\vecx)= \frac{1}{(2\pi\sigma^2)^{n/2}}e^{\frac{-\|\vecx -\muvec \|^2}{2\sigma^2}},
\end{equation}
which satisfies $\ds \int_{\R^n} f_{\muvec,\sigma}(\vecx)\, dA=1$.

\newpar
For a set of data points in $\R^n$, we can determine the mean, $\muvec$, and standard deviation, $\sigma$, using
a procedure known as \emph{Maximum Likelihood Estimation} (MLE). This tool determines parameter values that
most closely correspond to the given data. Assume first that this data is represented by $\ds \tilde{\vecx}_1,\tilde{\vecx}_2,\ldots,\tilde{\vecx}_N$, and
that each of these points is drawn from a common distribution having density $f_{\muvec,\sigma}(\vecx)$.
If the points are sampled in an independent
manner, then the joint probability distribution corresponding to the sampling process has a density function
given by the product,  $\ds \prod_{i=1}^{N}f_{\muvec,\sigma}(\vecx_i)$.  The \emph{Maximum Likelihood Function},
$L$, is obtained by substituting the data into this product and viewing the resulting expression as a function
of the original parameters, $\muvec$ and $\sigma$. In other words,
\begin{equation}\label{E:glaucoma_likelihood}
L(\muvec,\sigma)=  \prod_{i=1}^{N}f_{\muvec,\sigma}(\tilde{\vecx}_i).
\end{equation}
Note that $L$ is a function of $n+1$ variables consisting of the components of $\muvec$ together with $\sigma$.
MLE selects $\muvec$ and $\sigma$ that maximize the function $L$.

\be
\item[1.] The \emph{Negative Log-Likelihood Function}, $NLL$, corresponding to $L$ is given by $\ds -\ln L$. Note that maximizing
$L$ corresponds to minimizing the Negative Log-Likelihood Function. Calculate a formula for the $NLL$ when the distribution is Gaussian as in (\ref{E:glaucoma_gaussian}). Use logarithm properties to represent this expression as a sum containing no exponential terms.

\item[2.] Visual field readings for the eight nonglaucomic individuals are summarized in Table \ref{Ta:nonglaucomic_data}

\begin{table}[h]
\centering
\begin{tabular}{|c||c||c|}
\hline \hline
Individual & Field Location 1 & Field Location 2 \\ \hline
  1 & 6.6 & 2.5 \\  \hline
2 &-2.1& 4.7  \\  \hline
3 &-2.1& -4.3 \\  \hline
4 &-.5& -5.0\\  \hline
5 &1.3& 5.8\\  \hline
6 &1.4& -1 \\  \hline
7 &-.4& -.3\\  \hline
8 &-1.6& -1\\  \hline
\hline
\end{tabular}\caption{Ordered Pairs of Visual Field Readings for Eight Nonglaucomic Individuals.}
\label{Ta:nonglaucomic_data}
\end{table}

\noi
Calculate the Negative Log-Likelihood Function that corresponds
to this data, under the assumption that it is drawn from a Gaussian distribution in $\R^2$.
\item[3.] Use appropriate unconstrained optimization techniques to determine the components of $\muvec$ and the value
of $\sigma$ that
minimize the $NLL$. You may wish to use numeric techniques or to take advantage of the Maple \mf{GlobalSolve} command within
the \mf{GlobalOptimization} toolbox. Depending upon your computing platform, you may wish to utilize the \mf{timelimit} command, which takes the form \mf{timelimit(t=\emph{maximum allowed time in seconds},\emph{calculation})}. Here \mf{\emph{calculation}} can take a wide variety of forms, including the \mf{GlobalSolve} command or a call to a numeric optimization procedure.
\item[4.] Use the \mf{plot3d} command to plot the resulting density function, $f_{\muvec,\sigma}(\vecx)$. Your result should
depict a bell-shaped surface centered about $\muvec$.
\ee

\noi
Visual field readings for the eight glaucomic individuals is given in Table \ref{Ta:glaucoma_data}. These data points, along with those of the nonglaucomic individuals, are depicted
in Figure \ref{Fi:glaucoma_data_plot}.


\begin{table}[h]
\centering
\begin{tabular}{|c||c||c|}
\hline \hline
Individual & Field Location 1 & Field Location 2 \\ \hline
1& -7.7 & -12.9\\  \hline
2& -4& -18\\  \hline
3& -11.7& -14.6\\  \hline
4& -6.1& -15.7\\  \hline
5& -17& -4.8\\  \hline
6& -22& -8.4\\  \hline
7& -14.2& -8.0\\  \hline
8& -14.3& -9.5\\  \hline
\hline
\end{tabular}\caption{Ordered Pairs of Visual Field Readings for Eight Glaucomic Individuals.}
\label{Ta:glaucoma_data}
\end{table}



\begin{figure}[h]
\centering
\psfrag{field1}[Bc][Bc]{$\mbox{Location 1 (dB)}$}
\psfrag{field2}[Bc][Bc]{$\mbox{Location 1 (dB)}$}
\scalebox{.6}{\includegraphics{c:/images/OR/nonlinear/glaucoma_data_point_plot.eps} }
\caption{Plot of Field Data for 16 Individuals}
\label{Fi:glaucoma_data_plot}
\end{figure}

\newpar
Figure \ref{Fi:glaucoma_data_plot} suggests that the data for glaucomic individuals, unlike that of the nonglaucomic
individuals, forms two \emph{clusters}. Consequently, the use of a single Gaussian distribution is not appropriate.
Instead, we utilize a \emph{Gaussian Mixture Model} (GMM), a distribution formed by constructing a weighted sum, or linear combination, of two Gaussian distributions. The density function of our GMM is given by
\begin{equation}\label{E:glaucoma_gaussian_mixture}
f_{\muvec_1,\sigma_1,\muvec_2,\sigma_2,w_1,w_2}(\vecx)=
w_1\frac{1}{(2\pi\sigma_1^2)^{n/2}}e^{\frac{-\|\vecx -\muvec_1 \|^2}{2\sigma_1^2}}+ w_2 \frac{1}{(2\pi\sigma_2^2)^{n/2}}e^{\frac{-\|\vecx -\muvec_2 \|^2}{2\sigma_2^2}}.
\end{equation}

\noi
The conditions that $\ds f_{\muvec_1,\sigma_1,\muvec_2,\sigma_2}$ is nonnegative and satisfies
\begin{equation*}
\int_{\R^n} f_{\muvec_1,\sigma_1,\muvec_2,\sigma_2}(\vecx)\, dA=1,
\end{equation*}
implies that $w_1, w_2 \geq 0 \mbox{ and } w_1+w_2=1$.
We refer to the quantities $w_1$ as $w_2$ as \emph{mixing coefficients}. Because the plot of a Gaussian distribution in
$\R^2$ is bell-shaped, we might expect the plot of a GMM density function in $\R^2$ to resemble that given in
Figure \ref{Fi:GMM}.

\begin{figure}[h]
\centering
\psfrag{x1axis}[Bc][Bc]{$x_1$}
\psfrag{x2axis}[Bc][Bc]{$x_2$}
\scalebox{.75}{\includegraphics{c:/images/OR/nonlinear/project_glaucoma_gaussian_mixture_model.eps}}\caption{Sample Gaussian Mixture Model Density Function}\label{Fi:GMM}
\end{figure}

\newpar
To determine the density function $\ds f_{\muvec_1,\sigma_1,\muvec_2,\sigma_2,w_1,w_2}$ that best corresponds
to the data in Table \ref{Ta:glaucoma_data}, we again utilize Maximum Likelihood Estimation. In this situation, the
number of variables is considerably larger than for the simple Gaussian distribution. Namely, the estimation seeks
to determine the optimal values of the following quantities:
\bi
\item the components of $\muvec_1$, and $\muvec_2$;
\item the values of $\sigma_1$ and $\sigma_2$;
\item the mixing coefficients, $w_1$ and $w_2$.
\ei

\noi
But the greatest difference between the GMM as opposed to the simple Gaussian distribution, is that determining
the GMM becomes a {\underbar{constrained}} optimization problem. Specifically, the Maximum Likelihood Function,
\begin{equation}\label{E:glaucoma_GMM_likelihood}
L\left( \muvec_1,\sigma_1,\muvec_2,\sigma_2,w_1,w_2 \right)=  \prod_{i=1}^{n}f_{\muvec_1,\sigma_1,\muvec_2,\sigma_2,w_1,w_2}(\tilde{\vecx}_i).
\end{equation}
is maximized subject to the constraints $w_1, w_2 \geq 0 \mbox{ and } w_1+w_2=1$.

\be
\item[5.] Use Maximum Likelihood Estimation to determine the values of the
parameters, $\muvec_1$, $\sigma_1$, $\muvec_2$, $\sigma_2$, $w_1$, and $w_2$ that minimize
the Negative Log-Likelihood function, $NLL$ corresponding to
(\ref{E:glaucoma_GMM_likelihood}), subject to the mixing coefficient
conditions, $w_1, w_2 \geq 0 \mbox{ and } w_1+w_2=1$. Again, you may wish to
utilize the Maple \mf{GlobalSolve} and \mf{timelimit} commands in this process.
\item[6.] Using your resulting parameter values, plot the corresponding Gaussian Mixture Model corresponding to
the data from the glaucomic individuals.
\ee

\newpar
The process of determining the parameters $\muvec_1$, $\sigma_1$, $\muvec_2$, $\sigma_2$, $w_1$, and $w_2$ is referred to
as \emph{unsupervised learning}. The information from Table \ref{Ta:glaucoma_data} used in this process is
called the \emph{training data}.

\newpar
The resulting Gaussian Mixture Model, together with the results from problem (3) can then be used to assist in predicting whether an individual, other than the 16 tested thus far, has glaucoma
based solely upon his or her visual field results and not from an examination conducted by the expert. Essential to this prediction is Bayes' Theorem from probability.

\begin{theorem}\label{Th:BayesTheorem}
Suppose $A$ and $B$ are events from some sample space and that the probability of these events are denoted by
$P(A)$ and $P(B)$ respectively. Then the conditional probability, $P(A\mid B)$ can be written as
\begin{equation}\label{E:BayesStatement}
P(A\mid B)= \frac{P(B\mid A)P(A)}{P(B)}.
\end{equation}
\end{theorem}

\noi
We call the probability, $P(A)$, the \emph{prior probability}, as it takes no prior information about event $B$ into account. The conditional probability, $P(A\mid B)$ is the \emph{posterior probability}; it does factor in information regarding $B$.

\be
\item[7.]Suppose an individual, different from the 16 tested thus far, has a visual field exam resulting in
the data vector $\tilde{\vecx}$. Let $G$ denote the event that an individual has glaucoma; $G^c$ then is the event
he or she is not glaucomic. Justify each of the following steps:
\begin{align*}
P(G \mid \tilde{\vecx}) &= \frac{P(\tilde{\vecx} \mid G)P(G)}{P(\tilde{\vecx})} \\
&= \frac{P(\tilde{\vecx} \mid G)P(G)}{P(\tilde{\vecx} \mid G)+P(\tilde{\vecx} \mid G^c)}\\
\end{align*}
\item[8.] Assume that $50$\% of all individuals tested have glaucoma in general. Use this fact, along with the
Gaussian distribution for nonglaucomic individuals and the GMM for those with glaucoma, to determine the probability
that each of the individuals having field results in Table \ref{Ta:test_data} does in fact have glaucoma.

\begin{table}[h]
\centering
\begin{tabular}{|c||c||c|}
\hline \hline
Individual & Field Location 1 & Field Location 2 \\ \hline
1& -2 & -2\\  \hline
2& 0& 1\\  \hline
3& 5& 8\\  \hline
\hline
\end{tabular}\caption{Visual Field Results for Three New Individuals.}
\label{Ta:test_data}
\end{table}
\ee

\endinput
