\documentclass[11pt]{article}
\def\Z{\mathbb{Z}}
\def\Q{\mathbb{Q}}
\def\R{\mathbb{R}}
\def\C{\mathbb{C}}
\usepackage{latexsym,amsmath,amssymb}
\usepackage[all]{xy}
\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{0in}
\setlength{\textheight}{8.0in}
\begin{document}
\centerline{Mathematics 375 -- Probability Theory}
\centerline{Computer Lab Day 1 -- Introducing R}
\centerline{September 5, 2011}
\vskip 10pt
\noindent
{\it Goals}
\vskip 10pt
The goals of today's lab are:
\begin{itemize}
\item to introduce some basic features of the R statistics 
package that we will be using some this semester and much
more heavily next semester, and 
\item to introduce some additional descriptive statistics 
and graphics for understanding the ``shape'' or ``distribution'' 
of a data set.
\end{itemize}
\vskip 10pt
\noindent
{\it Background on R -- Getting Started}
\vskip 10pt
The instructions here are geared toward the Windows version
of R running in the Haberlin 136 lab.  R is also installed
on the Math/CS Linux network (Swords 219 lab, other machines
on Swords 3).  Since R is a free,
open source package, you can also obtain versions to run
on your own computers if you want, and there are Windows, 
Mac-OS, and Linux versions available.  These have the same 
functionality but slightly different interfaces.   

In HA 136, when you double-click the R shortcut from the
standard desktop, you will launch the R-GUI window.
Inside this is a workspace window with a white background.
This is \emph{not} entirely like a Maple worksheet, unfortunately.
It \emph{only} allows you to 
\begin{itemize}
\item enter commands, and 
\item view numerical output.
\end{itemize}
It does not let you create integrated documents, and
graphics are displayed in separate windows within the 
R-GUI, as you will see.

You can \emph{print} the workspace window or any graphics
window at any time.  Just highlight that window and
press the toolbar printer button.  

