\documentclass{article}
\usepackage{url}
\usepackage{natbib}
\usepackage{float}
% \usepackage{exercises}
\usepackage{nextpage}
\setcounter{tocdepth}{1}
\setlength{\voffset}{-3.5cm}
\setlength{\hoffset}{-2cm} \setlength{\textwidth}{15.5cm}
\setlength{\textheight}{24.75cm}
\title{Math3346 Assignment 1 -- Solutions and Commentary}
\author{Lecturer: John Maindonald}

\usepackage{Sweave}
\begin{document}
\maketitle
\thispagestyle{empty}








\section*{General Comments}
\begin{itemize}
\item Use tables to summarize information.  Use text to draw attention
  to major features in the tables. Place the main focus on features that 
  are clear and unequivocal, with more minor matters consigned to the
  later brief comment or to footnotes or endnotes.
\item In an assignment such as this, most code should be included,
  properly referenced, at the end of the document. Including limited
  amounts of code in the text is acceptable, providing it has a strong
  ``literate programming'' character!\newline [In a report for a
  manager, omit code from the main document.]
 \item Round proportions to (in these data) 3 or 4 decimal places.
 \item Be judicious in the choice of graphs. If you have 10 or 12
   graphs spread over 4 or 5 pages, you likely need to find a more
   succinct form of presentation.   Aim at a layout where it
   will be easy to identify the major differences.  Overlaid plots of
   density estimates sometimes work well, for an audience that has some
   modest level of technical sophistication.
 \item Some forms of presentation that can be effective with an
   audience that has the right level of technically sophistication are
   not suited to use in reports that are intended for a general
   audience. Thus in considering use of, for example, density plots
   and boxplots, keep the audience in mind.
\item Where a variable has missing values the total of zeros (or other
  such total) will be the total of those known to be zero.  On the
  (common) assumption that data are missing at random, division by the
  total of non-missing values gives an estimate of the proportion of
  zeros.  If data are not missing at random, then the number that
  results may be biased, perhaps so seriously biased that it is
  uninterpretable.  This seems an issue for the variable \texttt{re74}
  in the present data.
  \end{itemize}

\section*{Solutions}
\begin{enumerate}
\item Write brief notes on each of the
columns of data, noting whether columns should be treated as numeric
or categorical. \hfill [1 mark]\newline 
{\em NB: Describe, where relevant, what codes 0 and 1 mean.}
\newline
Document (briefly) each of 
  \texttt{countNA()}, \texttt{countzeros()} and \texttt{countlist()}
  \hfill [2 marks]
\item For each column of the data and for each of the four
  groups
  \begin{enumerate}
  \item Determine the number of missing values. \hfill [1 mark]
\begin{Schunk}
\begin{Soutput}
         age educ black hisp marr nodeg re74 re75 re78
cf-cps1    0    0     0    0    0     0    0    0    0
cf-psid1   0    0     0    0    0     0    0    0    0
exp-ctl    0    0     0    0    0     0  165    0    0
exp-trt    0    0     0    0    0     0  112    0    0
\end{Soutput}
\end{Schunk}
  \item For each of \texttt{re74}, \texttt{re75} and
    \texttt{re78}, what number and proportion are zero. \hfill [1 mark]
\begin{Schunk}
\begin{Soutput}
     cf-cps1 cf-psid1 exp-ctl exp-trt cf-cps1 cf-psid1 exp-ctl exp-trt
re74    1913      215     195     131   0.120    0.086   0.750   0.708
re75    1748      249     178     111   0.109    0.100   0.419   0.374
re78    2172      286     129      67   0.136    0.115   0.304   0.226
\end{Soutput}
\end{Schunk}
{\em Note: The proportion of (known) zeros for \texttt{re74} has been
  calculated relative to non-missing data.  This is the correct figure
  if missingness does not carry any information about income.  For the
  present experimental data, this implies a massive change ($>$ 30\%,
  from $\sim$ 0.4 to $>$ 0.7) in the proportion of non-earners,
  between 1974 and 1975.  As there is no comparable change in the
  \texttt{cps1} and \texttt{psid1} data, this seems unlikely.  Missing
  values for \texttt{re74} seem then unlikely to be missing at random.
  The variable \texttt{re74} has therefore been ignored in the subsequent 
  discussion.
  
  To calculate (known) zeros as a proportion of all data, change 
  \texttt{length(na.omit(x))} to \texttt{length(x)}. [Why would
  one want this?].
}
  \end{enumerate}

