\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}
\usepackage{Sweave}
\begin{document}
\maketitle
\thispagestyle{empty}
\begin{center}
  Due date: Sept 21 2009, at 5pm
\end{center}


Please keep code, apart perhaps from short code fragments, in an
appendix that is separate from the main body of the report. You will
be given marks for presentation and layout.
\vspace*{9pt}

This exercise will work with experimental control and treatment groups
in the data set \texttt{nswdemo}, and with non-experimental comparison
groups in the data sets \texttt{cps1} and \texttt{psdi1}.  All these
datasets are included in the \textit{DAAG} package.


\section*{Useful Functions}
The following functions will be useful.

\paragraph{Extract a random sample of rows:}
The following function is designed to extract a random sample of
\texttt{m} rows from a data frame (or matrix).
\begin{Schunk}
\begin{Sinput}
 samprows <- function(df, m)df[sample(1:nrow(df), m), ]
\end{Sinput}
\end{Schunk}

\paragraph{Function to evaluate predictive accuracy, given model predictions:}

\begin{Schunk}
\begin{Sinput}
 confusion <- function(actual, predicted, names=NULL, 
                        printit=TRUE, prior=NULL){
    if(is.null(names))names <- levels(actual)
    tab <- table(actual, predicted)
    acctab <- t(apply(tab, 1, function(x)x/sum(x)))
    dimnames(acctab) <- list(Actual=names,
                        "Predicted (cv)"=names) 
    if(is.null(prior)){
      relnum <- table(actual)
      prior <- relnum/sum(relnum)
      acc <- sum(tab[row(tab)==col(tab)])/sum(tab)
    } else
    {
      acc <- sum(prior*diag(acctab))
      names(prior) <- names
  }    
    if(printit)print(round(c("Overall accuracy"=acc, 
                             "Prior frequency"=prior),4))
    if(printit){
      cat("\nConfusion matrix", "\n")
      print(round(acctab,4))
    }
    invisible(list(accuracy=acc, confusion=acctab, prior=prior))
  }
\end{Sinput}
\end{Schunk}

\paragraph{Function to fit \texttt{lda()} or \texttt{qda()} model,
  and evaluate accuracy:}

This uses leave-one-out cross-validation, which is not ideal. However
it is a good way to get a quick rough initial assessment.


\begin{Schunk}
\begin{Sinput}
 check.model <- function(form, df=nswdem, dp=4, prior=NULL, discfun=lda){
      categ <- all.names(form)[2]
      if(is.null(prior)){
          df.disc <- discfun(form, data=df, CV=TRUE) 
          acc <- confusion(df[, categ], df.disc$class, printit=FALSE)$accuracy
      } else {
         df.disc <- discfun(form, data=df, CV=TRUE, prior=prior)
         acc <- confusion(df[, categ], df.disc$class, printit=FALSE, 
                          prior=prior)$accuracy   
     }
      print(paste("Model is", deparse(form)))
      print(c(Accuracy = round(acc, dp)))
      invisible(acc)
  }
\end{Sinput}
\end{Schunk}

\section*{Datasets}

We will want the following datasets

\paragraph{\texttt{nswdem} (Experimental data):}  Add the variable \texttt{logre75}
to the dataframe \texttt{nswdemo}
\begin{Schunk}
\begin{Sinput}
 library(DAAG); library(DAAGxtras)
 library(splines)
 dim(nswdemo)
\end{Sinput}
\begin{Soutput}
[1] 722  10
\end{Soutput}
\begin{Sinput}
 library(MASS)
 unique(sort(nswdemo$re75))[2]/2
\end{Sinput}
\begin{Soutput}
[1] 37.17173
\end{Soutput}
\begin{Sinput}
 nswdem <- nswdemo
 nswdem$logre75 <- log(nswdemo$re75+37)
\end{Sinput}
\end{Schunk}

\paragraph{\texttt{nswplus}: Experimental data, plus non-experimental controls} Create a dataset \texttt{nswplus} that combines the
experimental data with the two sets of non-experimental controls.
\begin{Schunk}
\begin{Sinput}
 library(DAAG); library(DAAGxtras)
 nswplus <- rbind(nswdemo, cps1, psid1)
 nswplus$gp <- factor(rep(c("nsw","cps","psid"), 
                          c(722, nrow(cps1), nrow(psid1))),          
                      levels=c("nsw","cps","psid"))
 nswplus$logre75 <- log(nswplus$re75+1)
 nswplus$gp1 <- 1-as.numeric(nswplus$gp=="nsw")
   # gp1 distinguishes between the experimental and other data
 ## Now check that the new dataset has been correctly created
 all.equal(nswplus[nswplus$gp=="nsw",1:10],nswdemo)
\end{Sinput}
\begin{Soutput}
[1] TRUE
\end{Soutput}
\begin{Sinput}
 all.equal(nswplus[nswplus$gp=="cps",1:10],cps1, check.attributes=FALSE)
\end{Sinput}
\begin{Soutput}
[1] TRUE
\end{Soutput}
\begin{Sinput}
 all.equal(nswplus[nswplus$gp=="psid",1:10],psid1, check.attributes=FALSE)
\end{Sinput}
\begin{Soutput}
[1] TRUE
\end{Soutput}
\end{Schunk}

\section*{Questions}

\subsection*{Function help page}


\paragraph{Question 1:} 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{Schunk}
\begin{Sinput}
 prompt(confusion)
\end{Sinput}
\end{Schunk}
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]}

\subsection*{Data transformation:}

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


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


Here is code that may be used to compare the control and treatment
groups in the experimental \texttt{nswdem} dataset. 
\begin{Schunk}
\begin{Sinput}
 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{Sinput}
\end{Schunk}

\paragraph{Question 3a:} With \texttt{form} taken to be
whichever of \texttt{form1}, \texttt{form2} and \texttt{form3}
gives the best accuracy, run the following code:
\begin{Schunk}
\begin{Sinput}
 nswdem1.lda <- lda(form, data=nswdem)
 score <- predict(nswdem.lda)$x
 plot(densityplot(~score, groups=nswdem$trt, auto.key=list(columns=2)))
\end{Sinput}
\end{Schunk}
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]}

\paragraph{3b:} 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]}

\subsection*{Pre-treatment comparison between
\texttt{nswdem},  \texttt{cps1} and  \texttt{psdi1}}
Check the number of rows in each of \texttt{nswdem}, \texttt{cps1}
and \texttt{psid1}.  

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{Schunk}
\begin{Sinput}
 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{Sinput}
\end{Schunk}

\paragraph{Question 4:} 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{Schunk}
\begin{Sinput}
 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{Sinput}
\end{Schunk}
\hfill \textbf{\em [2 marks]}

\subsection*{Comparison between the three groups, using \texttt{lda()}}
\paragraph{Question 5:} 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, implemeted in
\texttt{gam()}, for automatic choice of the degree of spline curve.)
In what follows, we try normal splines of degree 4.
\begin{Schunk}
\begin{Sinput}
 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{Sinput}
\end{Schunk}
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{Schunk}
\begin{Sinput}
 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{Sinput}
\end{Schunk}
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]}

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

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{Schunk}
\begin{Sinput}
 ## 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{Sinput}
\end{Schunk}
\paragraph{Question 6a:}  Annotate the above code, explaining what each
step does, and the choice of function arguments.\hfill \textbf{\em [3 marks]}

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

\section*{Overview}

\paragraph{Question 7:}  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]}

\subsection*{TOTAL MARKS = 20}

\end{document}
