\documentclass{article}
\usepackage{url}
\usepackage{natbib}
\usepackage{float}
% \usepackage{exercises}
\usepackage{nextpage}
\setcounter{tocdepth}{1}
\setlength{\voffset}{-2.5cm}
\setlength{\hoffset}{-0.5cm} \setlength{\textwidth}{14cm}
\setlength{\textheight}{23cm}
\title{Assignment 2, Math 3346, 2009}
\author{Lecturer: John Maindonald}
\begin{document}
\maketitle
\thispagestyle{empty}
\begin{center}
  Due date: Sept 21 2009, at 5pm
\end{center}

\section*{General Comments}
There were two issues that arise from the use of made of these data
in the econometrics literature:
\begin{enumerate}
\item {\em Is there any reason to suspect the randomisation in the NSW
    experimental study?}\\
  The question was designed to elicit consideration about what might
  be expected if the randomisation was properly done, and to consider
  clues that might suggest that it was not properly done. If there is
no reason to suspect the randomisation, then it ti reasonable to pool
the control and treatment groups, for purpose of the comparison with
the \texttt{cps1} and \texttt{psid1} control groups.  (For these purposes,
\texttt{re78} is of course ignored.)
\item {\em The \texttt{cps1} and \texttt{psid1} data, which was from
    observational studies individuals in neighbouring areas who had
    experienced similar employment difficulties, have been suggested
    for possible use as alternative controls.  Might that have yielded
    reasonable results?}\\
  No, not using the data as they stand.  The task set by the
  assignment had as a strong theme the assembling of the evidence, in
  a convincing manner.  The questions were, essentially, hand-holding.
\end{enumerate}


\section*{Questions}

\subsection*{Function help page}


\paragraph{Question 1:} {\em Test the function \texttt{confusion()} by 
using the function \texttt{lda()} to do discriminant calculations,
and then using this fucntion to determine the predictive accuracy.

By typing 
\begin{verbatim}
 prompt(confusion)
\end{verbatim}
you can get a skeleton for a help page for the function \texttt{confusion()}.
Fill in the details for the help page.}
\hfill \textbf{\em [3 marks]}

NB: help pages should be genuinely helpful!

\subsection*{Data transformation:}

\paragraph{Question 2:}
{\em Why is it preferable to work with \texttt{logre75}, rather than 
\texttt{re75}?\hfill }
\textbf{\em [1 mark]}

In the untransformed data, a small number of individuals with large
values of \texttt{re75} get large weight (or leverage).

\subsection*{Treated vs controls in  \texttt{nswdem},
  based on pre-treatment variables}

{\em
Here is code that may be used to compare the control and treatment
groups in the experimental \texttt{nswdem} dataset:}
\begin{verbatim}
 library(lattice)
 form1 <- trt ~ age+educ+black+hisp+marr+nodeg+logre75
 form2 <- trt ~ ns(age,3)+educ+black+hisp+marr+nodeg+logre75
   # form2 allows for the possibility that the effect of age mayh be nonlinear
 form3 <- trt ~ (ns(age,3)+logre75)*(educ+black+hisp+marr)+nodeg
   # form3 allows for interaction effects involving continuous variables
 ## Try also the equivalent models with form2 and form3.
 check.model(form1)
\end{verbatim}

\paragraph{Question 3a:} {\em With \texttt{form} taken to be whichever
  of \texttt{form1}, \texttt{form2} and \texttt{form3} gives the best
  accuracy, run the following code:
\begin{verbatim}
 nswdem1.lda <- lda(form, data=nswdem)
 score <- predict(nswdem.lda)$x
 plot(densityplot(~score, groups=nswdem$trt, auto.key=list(columns=2)))
\end{verbatim}
Given that subjects had been randomly
divided between treated and control (though with a greater number in
control), is the graph much what might be expected? Explain.}\newline
\hspace*{\fill} \textbf{\em [2 marks]}
Yes, it is. There are several ways that the answer could be justified:
\begin{itemize}
\item Bootstrap resampling;
\item Simulation;
\item Randomly permute labels between treatment and control groups.
\end{itemize}
No-one provided such a justification!