\item Now examine \texttt{re74}, comparing control and treatment data:
  \begin{enumerate}
  \item Compare the proportion of \texttt{NA}s between 
    the experimental control and treatment. \hfill [$\frac{1}{2}$ mark]

\begin{Schunk}
\begin{Soutput}
    0     1 
0.388 0.377 
\end{Soutput}
\end{Schunk}
  \item Compare the proportion of 0's in each of \texttt{re74}, \texttt{re75}
    and \texttt{re78} (obviously \texttt{NA}s have to 
    be excluded) between the four groups. \hfill [$\frac{1}{2}$ mark]
\begin{Schunk}
\begin{Soutput}
     cf-cps1 cf-psid1 exp-ctl exp-trt
re74   0.120    0.086   0.750   0.708
re75   0.109    0.100   0.419   0.374
re78   0.136    0.115   0.304   0.226
\end{Soutput}
\end{Schunk}
  \item For columns that are numeric with more than two unique values,
    determine for each of the four groups the range of values.
    \hfill [1 mark]
\begin{Schunk}
\begin{Soutput}
         age-l age-u educ-l educ-u re74-l re74-u re75-l re75-u re78-l re78-u
cf-cps1     16    55      0     18      0  25862      0  25244      0  25565
cf-psid1    18    55      0     17      0 137149      0 156653      0 121174
exp-ctl     17    55      3     14      0  39571      0  36941      0  39484
exp-trt     17    49      4     16      0  35040      0  37432      0  60308
\end{Soutput}
\end{Schunk}
\end{enumerate}    
\item Provide graphs that conveniently summarise differences between
  the four groups, with respect to \texttt{age} and \texttt{re75}.
\hfill  [4 marks]