To save your work, with comments and answers to questions
(e.g. for the lab reports you submit for the lab day 
assignments), I recommend 
\begin{itemize}
\item opening another editor window (either WordPad or NotePad
are fine for this on Windows systems), and
\item using it to create a text file that will be
your permanent record of your work.
\item When you have a result you want to save from
the R workspace, copy both the input command and
the output and paste into the text file.
\item Any time you want to add a comment to explain
a calculation or answer a question, enter that 
directly into the text file. 
\item Save that text file on your personal network
$P:$ drive and print it when you want a hardcopy
(from WordPad or NotePad).  
\end{itemize}
\vskip 10pt
\noindent
{\it A First R Session}
\vskip 10pt
In the lab today, you will basically work through some 
examples illustrating features of R and some of the 
descriptive statistics we have talked about so far in 
the course.  You will keep a record of what you do and
submit it as the first lab assignment.
\vskip 10pt
{\it Lab Exercises}
\vskip 10pt
\begin{enumerate} 
\item In its most basic form, R provides a command-driven
``statistical calculator'' in which you type in data and
commands to generate numerical and graphical output that
are displayed.  When you have the R window open, note the
$>$ input prompt after the opening message.  This is 
where you should start.  This exercise works through
problem 1.35 in our text (Wackerly, Mendenhall and Scheaffer).
First we need to get the systolic blood pressure readings
into R so that we can work with them.  The basic way is to 
enter (after the $>$ prompt):
\vskip 10pt
\begin{center}
{\tt bpdata <- c(172,140,123,130,115,148,108,129,137,161,123,152,133,128,142)}
\end{center}
\vskip 10pt
\noindent
and press ENTER.  Note:  the $<-$ is formed by typing the $<$ and $-$
characters right next to each other.  It is the {\it assignment operator}
in R.  This takes the data on the right and assigns it to the name on
the left.   The $c$, followed by the list
of numbers, separated by commas is R's way of constructing 
{\it an ordered list or vector}.  You should just get another
input prompt if everything went OK with this (error messages if not).
If you did get an error, try again, and try to follow the above exactly.
You can verify that the numbers are entered correctly by typing in the 
name {\tt bpdata} as a command; the numbers will be displayed in a row.
\begin{enumerate}
\item[a.]  Part a of the problem says to estimate the SD of the 
data using the range (max minus min).  The empirical rule we discussed
says that a \emph{very rough} estimate for the SD is 
$$SD \doteq \frac{1}{4}(range).$$
(This tends to work better for \emph{small} data sets without ``outliers.'')
In R we can do this computation very easily using built-in functions:
Type in 
\vskip 10pt
\begin{center}
{\tt (max(bpdata) - min(bpdata))/4}
\end{center}
\vskip 10pt
\noindent
and press ENTER to see the result.
\item[b.]  Part b of the problem says to compute the actual mean $\overline{y}$ and
$s$ of the data.  This can also be done in R (in several different ways!):
Type in 
\vskip 10pt
\begin{center}
{\tt ybar <- sum(bpdata)/length(bpdata)}\\
{\tt mean(bpdata)}\\
{\tt s <- sqrt(sum((bpdata - ybar)\^{}2)/(length(bpdata) - 1))}\\
{\tt sd(bpdata)}
\end{center}
\vskip 10pt
\noindent
and press ENTER after each line to see the results.  Note the way the sum function
works in a very nice way here -- in the computation of {\tt ybar}, we only need
to put in the name of the data list; the sum of all the elements is computed by 
default.  Then in the computation of {\tt s}, we can subtract \emph{each entry}
of the {\tt bpdata} vector from {\tt ybar} as shown!  Of course we will usually
want now to use the built-in functions {\tt mean} and {\tt sd}, but
the first of each pair of computations shows how we can also use more primitive functions
in R when we need them.  
\item[c.]  Tchebysheff's theorem is a result about the distribution of data 
sets that complements the ``empirical rule.''  It says in a particular case
that \emph{there will always be at least $75\%$ of the numbers in a data set 
within 2 standard deviations of the mean}.
We will prove this and more general statements later in the semester.  For 
now, let's use R to check that this is true here.  In its most basic
use, R can function as a simple numerical calculator(!)
\vskip 10pt
\begin{center}
{\tt upperb <- ybar + 2*s}\\
{\tt lowerb <- ybar - 2*s}
\end{center}
\vskip 10pt
\noindent
How many of the data values are in the range $[\overline{y} - 2 s, \overline{y} + 2s]$,
though?  We could just count by hand, of course, with a data set this small.
But R also contains a very rich collection of facilities for selecting, subsetting, 
and massaging data in numerous ways.  Here is a basic way to select the portion 
of the data set in the range we want:
\vskip 10pt
\begin{center}
{\tt middle <- subset(bpdata, bpdata <= upperb \& bpdata  >= lowerb)}\\
{\tt length(middle)/length(bpdata)}
\end{center}
\vskip 10pt
(This should be pretty self-explanatory; ask if you don't understand what this did!)
Did Tchebysheff's theorem ``work'' here?
\noindent
\end{enumerate}
\item Now let's look at a different data set.  Problem 1.36 in our text
deals with the number of parasites found in each of a sample of 100 individual foxes.
The information is given there in table format:
\begin{center}
\begin{tabular}{l|ccccccccc}
Parasites & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8\\ \hline
Foxes & 69 & 17 & 6 & 3 & 1 & 2 & 1 & 0 & 1\\
\end{tabular}
\end{center}
So we want to create a data set containing 69 0 entries, 17 1 entries, etc.
Clearly this is going to be tedious and error prone if we do it the most basic
way(!)  Fortunately, R lets us construct lists in more flexible ways
besides just listing the entries.  We can indicate how many times to repeat
a given value if we want, like this:
\vskip 10pt
\begin{center}
{\tt foxpar <- rep(c(0,1,2,3,4,5,6,7,8),c(69,17,6,3,1,2,1,0,1))}\\
\end{center}
\vskip 10pt
\noindent
Of course {\tt rep} stands for ``repetition'' here.
Note that the first list inside the {\tt rep} gives the values of the number of parasites -- the first row in the table format -- and the 
second list gives number of times each value appears in the data -- the second row.  The first list could also be abbreviated as
{\tt 0:8} (R's notation for a list of successive integer values).  
\begin{enumerate}
\item[a.]  Part a of the problem says to construct a relative frequency histogram
for the number of parasites per fox.  To generate this in R, let's begin by entering
\vskip 10pt
\begin{center}
{\tt hist(foxpar)}\\
\end{center}
\vskip 10pt
\noindent
The histogram will be generated in a separate graphics window.  
This is a ``default'' version of the command.  The {\it hist}
command takes a number of options that control how the ``bins''
are selected and how the histogram is drawn as well.  You 
will often want to modify the default results.   For
instance, looking at the output, how were the bins chosen here?
Is that optimal for showing the different numbers of foxes
with each number of parasites?  Suppose we want bins centered
on the integer values $0,1,2,3,4,5,6,7,8$ as we did in another
example in class.  We can put in an option to make
R do it that way too.  The {\tt breaks} option takes
a list of values and makes those the boundaries of the bins 
for the histogram:
\vskip 10pt
\begin{center}
{\tt hist(foxpar,breaks=c(-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5))}\\
\end{center}
\vskip 10pt
\noindent
How is this different and why is it better than the default
way that R drew the histogram?  Include a printout of this
histogram in your lab report.  (See above for printing directions.)
\item[b.] To practice using the commands we talked about
in problem 1, calculate $\overline{y}$ and $s$ for this data set.
How many of the numbers were within two standard deviations
of the mean in this case?  
\item[c.]  The list of breaks for the histogram can be
anything in R.  What happens if you change the command above
to 
\vskip 10pt
\begin{center}
{\tt hist(foxpar,breaks=c(-0.5,0.5,1.5,2.5,3.5,4.5,5.5,8.5))}\\
\end{center}
\vskip 10pt
\noindent
(deleting the 6.5 and 7.5 break points, so we're lumping
together the counts for the values 6,7,8 parasites)? 
What is different about the new histogram that R drew in this case?  
Look at the graph carefully: do the heights of the boxes still represent 
\emph{frequencies}?  (Hint: what
they do represent are \emph{densities} that are computed
by taking the frequency count and dividing by the 
width of the box in each case.  What is true if we 
add the areas of the boxes drawn this way?  We will see
other reasons why this is a preferable way to draw histograms
with unequal bin sizes later. )
If we \emph{really want frequencies}, we can make R do it that way:
\vskip 10pt
\begin{center}
{\tt hist(foxpar,breaks=c(-0.5,0.5,1.5,2.5,3.5,4.5,5.5,8.5),freq=TRUE)}\\
\end{center}
\vskip 10pt
\noindent
However note that now R generates a warning message saying
the AREAS in the plot are ``wrong.''  Think about this.
Why is the last box in this version of the histogram
\emph{misleading} if drawn this way?  
\end{enumerate}
\item  In almost all real-world analysis of data, the data
entry and analysis steps are separate. So, for instance,
someone might enter the data into a computer file in some
specified format.  Then a statistician would read
that data into a program like R to analyze the data and
look for patterns or try to draw inferences about the 
situation described by those measurements.  We will see
an example of that in this question.  Enter this command
to read in a modified version of the data set about lifetimes
of radio transmitter-receivers from 
problem 1.25 in the text.  
\vskip 10pt
\begin{center}
{\tt lifetimes <- read.table("http://mathcs.holycross.edu/\~{}little/\\
                  ProbStat1112/R/Lab1Data.dat",header=TRUE,sep=",")}\\
{\tt attach(lifetimes)}
\\
\end{center}
(Note: the first two lines here are part of one command, just keep typing until
the final close paren.)
\vskip 10pt
\noindent
The modification is that the numbers have been explicitly 
grouped into 8 columns of 11 numbers each, with a column
heading that is a character string.  The {\tt attach} 
command above sets up the headings as names for the different 
columns.  For the 
purposes of this question, we will think of them as 8 separate
samples of lifetimes, called {\tt s1, ..., s8}.   Each of those
8 data sets will allow for the same operations as the data sets from problems
1 and 2 above.  
\begin{enumerate}
\item[a.]  In addition to the mean and standard deviation, statisticians
use the median and other percentiles to understand how data sets
are distributed.  R has multiple ways to compute these, including
a default {\tt summary} command that computes the min, max,  median, mean, 
25th and 75th percentiles all together:
\vskip 10pt
\begin{center}
{\tt summary(s2)}\\
{\tt summary(s8)}
\\
\end{center}
\vskip 10pt
\noindent
Note that the 25th percentile is called the 1st quartile and the 
75th percentile is called the 3rd quartile in R.  What is your first
impression about the distributions of these two samples, given
these summary statistics?
\item[b.]  Plot the histograms for each of {\tt s2, s8} to see if
your impression is confirmed.
\item[c.]  There are other very useful graphical ways to represent data 
besides the frequency or density histogram.  For example, statisticians
may use what are called ``box and whisker'' plots for this purpose.  
Try 
\vskip 10pt
\begin{center}
{\tt boxplot(s2)}
\end{center}
\vskip 10pt
\noindent
to see one of these plots.  You can figure out what the plot represents
if you carefully look at the vertical scale together with the values
returned by the {\tt summary} command for this data set.  What is the 
``box;'' what are the ``whiskers?''  Can you guess what the circles represent?
\item[d.]  Often we want to compare the distributions of two or more
data sets.  We can do ``parallel box and whisker plots'' in R like this:
\vskip 10pt
\begin{center}
{\tt par(mfrow=c(1,2))}\\
{\tt boxplot(s2)}\\
{\tt boxplot(s8)}\\
{\tt par(mfrow=c(1,1))}\\
\end{center}
\vskip 10pt
\noindent
(Note:  The {\tt par} commands set up the graphing so that the box plots
come out side-by-side, then return the graphics format to the original
state for any graphing that comes later.)   Hand in a printout of this
plot with your lab report.
\item[e.]  Describe in words your conclusions about the distributions
of these two data sets.  (Comment:  There are apparent differences and
your job is to try to describe them.  Whether those are \emph{significant
enough} to justify a conclusion about a possible difference between
the two samples is the kind of question the statistical techniques
we will learn next semester are designed to help answer.)

\end{enumerate}
\end{enumerate}
\vskip 10pt
\noindent
{\it Assignment}
\vskip 10pt
Lab reports containing all input, output, and answers to the questions
posed above are due in class no later than Monday, September 12.  If you do not
finish during the hour today, you can return
to HA 136 any time it is not in use by another class.
 

\end{document}