\paragraph{3b:} {\em Compare the accuracies with the accuracy achieved by
assigning all observations to the most frequent category.
(In output from \texttt{rpart()}, this is the {\em root node accuracy}.)
Does the best accuracy that was achieved give any
reason to suspect the randomization? Explain.}\newline
\hspace*{\fill} \textbf{\em [2 marks]}

Thoughts that this was intended to stimulate include:
\begin{itemize}
\item If the randomisation was effective, what should an effective
dioscriminant analysis show?  
\item why the comparison the accuracy from asssignment to the most
frequent group?
\end{itemize}

\subsection*{Pre-treatment comparison between
\texttt{nswdem},  \texttt{cps1} and  \texttt{psdi1}}

{\em The following will use the smoothing spline methodology of the
function \texttt{gam()} in the \textit{mgcv} package to check
possible transformation of the variables  \texttt{age}, \texttt{educ}, 
and \texttt{logre75}.  The software automatically makes what should be
a sensible choice of the amount of smoothing, corresponding to each
of the terms \texttt{s(age)}, \texttt{s(educ)}, and \texttt{s(logre75)}.}
\begin{verbatim}
 library(mgcv)
 form <- gp1 ~ s(age)+s(educ)+s(logre75)+black+hisp+marr+nodeg
 nswcps.gam <- gam(form, data=subset(nswplus, gp!="psid"))
 nswpsid.gam <- gam(form, data=subset(nswplus, gp!="cps"))
\end{verbatim}

\paragraph{Question 4:} {\em Would the same transformations that are effective
for the comparison between the experimental data and the \texttt{cps1}
data be useful for the comparison between the experimental data and the 
\texttt{psid1} data.  Justify your answer by reference to the graphs given by
\texttt{plot(nswcps.gam)} and \texttt{plot(nswpsid.gam)}. You might
want to do:}
\begin{verbatim}
 opar <- par(mfrow=c(2,3), mar=c(3.6,3.6,0.6,0.6), mgp=c(2.25,.5,0))
 plot(nswcps.gam)
 plot(nswpsid.gam)
 par(opar)
\end{verbatim}
\hfill \textbf{\em [2 marks]}

The smoothing curves are broadly similar in their shape.  So, at least
as a starting point, use of the same transformation would probably be
satisfactory.

\subsection*{Comparison between the three groups, using \texttt{lda()}}
\paragraph{Question 5:} {\em Use \texttt{lda()} to obtain a graphical
comparison between the three groups. It may be desirable to account
for nonlinear effects from any or all of \texttt{age}, \texttt{educ}
and \texttt{logre75}.  The graphs from Question 4 suggest that spline
curves would be worth trying when we use \texttt{lda()}.  (The
function texttt{lda()} does not have the functionality, implemented in
\texttt{gam()}, for automatic choice of the degree of spline curve.)
In what follows, we try normal splines of degree 4.
\begin{verbatim}
 library(splines)
 fm1 <- gp ~ age+educ+black+hisp+marr+nodeg+logre75
 fm2 <- gp ~ ns(age,4)+ns(educ,4)+ns(logre75,4)+black+hisp+marr+nodeg
\end{verbatim}
Using whichever of \texttt{fm1} and \texttt{fm2} gives
the best cross-validation predictive accuracy, fit the model,
and graph the results.  For example:
\begin{verbatim}
 nswplus.lda <- lda(fm1, data=nswplus, prior=rep(1,3)/3)
 scores <- predict(nswplus.lda)$x
 library(lattice)
 xyplot(scores[,1] ~ scores[,2], groups=nswplus$gp,
         auto.key=list(columns=3), par.settings=simpleTheme(cex=0.5, alpha=0.2))
\end{verbatim}
The graph has a very striking feature.  Can you shed any light on it?
Hint: Consider looking at ranges of variables, within suitable splits of
the data.}\hfill \textbf{\em [2 marks]}