\setkeys{Gin}{width=1.15\linewidth}
\begin{figure}[H]
\includegraphics{assignment1sol-do-densities}
\caption{Densityplot comparisons, for Age, for Years of education
  and for \texttt{log(re75+37)}.
  The logarithmic transformation leads to a more symmetric
  distribution.  The plot for \texttt{log(re75+37} has a spike at zero
  that, in the density plot, is spread either side of zero.
}
\end{figure}

\setkeys{Gin}{width=1.15\linewidth}
\begin{figure}[H]
\includegraphics{assignment1sol-do-box}
\caption{Boxplot comparisons, for Age, for Years of education, 
  and for \texttt{log(re75+37)}.
  The logarithmic transformation leads to a more symmetric
  distribution.  The plot for \texttt{log(re75+37} has a spike at zero.}
\end{figure}

Density plots do not work well with \texttt{educ}. Hence we try
histograms.
\setkeys{Gin}{width=1.15\linewidth}
\begin{figure}[H]
\includegraphics{assignment1sol-do-hist}
\caption{Histograms, used to compare distributions of Years of education.}    
\end{figure}

{\em
\begin{itemize}
\item With one density plot for all data area under the curve gives an
indication the density of points, but the reader has to know or warned
of the spike at zero, and to work on interpreting the graph.  
Alternatively, zeros can be omitted, with information
provided on the proportion of zeros.  Histograms are also acceptable
(with a relatively large number of data points, their discreteness is
not too much of an issue).  They do not however readily lend themselves
to overlaying.

\item A logarithmic scale seems preferable for \texttt{re75}; otherwise the 
graph focuses unduly on a small number of very large values.  

\item Half the minimum non-zero value may be a reasonable offset
  when taking logarithms. Alternatively,
choose the offset that leads to a near symmetric density plot.
\end{itemize}
}
{\em Note:
  The assignment did not ask for a graph for years of education, 
      but it (or something like it) is needed for the comparison in Exercise
\ref{ex:major}.}
\item Use tables to summarise differences in categorical variables
between the two groups. \hfill [2 marks]
\begin{Schunk}
\begin{Soutput}
         black  hisp  marr nodeg
cf-cps1  0.074 0.072 0.712 0.296
cf-psid1 0.251 0.033 0.866 0.305
exp-ctl  0.800 0.113 0.158 0.814
exp-trt  0.801 0.094 0.168 0.731
\end{Soutput}
\end{Schunk}
\item What are the major differences between the four groups,
as evident from examining individual columns?  Comment especially
on any differences in the pre-training variables, i.e., all except
\texttt{re78}.\label{ex:major} \hfill [3 marks]
\newline
{\em The proportions in all categorizations except \texttt{hisp} are 
  very different. In the experimental groups there are many more blacks, 
  many more who had not  graduated from high school, and many fewer married.
  Here is the comparison between the proportion who had fewer than 12
  years education:
\begin{Schunk}
\begin{Soutput}
  cf-cps1  cf-psid1   exp-ctl   exp-trt 
0.2958354 0.3052209 0.8141176 0.7306397 
\end{Soutput}
\end{Schunk}
  Note also differences in the \texttt{age} distribution
  and in the proportion earning no income in 1975
  (those in the experimental groups were much younger, and
  more likely to have had no income in 1975. NB: You were asked to
draw attention to major differences.}
\item Do any differences in the pre-training variables have implications
  for the way that you might analyse the data, or the reliance that
  you might put on the results? \hfill [2 marks]
\newline
  {\em The experimental control and treatment groups should, because
    of the randomization, be similar on all pre-treatment variables.
    The various comparisons above suggest that they are acceptably
    similar.  There are however big differences between the
    experimental data and non-experimental controls. Unless there
    is some way to adjust for these differences, the non-experimental
    "controls" are clearly quite unsuitable for use as controls.}  
\item The aim of the study was to assess the effect of training.  Is
  it best to base the comparison between treatment and control group
  i) on \texttt{re78} alone, or ii) on \texttt{re78 - re75}?  Justify
  your answer.\hfill [2 marks]\newline
  {\em One might expect some correlation between income in 1978 and
    income in 1975. The plot
    \verb!plot(log(re78+37) ~ log(re75+37), data=nswdemo)!
    suggests that there is a small correlation.  Examination of the
    change in income might therefore give a fairer assessment.
    On the other hand, the contribution of a few large changes in income,
    from 1975 to 1978, may introduce undesirable noise. Some subtlety
    is required in deciding how best to proceed.
    }
\end{enumerate}
\vspace*{3pt}

\noindent [\textbf{TOTAL: 20 marks}]


\section*{Code}

\paragraph{Setup}
\begin{Schunk}
\begin{Sinput}
> opt <- options(SweaveHooks=list(fig=function(){par(cex.main=1.0, cex=0.8,
                           mar=c(4.1,3.6,1.6,1.1), pty="s", lwd=1,
                           mgp=c(2.25,0.5,0), tck=-0.025)}),
                           continue="  ")
> pdf.options(pointsize=8)
> library(DAAG)
\end{Sinput}
\end{Schunk}
\paragraph{Creation of a Suitable Data Objects}
\begin{Schunk}
\begin{Sinput}
> ## Note the coding for the experimental control & treatment groups
> table(nswdemo$trt)
> ## Code the non-experimental controls as 2 (cps1) and 3 (psid1)
> cps1$trt <- rep(2, dim(cps1)[1])
> psid1$trt <- rep(3, dim(psid1)[1])
> allsets <- rbind(nswdemo, cps1, psid1)
> allsets$trt <- factor(c("exp-ctl", "exp-trt", "cf-cps1", "cf-psid1")
                        [allsets$trt+1])
> ## Check that the numbers in the different groups make sense
> table(allsets$trt)
\end{Sinput}
\end{Schunk}

\paragraph{Counts of number of missing values, or zeros}
\begin{Schunk}
\begin{Sinput}
> ## Function that retrieves number of NAs for a single column
> countNA <- function(x)sum(is.na(x))
> ## Now use the function aggregate() to apply countNA() to
> ## allsets, indexed by trt
> aggregate(allsets, list(group = allsets$trt), FUN = countNA)
\end{Sinput}
\end{Schunk}

\begin{Schunk}
\begin{Sinput}
> countzeros <- function(x)sum(!is.na(x) & x==0)
> aggregate(allsets[, c("re74", "re75", "re78")], list(group=allsets$trt), 
            FUN=countzeros)
\end{Sinput}
\end{Schunk}

Alternative
\begin{Schunk}
\begin{Sinput}
> ## Function that retrieves number of NAs for each column of a data frame
> countlist <- function(z, statsfun=length)sapply(z, statsfun)
>   # Notice that statsfun has been given a default argument
> ## Apply this to the allsets data frame
> sapply(split(allsets, allsets$trt), FUN=countlist, statsfun=countNA)
\end{Sinput}
\end{Schunk}

Alternative data object
\begin{Schunk}
\begin{Sinput}
> listsets <- list("cf-cps1"=cps1[,-1], "cf-psid1"=psid1[,-1],
                   "exp-ctl"=subset(nswdemo, trt==0)[,-1],
                   "exp-trt"=subset(nswdemo, trt==1)[,-1])
> ## Check that the numbers in the different groups make sense
> sapply(listsets, nrow)
> ## Now do the calculations with listsets
> sapply(listsets, FUN=countlist, statsfun=countNA)
\end{Sinput}
\end{Schunk}

\subsection*{Code for exercises}

\begin{enumerate}
\item \ldots
\item
\begin{enumerate}
  \item
\begin{Schunk}
\begin{Sinput}
> t(sapply(listsets, FUN=countlist, statsfun=countNA))
\end{Sinput}
\end{Schunk}
\item 
\begin{Schunk}
\begin{Sinput}
> num <- sapply(split(allsets[, c("re74","re75","re78")], 
               allsets$trt), FUN=countlist, statsfun=countzeros)
> pr <- sapply(split(allsets[, c("re74","re75","re78")], 
               allsets$trt), FUN=countlist, 
         statsfun=function(x)countzeros(x)/length(na.omit(x)))
> cbind(num, round(pr,3))
\end{Sinput}
\end{Schunk}
\label{ex:numANDprop-zero}
\vspace*{3pt}

\fbox{\parbox[t]{\linewidth}{ Note: The above uses a ``put as much as
    possible onto one line'' style of code that is compact and easy to 
    manage, but which novices may find overly cryptic.
    The less condensed code that now follows is more explicit. It
    works with smaller units of data, which is a consideration if
    datasets are large relative to available memory. There is some
    redundancy, which should be removed if these calculations are
    repeated a large number of times. The more obvious ways to remove
    the redundancy lead to code that is, again, somewhat more cryptic.
      }}  
\begin{Schunk}
\begin{Sinput}
> ## Longer, less cryptic and clearer code, for exptl data only
> namre <- c("re74","re75","re78")
> expctlz <- sapply(subset(nswdemo[, namre], nswdemo$trt==0), countzeros)
> exptrtz <- sapply(subset(nswdemo[, namre], nswdemo$trt==1), countzeros)
> numz <- cbind(expctl=expctlz, exptrt=exptrtz)
> expctlp <- sapply(subset(nswdemo[, namre], nswdemo$trt==0), 
                    function(x)countzeros(x)/length(na.omit(x)))
> exptrtp <- sapply(subset(nswdemo[, namre], nswdemo$trt==1), 
                    function(x)countzeros(x)/length(na.omit(x)))
> propz <- cbind("expctl-p"=expctlp, "exptrt-p"=exptrtp)
> cbind(numz, propz)
\end{Sinput}
\begin{Soutput}
     expctl exptrt  expctl-p  exptrt-p
re74    195    131 0.7500000 0.7081081
re75    178    111 0.4188235 0.3737374
re78    129     67 0.3035294 0.2255892
\end{Soutput}
\end{Schunk}
  \end{enumerate}
\item
  \begin{enumerate}
    \item 
\begin{Schunk}
\begin{Sinput}
> round(with(nswdemo, sapply(split(re74,trt), 
                      function(x)countNA(x)/length(x))), 3)
\end{Sinput}
\end{Schunk}
\item
\begin{Schunk}
\begin{Sinput}
> round(sapply(split(allsets[, c("re74","re75","re78")], 
               allsets$trt), FUN=countlist, 
         statsfun=function(x)countzeros(x)/length(na.omit(x))), 3)
\end{Sinput}
\end{Schunk}
\item
\begin{Schunk}
\begin{Sinput}
> varnam <- c("age","educ","re74","re75","re78")
> ran <- sapply(split(allsets[, varnam], 
                allsets$trt), FUN=countlist, 
                statsfun=function(x)range(x,na.rm=TRUE))
> rownames(ran) <- paste(rep(varnam, rep(2,length(varnam))),"-",
                         rep(c("l","u"), length(varnam)), sep="")
> print(round(t(ran),0))
\end{Sinput}
\end{Schunk}
\end{enumerate}    
\item
\begin{Schunk}
\begin{Sinput}
> library(lattice)
> at1 <- pretty(allsets$age,5)
> labs1 <- paste(at1)
> at2 <- pretty(allsets$educ,5)
> labs2 <- paste(at1)
> at3 <- with(allsets, pretty(log(re75+37), 3))
> at3 <- c(log(37), at2[at2>log(75)])
> labs3 <- paste(round(exp(at2) - 37, 0))
> print(densityplot(~age+educ+log(re75+37), groups=trt, auto.key=list(columns=4),
               scales=list(x=list(at=list(at1,at2,at3), 
                           labels=list(labs1,labs2,labs3)), relation="free"), 
                    xlab=c("Age", "Education (yrs)", "Income in 1975"),
                    data=allsets))
\end{Sinput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
> print(bwplot(trt~age+educ+log(re75+37), outer=TRUE, 
               scales=list(x=list(at=list(at1,at2,at3), 
                           labels=list(labs1,labs2,labs3), relation="free")), 
               xlab=c("Age", "Education (yrs)", "Income in 1975"),
               data=allsets))
\end{Sinput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
> print(histogram(~educ | trt, data=allsets, layout=c(4,1)))
\end{Sinput}
\end{Schunk}
\item 
\begin{Schunk}
\begin{Sinput}
> cat2prop <- function(x) if(is.factor(x))sum(x==levels(x)[2])/length(x) else
                     sum(x==max(x))/length(x)       
> contnum <- match(c("trt",varnam), names(allsets))
> t(round(sapply(split(allsets[, -contnum], 
               allsets$trt), FUN=countlist, statsfun=cat2prop), 3))
\end{Sinput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
> ## Proportion with less than 12 years education
> with(allsets, sapply(split(educ,trt), function(x)sum(x<12)/length(x)))
\end{Sinput}
\end{Schunk}
\end{enumerate}  

\section*{Use of the {\em doBy} Package}

The \texttt{summaryBy()} function ({\em doBy} Package) has extensive
abilities for calculating summary information.  It returns a data
frame.  The following does the calculation required to calculate, for
each of the nominated splits of the data and for each of the nominated
columns, the number of zeros:
\begin{Schunk}
\begin{Sinput}
> library(doBy)
> ## The "0s" gets added to the names of the variables
> countzeros <- function(x)c("0s"=sum(!is.na(x) & x==0))
> num0 <- summaryBy(re74 + re75 + re78 ~ trt, data=allsets, FUN=countzeros)
\end{Sinput}
\end{Schunk}
The first column is a factor, with one row for each level of \texttt{trt}.
The level names are better used to label the rows. The first column
can then be omitted, leaving a data frame that has all columns numeric.
\begin{Schunk}
\begin{Sinput}
> row.names(num0) <- as.character(num0[,1])
\end{Sinput}
\end{Schunk}
Then \texttt{print(num0[, -1])} prints with the levels of
\texttt{trt} used as row names, while \texttt{print(t(num0[, -1]))}
prints without the quotes that would otherwise appear.\newline [The
object \texttt{t(num0)} is a matrix of character.  It has to be
character because its first row is character.  The command
\texttt{print(t(num0), quote=FALSE))} can be used to print it without
the quotes.]

\section*{Sweave markup}

To ensure that comments are retained in your code, place the
following in your document preamble (actually it does not have to be
in the preamble).
\begin{verbatim}
\\SweaveOpts{keep.source=TRUE}
\end{verbatim}

\paragraph{Setting options globally:}  The following is an example.

\begin{figure}[H]
  \fbox{
  \includegraphics{setup.pdf}
  }
\end{figure}

This reduces the margins around graphs (\texttt{mar=c(4.1,3.6,1.6,1.1)}
created using base graphics, gives plots that are square (\texttt{pty="s"}),
moves the axis labels and tick labels in somewhat (\texttt{mgp=c(2.25,0.5,0)}),
and reduces the length of tick marks (\texttt{tck=-0.025}).  The call to
\texttt{pdf.options()} sets the pointsize for postscript and pdf graphs.
Note that these settings have no effect for {\em lattice} and {\em ggplot2} 
graphs.

Also possible are \texttt{options(prompt="  ", continue="    ")},
or \texttt{options(continue="    ")}
Use the first of these options only once everything is working.  
Also, include near the end of the document:

\begin{figure}[H]
   \fbox{ 
  \includegraphics{prompt.pdf}
  }
\end{figure}

To restore the prompt and continuation prompt from the command line, perhaps
because the document has not processed properly through Sweave, type:
\begin{verbatim}
options(prompt="> ", continue="+  ")
\end{verbatim}

\subsection*{Create named code fragment, use, execute}
For example

\begin{figure}[H]
   \fbox{ 
  \includegraphics{namedfrag.pdf}
  }
\end{figure}

\begin{figure}[H]
   \fbox{ 
  \includegraphics{execfrag.pdf}
  }
\end{figure}

\begin{figure}[H]
   \fbox{ 
  \includegraphics{echofrag.pdf}
  }
\end{figure}

Use the argument \texttt{results=hide}, if you wish to hide
any results from the calculation.

\subsection*{Graphs}

Here is an example
\begin{figure}[H]
   \fbox{ 
  \includegraphics{histfrag.pdf}
  }
\end{figure}

\begin{figure}[H]
   \fbox{ 
  \includegraphics{creategph.pdf}
  }
\end{figure}

Note that use of \verb!@ %! rather than \verb!@! is optional.  The \verb!%!
tells editors to respect the line feed when reformatting.


\end{document}