The results (better accuracy, though not by a large margin) suggest
use of \texttt{fm2}.  Both graphs shows a separation into two groups,
weaker in the \texttt{fm2} graph than in the \texttt{fm1} graph.  The
seaparation is between blacks and non-blacks.  In part, this is
however a difference in ages.  Better modeling of the effects of age,
\texttt{re75} and \texttt{educ} (especially) reduces the apparent
effect of skin colour.

\subsection*{Comparison between the three groups, using \texttt{randomForest()}}

{\em Proximities from \texttt{randomForest} calculations can be used
to get a plot that reflects the way that the data have been
classified.  Points are close together or widely separated according
as they have, or have not, followed the same path down the trees to
the terminal node.  The proximity is, for each pair of points, the
proportion of trees in which they have come together at the same
terminal node.  With more than several thousand points, the number
gets so large that it creates problems for the computations.}
\begin{verbatim}
 ## Extract subset of data for plotting
 m <- cumsum(c(nrow(nswdem), nrow(cps1), nrow(psid1)))
 n1 <- sample(1:m[1], 500)
 n2 <- sample((m[1]+1):m[2], 500)
 n3 <- sample((m[2]+1):m[3], 500)
 take <- c(n1,n2,n3)
 nswplus.rf <- randomForest(fm1, data=nswplus, sampsize=rep(722,3))
 nswplus.pred <- predict(nswplus.rf, newdata=nswplus[take, ], proximity=TRUE)
 sim <- 1-nswplus.pred$proximity
 scores.rf <- cmdscale(sim)
 xyplot(scores.rf[,1] ~ scores.rf[,2], groups=nswplus$gp[take],
         auto.key=list(columns=3), par.settings=simpleTheme(cex=0.5, alpha=0.5))             
\end{verbatim}
\paragraph{Question 6a:}  {\em Annotate the above code, explaining what each
step does, and the choice of function arguments.}\hfill \textbf{\em [3 marks]}

\paragraph{Question 6b:} {\em Comment on the comparison between the 
 \texttt{randomForest()} analysis and the \texttt{lda()} analysis.}\hfill 
\textbf{\em [2 marks]}

In both plots, the three groups occupy different regions of the graph.
The plot based on \texttt{lda()} is a low-dimensional projection,
after rotation, of the space spanned by the columns of the model
matrix.  The plot from \texttt{randomForest()} reflects the way that
the algorithm has classified the observations.  The groupings do not
reflect the split into black and non-black. That split is, in the
\texttt{lda()} analysis, a proxy for other differences.

The overall accuracy for the random forest analysis was 79.9\%,
against an overall accuracy (leave one out CV) of 68.3\% for the
analysis that used \texttt{lda()}, with \texttt{fm2}.  (If 
\texttt{fm1} is used, the accuracy drops to 62.1\%.)
Note that, for fair comparison with the \texttt{randomForest()}
result, the prior must be specified as \texttt{c(1/3, 1/3, 1/3)},
both when calling \texttt{lda()} and when calling \texttt{confusion()}.
Clearly \texttt{randomForest()} does a much more effective job of
discrimination.

\section*{Overview}

\paragraph{Question 7:}  {\em Write brief notes on what can be deduced
from the several steps of the analysis.  NB: Keep in mind the reasons
for collecting and comparing these various sets of data.}\hfill
\textbf{\em [4 marks]}

See the ``General Comments'' at the beginning.  There are huge
differences between the three groups.  As they stand, the \texttt{cps1}
and \texttt{psid1} data are not suitable for use as controls.\\[6pt]
[Some possibilities that have been investigated are:
\begin{itemize}
\item Identify subgroups that are similar.\newline
[Indications are that such subgroups would be very small; look at the plots]
\item Propensity score methods try to identify a single (constructed)
covariate that can account for differences in explanatory variables as
a whole.  Under quite strong assumptions, which are open to limited
checking, the propensity score can be used to account for pre-treatment
group differences.
\end{itemize}
]

\subsection*{TOTAL MARKS = 21}

\end{document}
