\documentclass{article}
\usepackage{graphicx}
\oddsidemargin=0in \evensidemargin=0in \topmargin=-.2in
\textwidth=6.5in \textheight=8.5in
\pagestyle{headings}
\parskip=1ex
\begin{document}
\title{Notes on the use of R for psychology experiments and
questionnaires}
\author{Jonathan Baron\\
Department of Psychology, University of Pennsylvania \and
Yuelin Li\\
Department of Psychiatry \& Behavioral Sciences\\
Memorial Sloan-Kettering Cancer Center
\thanks{Copyright \copyright 2000, Jonathan Baron
and Yuelin Li. {\scriptsize Permission is granted to make
and distribute verbatim copies of this document provided
the copyright notice and this permission notice are
preserved on all copies. For other permissions, please
contact the first author at {\tt
baron@cattell.psych.upenn.edu.}} We thank Andrew Hochman,
Rashid Nassar, Christophe Pallier, and Hans-Rudolf Roth for
helpful comments.}}
\date{March 6, 2011}
\maketitle
\tableofcontents
\section{Introduction}
This is a set of notes and annotated examples of the use of the
statistical package R. It is ``for psychology experiments and
questionnaires'' because we cover the main statistical methods
used by psychologists who do research on human subjects, but of
course it this is also relevant to researchers in others fields
that do similar kinds of research.
R, like S--Plus, is based on the S language invented at Bell
Labs. Most of this should also work with S--Plus. Because R is
open-source (hence also free), it has benefitted from the work of
many contributors and bug finders. R is a complete package. You
can do with it whatever you can do with Systat, SPSS, Stata, or
SAS, including graphics. Contributed packages are added or
updated almost weekly; in some cases these are at the cutting
edge of statistical practice.
Some things are more difficult with R, especially if you are used
to using menus. With R, it helps to have a list of commands in
front of you. There are lists in the on-line help and in the
index of \textit{An introduction to R} by the R Core Development
Team, and in the reference cards listed in
\verb"http://finzi.psych.upenn.edu/".
Some things turn out to be easier in R. Although there are no
menus, the on-line help files are very easy to use, and quite
complete. The elegance of the language helps too, particularly
those tasks involving the manipulation of data.
The purpose of this document is to reduce the difficulty of the
things that are more difficult at first. We assume that you have
read the relevant parts of \textit{An introduction to R}, but we
do not assume that you have mastered its contents. We assume that
you have gotten to the point of installing R and trying a couple
of examples.
\section{A few useful concepts and commands}
\subsection{Concepts}
In R, most commands are functions. That is, the command is
written as the name of the function, followed by parentheses,
with the arguments of the function in parentheses, separated by
commas when there is more than one, e.g., {\tt plot(mydata1)}.
When there is no argument, the parentheses are still needed,
e.g., {\tt q()} to exit the program.
In this document, we use names such as {\tt x1} or {\tt file1},
that is, names containing both letters and a digit, to indicate
variable names that the user makes up. Really these can be of
any form. We use the number simply to clarify the distinction
between a made up name and a key word with a pre-determined
meaning in R. R is case sensitive.
Although most commands are functions with the arguments in
parentheses, some arguments require specification of a key word
with an equal sign and a value for that key word, such as {\tt
source("myfile1.R",echo=T)}, which means read in {\tt
myfile1.R} and echo the commands on the screen. Key words can
be abbreviated (e.g., \texttt{e=T}).
In addition to the idea of a function, R has \textit{objects} and
\textit{modes}. Objects are anything that you can give a name.
There are many different \textit{classes} of objects. The main
classes of interest here are {\it vector, matrix, factor, list},
and {\it data frame}. The mode of an object tells what kind of
things are in it. The main modes of interest here are {\tt
logical, numeric,} and {\tt character}.
We sometimes indicate the class of object (vector, matrix, factor,
etc.) by using {\tt v1} for a vector, {\tt m1} for a matrix, and
so on. Most R functions, however, will either accept more than
one type of object or will ``coerce'' a type into the form that
it needs.
The most interesting object is a data frame. It is useful to
think about data frames in terms of rows and columns. The rows
are subjects or observations. The columns are variables, but a
matrix can be a column too. The variables can be of different
classes.
The behavior of any given function, such as {\tt plot(), aov()}
(analysis of variance) or {\tt summary()} depends on the object
class and mode to which it is applied. A nice thing about R is
that you almost don't need to know this, because the default
behavior of functions is usually what you want. One way to use R
is just to ignore completely the distinction among classes and
modes, but \textit{check} every step (by typing the name of the
object it creates or modifies). If you proceed this way, you
will also get error messages, which you must learn to interpret.
Most of the time, again, you can find the problem by looking at
the objects involved, one by one, typing the name of each object.
Sometimes, however, you must know the distinctions. For example,
a factor is treated differently from an ordinary vector in an
analysis of variance or regression. A factor is what is often
called a categorical variable. Even if numbers are used to
represent categories, they are not treated as ordered. If you
use a vector and think you are using a factor, you can be misled.
\subsection{Commands}
As a reminder, here is a list of some of the useful commands that
you should be familiar with, and some more advanced ones that are
worth knowing about. We discuss graphics in a later section.
\subsubsection{Getting help}
\noindent
{\tt help.start()} starts the browser version of the help files.
(But you can use {\tt help()} without it.) With a fast computer
and a good browser, it is often simpler to open the html
documents in a browser while you work and just use the browser's
capabilities.
\noindent
{\tt help(command1)} prints the help available about {\tt
command1}. {\tt help.search("keyword1")} searches keywords for
help on this topic.
\noindent
{\tt apropos(topic1)} or {\tt apropos("topic1")} finds commands
relevant to topic1, whatever it is.
\noindent
{\tt example(command1)} prints an example of the use of the
command. This is especially useful for graphics commands. Try,
for example, {\tt example(contour)}, {\tt example(dotchart)}, {\tt
example(image)}, and {\tt example(persp)}.
\subsubsection{Installing packages}
\label{sec:bioc}
\noindent
\texttt{install.packages(c("package1","package2"))} will install
these two packages from CRAN (the main archive), if your computer
is connected to the Internet. You don't need the \texttt{c()} if
you just want one package. You should, at some point, make sure
that you are using the CRAN mirror page that is closest to you.
For example, if you live in the U.S., you should have a
\texttt{.Rprofile} file with
\texttt{options(CRAN = "http://cran.us.r-project.org")} in it.
(It may work slightly differently on Windows.)
\noindent
\texttt{CRAN.packages(), installed.packages(), and
update.packages()} are also useful. The first tells you what
is available. The second tells you what is installed. The third
updates the packages that you have installed, to their latest
version.
To install packages from the Bioconductor set, see the
instructions in\\
\verb"http://www.bioconductor.org/reposToolsDesc.html".
When packages are not on CRAN, you can download them and use
\texttt{R CMD INSTALL package1.tar.gz} from a Unix/Linux command
line. (Again, this may be different on Windows.)
\subsubsection{Assignment, logic, and arithmetic}
\noindent
{\tt <-} assigns what is on the right of the arrow to what is on
the left. (If you use ESS, the \verb"_" key will produce this
arrow with spaces, a great convenience.)
\noindent
Typing the name of the object prints the object. For example, if
you say:
\begin{verbatim}
t1 <- c(1,2,3,4,5)
t1
\end{verbatim}
you will see {\tt 1 2 3 4 5}.
\noindent
Logical objects can be true or false. Some functions and
operators return TRUE or FALSE. For example, \texttt{1==1}, is
TRUE because 1 does equal 1. Likewise, \texttt{1==2} is FALSE,
and \texttt{1<2} is TRUE. Use \texttt{all()},
\verb"any(), |, ||, &", and \verb"&&" to combine logical expressions, and use
\texttt{!} to negate them. The difference between the \texttt{|}
and \texttt{||} form is that the shorter form, when applied to
vectors, etc., returns a vector, while the longer form stops when
the result is determined and returns a single TRUE or FALSE.
\noindent
Set functions operate on the elements of vectors:
\texttt{union(v1,v2), intersect(v1,v2), setdiff(v1,v2),
setequal(v1,v2), is.element(element1,v1)} (or, \texttt{element1
\%in\% v1}).
\noindent
Arithmetic works. For example, {\tt -t1} yields {\tt -1 -2 -3 -4
-5}. It works on matrices and data frames too. For example,
suppose {\tt m1} is the matrix
\begin{verbatim}
1 2 3
4 5 6
\end{verbatim}
Then {\tt m1 * 2} is
\begin{verbatim}
2 4 6
8 10 12
\end{verbatim}
\noindent
Matrix multiplication works too. Suppose {\tt m2} is the matrix
\begin{verbatim}
1 2
1 2
1 2
\end{verbatim}
then \verb"m1 %*% m2" is
\begin{verbatim}
6 12
15 30
\end{verbatim}
and \verb"m2 %*% m1" is
\begin{verbatim}
9 12 15
9 12 15
9 12 15
\end{verbatim}
You can also multiply a matrix by a vector using matrix
multiplication, vectors are aligned
vertically when they come after the {\verb"%*%"} sign and
horizontally when they come before it. This is a good way to
find weighted sums, as we shall explain.
For ordinary multiplication of a matrix times a vector, the
vector is vertical and is repeated as many times as needed. For
example {\tt m2 * 1:2} yields
\begin{verbatim}
1 4
2 2
1 4
\end{verbatim}
Ordinarily, you would multiply a matrix by a vector when the
length of the vector is equal to the number of rows in the
matrix.
\subsubsection{Vectors, matrices, lists, arrays, and data frames}
\noindent
{\tt :} is a way to abbreviate a sequence of numbers, e.g., {\tt
1:5} is equivalent to 1,2,3,4,5.
\noindent
{\tt c(number.list1)} makes the list of numbers (separated by
commas) into a vector object. For example, {\tt c(1,2,3,4,5)}
(but {\tt 1:5} is already a vector, so you do not need to say
{\tt c(1:5)}).
\noindent
{\tt rep(v1,n1)} repeats the vector {\tt v1 n1} times. For example,
{\tt rep(c(1:5),2)} is {\tt 1,2,3,4,5,1,2,3,4,5}.
\noindent
{\tt rep(v1,v2)} repeats each element of the vector {\tt v1} a
number of times indicated by the corresponding element of the
vector {\tt v2}. The vectors {\tt v1} and {\tt v2} must have the
same length. For example, {\tt rep(c(1,2,3),c(2,2,2))} is {\tt
1,1,2,2,3,3}. Notice that this can also be written as {\tt
rep(c(1,2,3),rep(2,3))}. (See also the function {\tt gl()} for
generating factors according to a pattern.)
\noindent
{\tt cbind(v1,v2,v3)} puts vectors {\tt v1, v2,} and {\tt v3}
(all of the same length) together as columns of a matrix. You
can of course give this a name, such as {\tt mat1 <-
cbind(v1,v2,v2)}.
\noindent
{\tt matrix(v1,rows1,colmuns1)} makes the vector {\tt v1} into a
matrix with the given number of rows and columns. You don't need
to specify both rows and columns, but you do need to put in both
commas. You can also use key words instead of using position to
indicate which argument is which, and then you do not need the
commas. For example, {\tt matrix(1:10, ncol=5)} represents the
matrix
\begin{verbatim}
1 3 5 7 9
2 4 6 8 10
\end{verbatim}
Notice that the matrix is filled column by column.
\noindent
{\tt data.frame(vector.list1)} takes a list of vectors, all of the
same length (error message if they aren't) and makes them into a
data frame. It can also include factors as well as vectors.
\noindent
{\tt dim(obj1)} prints the dimensions of a matrix, array or data frame.
\noindent
{\tt length(vector1)} prints the length of {\tt vector1}.
\noindent
You can refer to parts of objects. {\tt m1[,3]} is the third
column of matrix {\tt m1}. {\tt m1[,-3]} is all the columns
except the third. {\tt m1[m1[,1]>3,]} is all the rows for which
the first column is greater than 3. {\tt v1[2]} is the second
element of vector {\tt v1}. If {\tt df1} is a data frame with
columns {\tt a, b}, and {\tt c}, you can refer to the third
column as \verb+df1$c+. % $
\noindent
Most functions return lists. You can see the elements of a list
with \texttt{unlist()}. For example, try
\texttt{unlist(t.test(1:5))} to see what the \texttt{t.test()}
function returns. This is also listed in the section of help
pages called ``Value.''
\noindent
\texttt{array()} seems very complicated at first, but it is
extremely useful when you have a three-way classification, e.g.,
subjects, cases, and questions, with each question asked about
each case. We give an example later.
\noindent
\texttt{outer(m1,m2,"fun1")} applies \texttt{fun1}, a function of
two variables, to each combination of \texttt{m1} and
\texttt{m2}. The default is to multiply them.
\noindent
\texttt{mapply("fun1",o1,o2)}, another very powerful function,
applies \texttt{fun1} to the elements of \texttt{o1} and
\texttt{o2}. For example, if these are data frames, and
\texttt{fun1} is \texttt{"t.test"}, you will get a list of t
tests comparing the first column of \texttt{o1} with the first
column of \texttt{o2}, the second with the second, and so on.
This is because the basic elements of a data frame are the
columns.
\subsubsection{String functions}
R is not intended as a language for manipulating text, but it is
surprisingly powerful. If you know R, you might not need to
learn Perl. Strings are character variables that consist of
letters, numbers, and symbols.
\texttt{strsplit()} splits a string, and \texttt{paste()} puts a
string together out of components.
\texttt{grep(), sub(), gsub(),} and \texttt{regexpr()} allow you
to search for, and replace, parts of strings.
The set functions such as \texttt{union(), intersect(),
setdiff()}, and \texttt{\%in\%} are also useful for dealing
with databases that consist of strings such as names and email
addresses.
You can even use these functions to write new R commands as
strings, so that R can program itself! Just to see an example of
how this works, try \texttt{eval(parse(text="t.test(1:5)"))}.
The \texttt{parse()} function turns the text into an expression,
and \texttt{eval()} evaluates the expression. So this is
equivalent to \texttt{t.test(1:5)}. But you could replace
\texttt{t.test(1:5)} with any string constructed by R itself.
\subsubsection{Loading and saving}
\noindent
{\tt library(xx1)} loads the extra library. A useful library for
psychology is and {\tt mva} (multivariate analysis). To find the
contents of a library such as {\tt mva} before you load it, say
{\tt library(help=mva)}. The {\tt ctest} library is already
loaded when you start R.
\noindent
{\tt source("file1")} runs the commands in {\tt file1}.
\noindent
{\tt sink("file1")} diverts output to {\tt file1} until you say {\tt
sink()}.
\noindent
{\tt save(x1,file="file1")} saves object {\tt x1} to file {\tt file1}.
To read in the file, use {\tt load("file1")}.
\noindent
{\tt q()} quits the program. \verb'q("yes")' saves everything.
\noindent
{\tt write(object, "file1")} writes a matrix or some other object to
{\tt file1}.
\noindent
{\tt write.table(object1,"file1")} writes a table and has an option
to make it comma delimited, so that (for example) Excel can read
it. See the help file, but to make it comma delimited, say\\
{\tt write.table(object1,"file1",sep=",")}.
\noindent
{\tt round()} produces output rounded off, which is useful when you
are cutting and pasting R output into a manuscript. For example,
\verb"round(t.test(v1)$statistic,2)" rounds off the value of {\tt
t} to two places. Other useful functions are {\tt format} and
{\tt formatC}. For example, if we assign {\tt t1 <- t.test(v1},
then the following command prints out a nicely formatted result,
suitable for dumping into a paper:
\begin{verbatim}
print(paste("(t_{",t1[[2]],"}=",formatC(t1[[1]],format="f",digits=2),
", p=",formatC(t1[[3]],format="f"),")",sep=""),quote=FALSE)
\end{verbatim}%$
This works because {\tt t1} is actually a list, and the numbers
in the double brackets refer to the elements of the list.
\noindent
{\tt read.table("file1")} reads in data from a file. The first line
of the file can (but need not) contain the names of the variables
in each column.
\subsubsection{Dealing with objects}
\noindent
{\tt ls()} lists all the active objects.
\noindent
{\tt rm(object1)} removes object1. To remove all objects, say
{\tt rm(list=ls())}.
\noindent
{\tt attach(data.frame1)} makes the variables in {\tt
data.frame1} active and available generally.
\noindent
\texttt{names(obj1)} prints the names, e.g., of a matrix or data
frame.
\noindent
\texttt{typeof()}, \texttt{mode())}, and \texttt{class()} tell
you about the properties of an object.
\subsubsection{Summaries and calculations by row, column, or
group}
\noindent
{\tt summary(x1)} prints statistics for the variables
(columns) in {\tt x1}, which may be a vector, matrix, or data
frame. See also the {\tt str()} function, which is similar, and
\texttt{aggregate()}, which summarizes by groups.
\noindent
{\tt table(x1)} prints a table of the number of times each value
occurs in {\tt x1}. {\tt table(x1,y1)} prints a cross-tabulation
of the two variables. The {\tt table} function can do a lot
more. Use {\tt prop.table()} when you want proportions rather
than counts.
\noindent
{\tt ave(v1,v2)} yields averages of vector {\tt v1} grouped by
the factor {\tt v2}.
\noindent
{\tt cumsum(v1)} is the cumulative sum of vector {\tt v1}.
\noindent
You can do calculations on rows or columns of a matrix and get
the result as a vector. {\tt apply(x1,2,mean)} yields just the
means of the columns. Use {\tt apply(x1,1,mean)} for the rows.
You can use other functions aside from {\tt mean}, such as {\tt
sd}, {\tt max}, {\tt min} or {\tt sum}. To ignore missing
data, use {\tt apply(x1,2,mean,na.rm=T)}, etc. For sums and
means, it is easier to use \texttt{rowSums(), colSums(),
rowMeans(),} and \texttt{colMeans} instead of
\texttt{apply()}. Note that you can use apply with a function,
e.g., \texttt{apply(x1,1,function(x) exp(sum(log(x)))} (which is
a roundabout way to write \texttt{apply(x1,1,prod)}). The same
thing can be written in two steps, e.g.:
\begin{verbatim}
newprod <- function(x) {exp(sum(log(x)))
apply(x1,1,newprod)
\end{verbatim}
You can refer to a subset of an object in many other ways. One
way is to use a square bracket at the end, e.g.,
\texttt{matrix1[,1:5]} refers to columns 1 through 5 of the
matrix. You can also use this method for new objects, e.g.,
\texttt{(matrix1+matrix2)[,1:5]}, which refers to the first five
columns of the sum of the two matrices. Another important method
is the use of \texttt{by()} or \texttt{aggregate()} to compute
statistics for subgroups defined by vectors or factors. You can
also use \texttt{split()} to get a list of subgroups. Finally,
many functions allow you to use a \texttt{subset} argument.
\subsubsection{Functions and debugging}
\noindent
\texttt{function()} allows you to write your own functions.
\noindent
Several functions are useful for debugging your own functions or
scripts: \texttt{traceback(), debug(), browser(), recover()}.
\section{Basic method}
The following basic method is assumed here. You have a command
file and then submit it, for each data set. Thus, for each
experiment or study, you have two files. One consists of the
data. Call it {\tt exp1.data}. The other is a list of commands
to be executed by R, {\tt exp1.R}. (Any suffixes will do,
although ESS recognizes the R suffix and loads special features
for editing.) The advantage of this approach is that you have a
complete record of what your transformed variables mean. If your
data set is small relative to the speed of your computer, it is a
good idea to revise exp1.R and re-run it each time you make a
change that you want to keep. So you could have exp1.R open in
the window of an editor while you have R in another
window.\footnote{If you want, you can put your data in the same
file as the commands. The simplest way to do this is to put
commas between the numbers and then use a command like {\tt t1
<- c(1,2,3, \dots 26,90)}, possibly over several lines, where
the numbers are your data.}
To analyze a data set, you start R in the directory where the
data and command file are. Then, at the R prompt, you type
\begin{verbatim}
source("exp1.R")
\end{verbatim}
and the command file runs. The first line of the command file
usually reads in the data. You may include statistics and
graphics commands in the source file. You will not see the
output if you say \verb'source("data1.R")', although they will
still run. If you want to see the output, say
\begin{verbatim}
source("data1.R",echo=T)
\end{verbatim}
Command files can and should be annotated. R ignores everything
after a \#. In this document, the examples are not meant to be
run.
We have mentioned ESS, which stands for ``Emacs Speaks
Statistics.'' This is an add-on for the Emacs editor, making
Emacs more useful with several different statistical programs,
including R, S--Plus, and SAS. \footnote{ESS is wonderful, but
Emacs will cause trouble for you if you use a word processor
like Word, and if you are used to shortcut keys such as ctrl-x
for cutting text. The shortcut keys in Emacs are all
different, and this leads to serious mind-boggle. One
solution, adopted by the first author, is to give up word
processors and editors that use the same shortcuts, such as
Winedt. A side effect of this solution is that you have very
little reason to use Microsoft Windows.} If you use ESS, then
you will want to run R as a process in Emacs, so, to start R, say
{\tt emacs -f R}. You will want {\tt exp1.R} in another window,
so also say {\tt emacs exp1.R}. With ESS, you can easily cut and
paste blocks (or lines) of commands from one window to another.
Here are some tips for debugging:
\begin{itemize}
\item If you use the \verb'source("exp1.R")' method described here,
use \verb'source("exp1.R",echo=T)' to echo the input and see
how far the commands get before they ``bomb.''
\item Use {\tt ls()} to see which objects have been created.
\item Often the problem is with a particular function, often
because it has been applied to the wrong type or size of
object. Check the sizes of objects with {\tt dim()} or (for
vectors) {\tt length()}.
\item Look at the {\tt help()} for the function in question. (If
you use {\tt help.start()} at the beginning, the output will
appear in your browser. The main advantage of this is that you
can follow links to related functions very easily.)
\item Type the names of the objects to make sure they are what
you think they are.
\item If the help is not helpful enough, make up a little example
and try it. For example, you can get a matrix by saying {\tt
m1 <- matrix(1:12,,3)}.
\item For debugging functions, try {\tt debug()}, {\tt
browser()}, and {\tt traceback()}. (See their help pages.
We do very little with functions here.)
\end{itemize}
\section{Reading and transforming data}
\subsection{Data layout}
R, like Splus and S, represents an entire conceptual system for
thinking about data. You may need to learn some new ways of
thinking. One way that is new for users of Systat in particular
(but perhaps more familiar to users of SAS) concerns two
different ways of laying out a data set. In the Systat way, each
subject is a row (which may be continued on the next row if too
long, but still conceptually a row) and each variable is a column.
You can do this in R too, and most of the time it is sufficient.
But some the features of R will not work with this kind of
representation, in particular, repeated-measures analysis of
variance. So you need a second way of representing data, which
is that each row represents a single datum, e.g., one subject's
answer to one question. The row also contains an identifier for
all the relevant classifications, such as the question number,
the subscale that the question is part of, AND the subject.
Thus, ``subject'' becomes a category with no special status,
technically a factor (and remember to make sure it is a factor,
lest you find yourself studying the effect of the subject's
number).
\subsection{A simple questionnaire example}
Let us start with an example of the old-fashioned way. In the
file {\tt ctest3.data}, each subject is a row, and there are 134
columns. The first four are age, sex, student status, and time
to complete the study. The rest are the responses to four
questions about each of 32 cases. Each group of four is preceded
by the trial order, but this is ignored for now.
\begin{verbatim}
c0 <- read.table("ctest3.data")
\end{verbatim}
The data file has no labels, so we can read it with {\tt
read.table}.
\begin{verbatim}
age1 <- c0[,1]
sex1 <- c0[,2]
student1 <- c0[,3]
time1 <- c0[,4]
nsub1 <- nrow(c0)
\end{verbatim}
We can refer to elements of {\tt c0} by {\tt c0[row,column]}.
For example, {\tt c0[1,2]} is the sex of the first subject. We
can leave one part blank and get all of it, e.g., {\tt c0[,2]} is
a vector (column of numbers) representing the sex of all the
subjects. The last line defines {\tt nsub1} as the number of
subjects.
\begin{verbatim}
c1 <- as.matrix(c0[,4+1:128])
\end{verbatim}
Now {\tt c1} is the main part of the data, the matrix of responses.
The expression 1:128 is a vector, which expands to 1 2 3 \dots 128.
By adding 4, it becomes 5 6 7 \dots 132.
\subsubsection{Extracting subsets of data}
\begin{verbatim}
rsp1 <- c1[,4*c(1:32)-2]
rsp2 <- c1[,4*c(1:32)-1]
\end{verbatim}
The above two lines illustrate the extraction of sub-matrices
representing answers to two of the four questions making up each
item. The matrix {\tt rsp1} has 32 columns, corresponding to
columns 2 6 10 ... 126 of the original 128-column matrix {\tt
c1}. The matrix {\tt rsp2} corresponds to 3 7 11 ...
127.
Another way to do this is to use an array. We could say {\tt a1
<- array(c1,c(ns,4,32))}. Then {\tt a1[,1,]} is the equivalent
of {\tt rsp1}, and {\tt a1[20,1,]} is {\tt rsp1} for subject 20.
To see how arrays print out, try the following:
\begin{verbatim}
m1 <- matrix(1:60,5,)
a1 <- array(m1,c(5,2,6))
m1
a1
\end{verbatim}
You will see that the rows of each table are the first index and
the columns are the second index. Arrays seem difficult at
first, but they are very useful for this sort of analysis.
\subsubsection{Finding means (or other things) of sets of variables}
\begin{verbatim}
r1mean <- apply(rsp1,1,mean)
r2mean <- apply(rsp2,1,mean)
\end{verbatim}
The above lines illustrate the use of {\tt apply} for getting
means of subscales. In particular, {\tt abrmean} is the mean of
the subscale consisting of the answers to the second question in
each group. The {\tt apply} function works on the data in its
first argument, then applies the function in its third argument,
which, in this case, is {\tt mean}. (It can be {\tt max} or {\tt
min} or any defined function.) The second argument is 1 for
rows, 2 for columns (and so on, for arrays). We want the
function applied to rows.
\begin{verbatim}
r4mean <- apply(c1[,4*c(1:32)],1,mean)
\end{verbatim}
The expression here represents the matrix for the last item in
each group of four. The first argument can be any matrix or data
frame. (The output for a data frame will be labeled with row or
column names.) For example, suppose you have a list of variables
such as {\tt q1, q2, q3}, etc. Each is a vector, whose length is
the number of subjects. The average of the first three variables
for each subject is, {\tt apply(cbind(q1,q2,q3),1,mean)}. (This
is the equivalent of the Systat expression {\tt avg(q1,q2,q3)}.
A little more verbose, to be sure, but much more flexible.)
You can use apply is to tabulate the values of each column of a
matrix {\tt m1}: {\tt apply(m1,2,table)}. Or, to find column
means, {\tt apply(m1,2,mean)}.
There are many other ways to make tables. Some of the relevant
functions are {\tt table}, {\tt tapply}, {\tt sapply}, {\tt ave},
and {\tt by}. Here is an illustration of the use of {\tt by}.
Suppose you have a matrix {\tt m1} like this:
\begin{verbatim}
1 2 3 4
4 4 5 5
5 6 4 5
\end{verbatim}
The columns represent the combination of two variables, {\tt y1}
is {\tt 0 0 1 1}, for the four columns, respectively, and {\tt
y2} is {\tt 0 1 0 1}. To get the means of the columns for the
two values of {\tt y1}, say {\tt by(t(m1), y1, mean)}. You get
3.67 and 4.33 (labeled appropriately by values of y1). You need
to use {\tt t(m1)} because {\tt by} works by rows. If you say
{\tt by(t(m1), data.frame(y1,y2), mean)}, you get a cross
tabulation of the means by both factors. (This is, of course,
the means of the four columns of the original matrix.)
Of course, you can also use {\tt by} to classify rows; in the
usual examples, this would be groups of subjects rather than
classifications of variables.
\subsubsection{One row per observation}
The next subsection shows how to transform the data from a layout
from ``one row per subject'' to ``one row per observation.''
We're going to use the matrix {\tt rsp1}, which has 32 columns
and one row per subject. Here are five subjects:
\begin{verbatim}
1 1 2 2 1 2 3 5 2 3 2 4 2 5 7 7 6 6 7 5 7 8 7 9 8 8 9 9 8 9 9 9
1 2 3 2 1 3 2 3 2 3 2 3 2 3 2 4 1 2 4 5 4 5 5 6 5 6 6 7 6 7 7 8
1 1 2 3 1 2 3 4 2 3 3 4 2 4 3 4 4 4 5 5 5 6 6 7 6 7 7 8 7 7 8 8
1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6 7 7 7 7 8 8 8 8 9 9 9
1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6 7 7 7 7 8 8 8 8 9
\end{verbatim}
We'll create a matrix with one row per observation. The first
column will contain the observations, one variable at a time, and
the remaining columns will contain numbers representing the
subject and the level of the observation on each variable of
interest. There are two such variables here, {\tt r2} and {\tt
r1}. The variable {\tt r2} has four levels, {\tt 1 2 3 4}, and
it cycles through the 32 columns as {\tt 1 2 3 4 1 2 3 4 ...}
The variable {\tt r1} has the values (for successive columns)
{\tt 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 1 1 1 1 2 2 2 2 3 3 3 3 4 4
4 4}. These levels are ordered. They are not just arbitrary
labels. (For that, we would need the {\tt factor} function.)
\begin{verbatim}
r2 <- rep(1:4,8)
r1 <- rep(rep(1:4,rep(4,4)),2)
\end{verbatim}
The above two lines create vectors representing the levels of
each variable for each subject. The {\tt rep} command for r2
says to repeat the sequence {\tt 1 2 3 4}, 8 times. The {\tt
rep} command for {\tt r1} says take the sequence {\tt 1 2 3 4},
then repeat the first element 4 times, the second element 4
times, etc. It does this by using a vector as its second
argument. That vector is {\tt rep(4,4)}, which means repeat the
number 4, 4 times. So {\tt rep(4,4)} is equivalent to {\tt c(4 4
4 4)}. The last argument, 2, in the command for {\tt r1} means
that the whole sequence is repeated twice. Notice that {\tt r1}
and {\tt r2} are the codes for one row of the matrix {\tt rsp1}.
\begin{verbatim}
nsub1 <- nrow(rsp1)
subj1 <- as.factor(rep(1:nsub1,32))
\end{verbatim}
{\tt nsub1} is just the number of subjects (5 in the example), the
number of rows in the matrix {\tt rsp1}. The vector {\tt subj1}
is what we will need to assign a subject number to each
observation. It consists of the sequence {\tt 1 2 3 4 5},
repeated 32 times. It corresponds to the columns of {\tt
rsp1}.
\begin{verbatim}
abr1 <- data.frame(ab1=as.vector(rsp1),sub1=subj1,
dcost1=rep(r1,rep(nsub1,32)),abcost1=rep(r2,rep(nsub1,32)))
\end{verbatim}
The {\tt data.frame} function puts together several vectors into
a data.frame, which has rows and columns like a
matrix.\footnote{The {\tt cbind} function does the same thing but
makes a matrix instead of a data frame.} Each vector becomes a
column. The {\tt as.vector} function reads down by columns, that
is, the first column, then the second, and so on. So {\tt ab} is
now a vector in which the first {\tt nsub1} elements are the same
as the first column of {\tt rsp1}, that is, {\tt 1 1 1 1 1}. The
first 15 elements of ab are: {\tt 1 1 1 1 1 1 2 1 2 1 2 3 2 2 1}.
Notice how we can define names within the arguments to the {\tt
data.frame} function. Of course, {\tt sub} now represents the
subject number of each observation. The first 10 elements of
{\tt sub1} are {\tt 1 2 3 4 5 1 2 3 4 5}. The variable {\tt
abcost} now refers to the value of {\tt r2}. Notice that each
of the 32 elements of {\tt r2} is repeated {\tt nsub} times.
Thus the first 15 values of {\tt abcost1} are {\tt 1 1 1 1 1 2 2
2 2 2 3 3 3 3 3}. Here are the first 10 rows of {\tt abr1}:
\begin{verbatim}
ab1 sub1 dcost1 abcost1
1 1 1 1 1
2 1 2 1 1
3 1 3 1 1
4 1 4 1 1
5 1 5 1 1
6 1 1 1 2
7 2 2 1 2
8 1 3 1 2
9 2 4 1 2
10 1 5 1 2
\end{verbatim}
The following line makes a table of the means of {\tt abr1},
according to the values of {\tt dcost1} (rows) and {\tt abcost1}
(columns).
\begin{verbatim}
ctab1 <- tapply(abr1[,1],list(abr1[,3],abr1[,4]),mean)
\end{verbatim}
It uses the function {\tt tapply}, which is like the {\tt apply}
function except that the output is a table. The first argument
is the vector of data to be used. The second argument is a list
supplying the classification in the table. This list has two
columns corresponding to the columns of abr representing the
classification. The third argument is the function to be applied
to each grouping, which in this case is the mean. Here is the
resulting table:
\begin{verbatim}
1 2 3 4
1 2.6 3.0 3.7 3.8
2 3.5 4.4 4.4 5.4
3 4.5 5.2 5.1 5.9
4 5.1 6.1 6.2 6.8
\end{verbatim}
The following line provides a plot corresponding to the table.
\begin{verbatim}
matplot(ctab1,type="l")
\end{verbatim}
Type {\tt l} means lines. Each line plots the four points in a
column of the table. If you want it to go by rows, use {\tt
t(ctab1)} instead of {\tt ctab1}. The function {\tt t()}
transposes rows and columns.
Finally, the following line does a regression of the response on
the two classifiers, actually an analysis of variance.
\begin{verbatim}
summary(aov(ab1 ~ dcost1 + abcost1 + Error(sub1/(dcost1 + abcost1)),data=abr))
\end{verbatim}
The function {\tt aov}, like {\tt lm}, fits a linear model,
because {\tt dcost1} and {\tt abcost1} are numerical variables, not
factors (although {\tt sub1} is a factor). The model is defined
by its first argument (to the left of the comma), where \verb"~"
separates the dependent variable from the predictors. The second
element defines the data frame to be used. The {\tt summary}
function prints a summary of the regression. (The {\tt lm} and
{\tt aov} objects themselves contains other things, such as
residuals, many of which are not automatically printed.) We
explain the {\tt Error} term later, but the point of it is to
make sure that we test against random variation due to subjects,
that is, test ``across subjects.'' Here is some of the output,
which shows significant effects of both predictors:
\begin{verbatim}
Error: sub1
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 52.975 13.244
Error: sub1:dcost1
Df Sum Sq Mean Sq F value Pr(>F)
dcost1 1 164.711 164.711 233.63 0.0001069 ***
Residuals 4 2.820 0.705
---
Error: sub1:abcost1
Df Sum Sq Mean Sq F value Pr(>F)
abcost1 1 46.561 46.561 41.9 0.002935 **
Residuals 4 4.445 1.111
---
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 145 665.93 4.59
\end{verbatim}
\subsection{Other ways to read in data}
\paragraph{First example.}
Here is another example of creating a matrix with one row per
observation.
\begin{verbatim}
symp1 <- read.table("symp1.data",header=T)
sy1 <- as.matrix(symp1[,c(1:17)])
\end{verbatim}
The first 17 columns of {\tt symp1} are of interest. The file
{\tt symp1.data} contains the names of the variables in its first
line. The {\tt header=T} (an abbreviation for {\tt header=TRUE})
makes sure that the names are used; otherwise the variables will
be names {\tt V1, V2,} etc.
\begin{verbatim}
gr1 <- factor(symp1$group1)
\end{verbatim}%$
The variable {\tt group1}, which is in the original data, is a
factor that is unordered.
The next four lines create the new matrix, defining identifiers
for subjects and items in a questionnaire.
\begin{verbatim}
syv1 <- as.vector(sy1)
subj1 <- factor(rep(1:nrow(sy1),ncol(sy1)))
item <- factor(rep(1:ncol(sy1),rep(nrow(sy1),ncol(sy1))))
grp <- rep(gr1,ncol(sy1))
cgrp <- ((grp==2) | (grp==3))+0
\end{verbatim}
The variable {\tt cgrp} is a code for being in {\tt grp} 2 or 3.
The reason for adding 0 is to make the logical vector of {\tt T}
and {\tt F} into a numeric vector of 1 and 0.
The following three lines create a table from the new matrix,
plot the results, and report the results of an analysis of
variance.
\begin{verbatim}
sytab <- tapply(syv,list(item,grp),mean)
matplot(sytab,type="l")
svlm <- aov(syv ~ item + grp + item*grp)
\end{verbatim}
\paragraph{Second example.}
In the next example, the data file has labels. We want to refer
to the labels as if they were variables we had defined, so we use
the {\tt attach} function.
\begin{verbatim}
t9 <- read.table("tax9.data",header=T)
attach(t9)
\end{verbatim}
\paragraph{Third example.}
In the next example, the data file has no labels, so we can read
it with {\tt scan}. The {\tt scan} function just reads in the
numbers and makes them into a vector, that is, a single column of
numbers.
\begin{verbatim}
abh1 <- matrix(scan("abh1.data"),,224,byrow=T))
\end{verbatim}
We then apply the {\tt matrix} command to make it into a matrix.
(There are many other ways to do this.) We know that the matrix
should have 224 columns, the number of variables, so we should
specify the number of columns. If you say help(matrix) you will
see that the matrix command requires several arguments, separated
by commas. The first is the vector that is to be made into a
matrix, which in this case is \verb'scan("abh1.data")'. We could have
given this vector a name, and then used its name, but there is no
point. The second and third arguments are the number of rows and
the number of columns. We can leave the number of rows blank.
(That way, if we add or delete subjects, we don't need to change
anything.) The number of columns is 224. By default, the matrix
command fills the matrix by columns, so we need to say {\tt
byrow=TRUE} or {\tt byrow=T} to get it to fill by rows, which
is what we want. (Otherwise, we could just leave that field
blank.)
We can refer to elements of abh1 by {\tt abh1[row,column]}. For
example, {\tt abh[1,2]} is the sex of the first subject. We can
leave one part blank and get all of it, e.g., {\tt abh1[,2]} is a
vector (column of numbers) representing the sex of all the
subjects.
\subsection{Other ways to transform variables}
\subsubsection{Contrasts}
Suppose you have a matrix {\tt t1} with 4 columns. Each row is a
subject. You want to contrast the mean of columns 1 and 3 with
the mean of columns 2 and 4. A t-test would be fine.
(Otherwise, this is the equivalent of the {\tt cmatrix} command
in Systat.) Here are three ways to do it. The first way
calculates the mean of the columns 1 and 3 and subtracts the mean
of columns 2 and 4. The result is a vector. When we apply {\tt
t.test()} to a vector, it tests whether the mean of the values
is different from 0.
\begin{verbatim}
t.test(apply(t1[c(1,3),],2,mean)-apply(t1[c(2,4),],2,mean))
\end{verbatim}
The second way multiplies the matrix by a vector representing the
contrast weights, {\tt 1, -1, 1, -1}. Ordinary multiplication of
a matrix by a vector multiplies the rows, but we want the
columns, so we must apply {\tt t()} to transform the matrix, and
then transform it back.
\begin{verbatim}
t.test(t(t(t1)*c(1,-1,1,-1)))
\end{verbatim}
or
\begin{verbatim}
contr1 <- c(1,-1,1,-1)
t.test(t(t(t1)*contr1))
\end{verbatim}
The third way is the most elegant. It uses matrix multiplication
to accomplish the same thing.
\begin{verbatim}
contr1 <- c(1,-1,1,-1)
t.test(t1 %*% contr1)
\end{verbatim}
\subsubsection{Averaging items in a within-subject design}
Suppose we have a matrix {\tt t2}, with 32 columns. Each row is
a subject. The 32 columns represent a 8x4 design. The first 8
columns represent 8 different levels of the first variable, at
the first level of the second variable. The next 8 columns are
the second level of the second variable, etc. Suppose we want a
matrix in which the columns represent the 8 different levels of
the first variable, averaged across the second variable.
\paragraph{First method: loop.}
One way to do it --- inelegantly but effectively --- is with a
loop. First, we set up the resulting matrix. (We can't put
anything in it this way if it doesn't exist yet.)
\begin{verbatim}
m2 <- t2[,c(1:8)]*0
\end{verbatim}
The idea here is just to make sure that the matrix has the right
number of rows, and all 0's. Now here is the loop:
\begin{verbatim}
for (i in 1:8) m2[,i] <- apply(t2[,i+c(8*0:3)],1,mean)
\end{verbatim}
Here, the index {\tt i} is stepped through the columns of {\tt
m2}, filling each one with the mean of four columns of {\tt
t2}. For example, the first column of {\tt m2} is the mean of
columns 1, 9, 17, and 25 of {\tt t2}. This is because the vector
{\tt c(8*0:3)} is 0, 8, 16, 24. The {\tt apply} function uses
{\tt 1} as its second argument, which means to apply the function
{\tt mean} across \textit{rows}.
\paragraph{Second method: matrix multiplication.}
Now here is a more elegant way, but one that requires an
auxiliary matrix, which may use memory if that is a problem.
This time we want the means according to the second variable,
which has four levels, so we want a matrix with four columns. We
will multiply the matrix {\tt t2} by an auxiliary matrix {\tt
c0}.
The matrix {\tt c0} has 32 rows and four columns. The first
column is 1,1,1,1,1,1,1,1 followed by 24 0's. This is the result
of {\tt rep(c(1,0,0,0),rep(8,4))}, which repeats each of the
elements of 1,0,0,0 eight times (since {\tt rep(8,4)} means
8,8,8,8). The second column is 8 0's, 8 1's, and 16 0's.
\begin{verbatim}
c0 <- cbind(rep(c(1,0,0,0),rep(8,4)),rep(c(0,1,0,0),rep(8,4)),
rep(c(0,0,1,0),rep(8,4)),rep(c(0,0,0,1),rep(8,4)))
c2 <- t2 %*% c0
\end{verbatim}
The last line above uses matrix multiplication to create the
matrix {\tt c2}, which has 4 columns and one row per subject.
Note that the order here is important; switching {\tt t2} and
{\tt c0} will not work.
\subsubsection{Selecting cases or variables}
There are several other ways for defining new matrices or data
frames as subsets of other matrices or data frames.
One very useful function is {\tt which()}, which yields the
indices for which its argument is true. For example, the output
of {\tt which(3:10 > 4)} is the vector {\tt 3 4 5 6 7 8}, because
the vector {\tt 3:10} has a length of 8, and the first two places
in it do not meet the criterion that their value is greater than
4. With {\tt which()}, you can use a vector to select rows or
columns from a matrix (or data frame). For example, suppose you
have nine variables in a matrix {\tt m9} and you want to select
three sub-matrices, one consisting of variables 1, 4, 7, another
with 2, 5, 8, and another with 3, 6, 9. Define {\tt mvec} so
that it is the vector {\tt 1 2 3 1 2 3 1 2 3}.
\begin{verbatim}
mvec9 <- rep(1:3,3)
m9a <- m9[,which(mvec9 == 1)]
m9b <- m9[,which(mvec9 == 2)]
m9c <- m9[,which(mvec9 == 3)]
\end{verbatim}
You can use the same method to select subjects by any criterion,
putting the {\tt which()} expression before the comma rather than
after it, so that it indicates rows.
\subsubsection{Recoding and replacing data}
Suppose you have {\tt m1} a matrix of data in which 99 represents
missing data, and you want to replace each 99 with {\tt NA}.
Simply say {\tt m1[m1==99] <- NA}. Note that this will work only
if {\tt m1} is a matrix (or vector), not a data frame (which
could result from a {\tt read.table()} command). You might need
to use the {\tt as.matrix()} function first.
Sometimes you want to recode a variable, e.g., a column in a
matrix. If {\tt q1[,3]} is a 7-point scale and you want to
reverse it, you can say
\begin{verbatim}
q1[,3] <- 8 - q1[,3]
\end{verbatim}
In general, suppose you want to recode the numbers 1,2,3,4,5 so
that they come out as 1,5,3,2,4, respectively. You have a
matrix \texttt{m1}, containing just the numbers 1 through 5. You
can say
\begin{verbatim}
c1 <- c(1,5,3,2,4)
apply(m1,1:2,function(x) c1[x])
\end{verbatim}
In this case \texttt{c1[x]} is just the value at the position
indicated by \texttt{x}.
Suppose that, instead of 1 through 5, you have A through E, so
that you cannot use numerical positions. You want to convert
A,B,C,D,E to 1,5,3,2,4, respectively. You can use two vectors:
\begin{verbatim}
c1 <- c(1,5,3,2,4)
n1 <- c("A","B","C","D","E")
apply(m1,1:2,function(x) c1[which(n1)==x])
\end{verbatim}
Or, alternatively, you can give names to \texttt{c1} instead of
using a second vector:
\begin{verbatim}
c1 <- c(1,5,3,2,4)
names(c1) <- c("A","B","C","D","E")
apply(m1,1:2,function(x) c1[x])
\end{verbatim}
The same general idea will work for arrays, vectors, etc.,
instead of matrices.
Here are some other examples, which may be useful in simple
cases, or as illustrations of various tricks.
In this example, {\tt q2[,c(2,4)]} are two columns that must be recoded
by switching 1 and 2 but leaving responses of 3 or more intact.
To do this, say
\begin{verbatim}
q2[,c(2,4)] <- (q2[,c(2,4)] < 3) * (3 - q2[,c(2,4)]) +
(q2[,c(2,4)] >= 3) * q2[,c(2,4)]
\end{verbatim}
Here the expression {\tt q2[,c(2,4)] < 3} is a two-column matrix
full of {\tt TRUE} and {\tt FALSE}. By putting it in
parentheses, you can multiply it by numbers, and {\tt TRUE} and
{\tt FALSE} are treated as 1 and 0, respectively. Thus, {\tt
(q2[,c(2,4)] < 3) * (3 - q2[,c(2,4)])} switches 1 and 2, for
all entries less than 3. The expression {\tt (q2[,c(2,4)] >= 3)
* q2[,c(2,4)]} replaces all the other values, those greater
than or equal to 3, with themselves.
Here is an example that will switch 1 and 3, 2 and 4,
but leave 5 unchanged, for columns 7 and 9
\begin{verbatim}
q3[,c(7,9)] <- (q3[,c(7,9)]==1)*3 + (q3[,c(7,9)]==2)*4 +
(q3[,c(7,9)]==3)*1 + (q3[,c(7,9)]==4)*2 + (q3[,c(7,9)]==5)*5
\end{verbatim}
Notice that this works because everything on the right of {\tt
<-} is computed on the values in {\tt q3} before any of these
values are replaced.
\subsubsection{Replacing characters with numbers}
Sometimes you have questionnaire data in which the responses are
represented as (for example) ``y'' and ``n'' (for yes and no).
Suppose you want to convert these to numbers so that you can
average them. The following command does this for a matrix {\tt
q1}, whose entries are {\tt y}, {\tt n}, or some other
character for ``unsure.'' It converts {\tt y} to 1 and {\tt n}
to -1, leaving 0 for the ``unsure'' category.
\begin{verbatim}
q1 <- (q1[,]=="y") - (q1[,]=="n")
\end{verbatim}
In essence, this works by creating two new matrices and then
subtracting one from the other, element by element.
\subsection{Using R to compute course grades}
Here is an example that might be useful as well as instructive.
Suppose you have a set of grades including a midterm with two
parts {\tt m1} and {\tt m2}, a final with two parts, and two
assignments. You told the students that you would standardize
the midterm scores, the final scores, and each of the assignment
scores, then compute a weighted sum to determine the grade.
Here, with comments, is an R file that does this. The critical
line is the one that standardizes and computes a weighted sum,
all in one command.
\begin{verbatim}
g1 <- read.csv("grades.csv",header=F) # get the list of scores
a1 <- as.vector(g1[,4])
m1 <- as.vector(g1[,5])
m2 <- as.vector(g1[,6])
a2 <- as.vector(g1[,7])
f1 <- as.vector(g1[,8])
f2 <- as.vector(g1[,9])
a1[a1=="NA"] <- 0 # missing assignment 1 gets a 0
m <- 2*m1+m2 # compute midterm score from the parts
f <- f1+f2
gdf <- data.frame(a1,a2,m,f)
gr <- apply(t(scale(gdf))*c(.10,.10,.30,.50),2,sum)
# The last line standardizes the scores and computes their weighted sum
# The weights are .10, .10, .30, and .50 for a1, a2, m, and f
gcut <- c(-2,-1.7,-1.4,-1.1,-.80,-.62,-.35,-.08,.16,.40,.72,1.1,2)
# The last line defines cutoffs for letter grades.
glabels <- c("f","d","d+","c-","c","c+","b-","b","b+","a-","a","a+")
gletter <- cut(gr,gcut,glabels) # creates a vector of letter grades
grd <- cbind(g1[,1:2],round(gr,digits=4),gletter) # makes a matrix
# gl[,1:2] are students' names
grd[order(gr),] # sorts the matrix in rank order and prints it
round(table(gletter)/.83,1) # prints, with rounding
# the .83 is because there are 83 students, also gets percent
gcum <- as.vector(round(cumsum(table(gletter)/.83),1))
names(gcum) <- glabels
gcum # prints cumulative sum of students with different grades
\end{verbatim}
\section{Graphics}
R can do presentation-quality and publication-quality graphics.
These often require some trial-and-error manipulation of labels,
line styles, axes, fonts, etc., and you should consult the help
pages and the \textit{Introduction} for more details. The
emphasis in this section is on the use of graphics for data
exploration, but we provide some leads into the more advanced
uses.
One trick with graphics is to know how each of the various
graphics commands responds (or fails to respond) to each kind of
data object: {\tt data.frame, matrix,} and {\tt vector}. Often,
you can be surprised.
\subsection{Default behavior of basic commands}
Here is the default behavior for each object for each of some of
the plotting commands, e.g., {\tt plot(x1)} where {\tt x1} is a
vector, matrix, or data frame.
\bigskip\noindent
\begin{tabular}{|l|p{1.75in}|p{1.75in}|p{1.75in}|} \hline
& vector & matrix & data.frame \\ \hline
{\tt plot}
& values as function of position
& 2nd column as function of 1st
& plots of each column as function of others \\ \hline
{\tt boxplot}
& one box for whole vector
& one box for all values in matrix
& one box for each column (variable) \\ \hline
{\tt barplot}
& one bar for each position, height is value
& one bar for each column, summing successive values in
colors
& error \\ \hline
{\tt matplot}
& one labeled point for each position, height is value
& X axis is row, Y is value, label is column
& X axis is row, Y is value, label is column \\ \hline
\end{tabular}\bigskip
The {\tt barplot} of a matrix is an interesting display worth
studying. Each bar is stack of smaller bars in different
colors. Each smaller bar is a single entry in the matrix. The
colors represent the row. Adjacent negative and positive values
are combined. (It is easier to understand this plot if all
values have the same sign.)
\subsection{Other graphics}
To get a bar plot of the column means in a data frame {\tt df1},
you need to say\\
{\tt barplot(height=apply(df1),2,mean))}.
To get a nice parallel coordinate display like that in Systat,
use {\tt matplot} but transform the matrix and use lines instead
of points, that is: {\tt matplot(t(mat1),type="l")}. You can
abbreviate {\tt type} with {\tt t}.
{\tt matplot(v1, m1, type="l")} also plots the columns of the
matrix {\tt m1} on one graph, with {\tt v1} as the horizontal
axis. This is a good way to get plots of two functions on one
graph.
To get scatterplots of the columns of a matrix against each
other, use {\tt pairs(x1)}, where {\tt x1} is a matrix or data
frame. (This is like ``splom'' in Systat, which is the default
graph for correlation matrices.)
Suppose you have a measure {\tt y1} that takes several different
values, and you want to plot histograms of {\tt y1} for different
values of {\tt x1}, next to each other for easy comparison. The
variable {\tt x1} has only two or three values. A good plot is
\verb"stripchart(y1 ~ x1, method='stack')". When {\tt y1} is more
continuous, try \verb"stripchart(y1 ~ x1, method='jitter')".
Here are some other commands in their basic form. There are
several others, and each of these has several variants. You need
to consult the help pages for details.
{\tt plot(v1,v2)} makes a scatterplot of {\tt v2} as a function of
{\tt v1}. If {\tt v1} and {\tt v2} take only a small number of
values, so that the plot has many points plotted on top of each
other, try {\tt plot(jitter(v1),jitter(v2))}.
{\tt hist(x1)} gives a histogram of vector {\tt x1}.
\verb"coplot(y1 ~ x1 | z1)" makes several plots
of {\tt y1} as a function of {\tt x1}, each for a different range
of values of {\tt z1}.
\texttt{interaction.plot(factor1,factor2,v1)} shows how
\texttt{v1} depends on the interaction of the two factors.
Many wonderful graphics functions are available in the Grid and
Lattice packages. Many of these are illustrated and explained in
Venables and Ripley (1999).
\subsection{Saving graphics}
To save a graph as a {\tt png} file, say {\tt png("file1.png")}.
Then run the command to draw the graph, such as {\tt
plot(x1,y1)}. Then say {\tt dev.off()}. You can change the
width and height with arguments to the function. There are many
other formats aside from png, such as pdf, and postscript. See
{\tt help(Devices)}.
There are also some functions for saving graphics already made,
which you can use after the graphic is plotted:
\texttt{dev.copy2eps("file1.eps")} and \texttt{dev2bitmap()}.
\subsection{Multiple figures on one screen}
The {\tt par()} function sets graphics parameters. One type of
parameter specifies the number and layout of multiple figures on
a page or screen. This has two versions, {\tt mfrow} and {\tt
mfcol}. The command {\tt par(mfrow=c(3,2))}, sets the display
for 3 rows and 2 columns, filled one row at a time. The command
f{\tt par(mfcol=c(3,2))} also specifies 3 rows and 2 columns, but
they are filled one column at a time as figures are plotted by
other commands.
Here is an example in which three histograms are printed one
above the other, with the same horizontal and vertical axes and
the same bar widths. The breaks are every 10 units. The {\tt
freq=FALSE} command means that densities are specified rather
than frequencies. The {\tt ylim} commands set the range of the
vertical axis. The {\tt dev.print} line prints the result to a
file. The next three lines print out the histogram as numbers
rather than a plot; this is accomplished with {\tt print=FALSE}.
These are then saved to {\tt hfile1}.
\begin{verbatim}
par(mfrow=c(3,1))
hist(vector1,breaks=10*1:10,freq=FALSE,ylim=c(0,.1))
hist(vector2,breaks=10*1:10,freq=FALSE,ylim=c(0,.1))
hist(vector3,breaks=10*1:10,freq=FALSE,ylim=c(0,.1))
dev.print(png,file="file1.png",width=480,height=640)
h1 <- hist(vector1,breaks=10*1:10,freq=FALSE,ylim=c(0,.1),plot=FALSE)
h2 <- hist(vector2,breaks=10*1:10,freq=FALSE,ylim=c(0,.1),plot=FALSE)
h3 <- hist(vector3,breaks=10*1:10,freq=FALSE,ylim=c(0,.1),plot=FALSE)
sink("hfile1")
h1
h2
h3
sink()
\end{verbatim}
For simple over-plotting, use \texttt{par(new=T)}. Of course,
this will also plot axis labels, etc. To avoid that, you might
say \texttt{par(new=T,ann=F)}. (Apparent undocumented feature:
this setting conveniently disappears after it is used once.) To
plot several graphs of the same type, you can also use
\texttt{points()}, \texttt{lines()}, or \texttt{matplot()}.
\subsection{Other graphics tricks}
When you use \texttt{plot()} with course data (e.g., integers),
it often happens that points fall on top of each other. There
are at least three ways to deal with this. One is to use
\texttt{stripchart()} (see above). Another is to apply
\texttt{jitter()} to one or both of the vectors plotted against
each other, e.g., \texttt{plot(jitter(v1),v2)}. A third is to
use \texttt{sunflowerplot(v1,v2)}, which uses symbols that
indicated how many points fall in the same location.
Use \texttt{identify()} to find the location and index of a point
in a scatterplot made with \texttt{plot()}. Indicate the point
you want by clicking the mouse on it. The function
\texttt{locator()} just gives the coordinates of the point. This
is useful for figuring out where you want to add things to a
plot, such as a legend.
\texttt{text()} uses a vector of strings instead of points in a
plot. If you want a scatterplot with just these name, first make
an empty plot (with \texttt{type="n"}) to get the size of the
plot correct) and then use the text command, e.g.:
\begin{verbatim}
x <- 1:5
plot(x,x^2,type="n")
text(x,x^2,labels=c("one","two","three","four","five"),col=x)
\end{verbatim}
In this case, the \texttt{col=x} argument plots each word in a
different color.
To put a legend on a plot, you can use the ``\texttt{legend=}''
argument of the plotting function, or the \texttt{legend()}
function, e.g.,
\texttt{legend(3,4,legend=c("Self","Trust"),fill=c("gray25","gray75"))}.
This example illustrates the use of gray colors indicated by
number, which is convenient for making graphics for publication.
(For presentation or data exploration, the default colors are
usually excellent.)
Several functions will draw various things on graphs.
\texttt{polygon()} and \texttt{segments()} draw lines. They
differ in the kind of input they want, and the first one closes
the polygon it draws.
\section{Examples of simple graphs in publications}
This section illustrates some practical techniques for making
publication-quality graphs with very basic graphics commands.k
The second author, as editor of the open-access journal
\textit{Judgment and Decision Making}
(\texttt{http://journal.sjdm.org}), has found it necessary to re-draw
some graphs. Usually the originals were made with expensive
proprietary software, most of which is designed for printing on paper
but sometimes is difficult to use for publication graphics, which
usually must be re-sized to fit the journal's format. For this
purpose, the best format is encapsulated PostScript (eps).
But the eps format itself is not enough because of the two types of
graphics formats. \textit{Vector} graphics describe images in terms
of commands, of the form ``draw a black line from point x1,y1 to point
x2,y2.'' Of course, these commands are abbreviated in different ways
for different formats. The details of which points should be black
are left to other software (or to a printer), which is usually
designed to do the best possible job of displaying the element in
question. On the other hand \textit{raster} (or ``bitmap'') images
specify all the points. This works fine if the result is printed on
paper or if the computer software plots the image point by point on
the user's display. If the image must be re-sized, however, the
display program cannot fully recover the original information, and
plots are usually somewhat messy. Common raster formats are tiff,
png, gif, and bmp. Common vector formats are eps, svg (scalable
vector graphics), wmf (Windows meta-file), and emf (extended wmf).
But \textit{all} of these ``vector formats'' can include raster images
within them. Unfortunately, software that claims to produce eps often
does it simply by making a raster image and including it in an eps
file. R is one of the few that makes true vector images correctly.
(Others are Stata and SigmaPlot.)\footnote{Note that the jpeg (or jpg)
format, commonly used for photographs, is a third category. It is
closer to a raster image and is definitely not vector format. But
it re-codes the images into a more compact format by using
interesting mathematical transformations (involving singular value
decomposition), which preserve the useful information found in most
photographs. It is called a ``lossy'' format because it loses some
information. It should not be used for ``line art,'' that is,
graphs, but it is practically necessary for photos. Although most
digital cameras now have a tiff (raster) option, it is rarely used
because it requires much more storage space than jpeg.}
To produce good eps with R, I generally use the following format:
\begin{verbatim}
postscript(file="fig1.eps",width=8,height=8,
horiz=F,onefile=F,pointsize=16,paper="special")
[plotting commands here]
dev.off()
system("bbox fig1.eps")
\end{verbatim}
\noindent
All of these options are described in the help file for
\texttt{postscript()} (in the graphics package), but some comments are
needed. First, \texttt{pointsize} controls the size of the font. A
setting of 14 or 16 is good for a large figure that covers the width
of a page, but usually 18 or 20 is better for a figure that will go in
a single column of a two-column text. Note that the advantage of eps
is that you can resize it, without loss, to fit wherever it
goes.%\footnote{In \LaTeX, the usual way to do that is with
% \verb"\includegraphics[width=\textwidth]{fig1.eps}", or
% \verb"\includegraphics[width=\columnwidth]{fig1.eps}". The latter is
% for inclusion in a single column of a two-column text.}
The \texttt{dev.off()} command is necessary to save the result to
disk.
For all the niceties of R , there is one thing it doesn't do, which is
to make a ``tight bounding box'' around a plot. The difference
between eps and ordinary Postscript is the specification of a bounding
box, which is a description of the box containing the plot. You can
see these commands by looking at the top of almost any eps file, which
is just text. The problem is that R 's default bounding box includes
a blank margin, which you usually do not want. To remove the margin,
I use the following script (which requires Ghostscript and uses
Linux shell commands, which can probably be translated for other
operating systems):
\begin{verbatim}
#!/bin/bash
cat $1 | sed -r -e "s/BoundingBox:[\ ]+[0-9]+[\ ]+[0-9]+[\ ]+[0-9]+
[\ ]+\[0-9]+/`gs -sDEVICE=bbox -dBATCH -dNOPAUSE -q`/" > temp.eps
gs -sDEVICE=bbox -sNOPAUSE -q $1 $showpage -c quit 2> bb.out
sed -e"1 r bb.out" temp.eps > $1
/bin/rm bb.out
/bin/rm temp.eps
\end{verbatim}
\noindent
Note that the first line of this script is folded to make it easier to
read here. It should be unfolded. This script removes the original
bounding box and replaces it with the smallest possible box. The
\texttt{system()} command above simply calls the script.
Each of the following examples is listed according to the URL of the
paper in which it appears. The complete R scripts for these and
other figures are linked from \texttt{http://journal.sjdm.org/RX.html}
\subsection{\tt http://journal.sjdm.org/8827/jdm8827.pdf}
\begin{center}
\includegraphics[width=4in]{8827.fig1.eps}
\end{center}
\begin{verbatim}
postscript(file="fig1.eps",width=8,height=8,
horiz=F,onefile=F,pointsize=16,paper="special")
c1 <- c(683,605)
c2 <- c(648,594)
c3 <- c(619,577)
c4 <- c(520,489)
c5 <- c(525,507)
plot(c1,xlab=expression(bold("Distance from target")),
ylab=expression(bold("Mean RT (milliseconds)")),
ylim=c(475,700),xlim=c(.8,2.2),type="b",xaxt="n")
axis(1,at=c(1,2),labels=c("Near","Far"))
lines(c2,type="b")
lines(c3,type="b")
lines(c4,type="b")
lines(c5,type="b")
text(c(.92,2.08),c1,labels=c1)
text(c(.92,2.08),c2,labels=c2)
text(c(.92,2.08),c3,labels=c3)
text(c(.92,2.08),c4,labels=c4)
text(c(.92,2.08),c5,labels=c5)
text(1.5,mean(c1)+8.6,srt=-23,labels="First quintile: Largest distance-effect slope")
text(1.5,mean(c2)+7,srt=-16,labels="Second quintile")
text(1.5,mean(c3)+7,srt=-12,labels="Third quintile")
text(1.5,mean(c4)-7.4,srt=-9.2,labels="Fourth quintile")
text(1.5,mean(c5)+8.8,srt=-5.4,labels="Fifth quintile: Smallest distance-effect slope")
dev.off()
system("bbox fig1.eps")
\end{verbatim}
This figure illustrates the use of the \texttt{text()} function. Here
I adjusted the slope with \texttt{srt} by trial and error, although
the initial errors got smaller after the first one.
\begin{center}
\includegraphics[width=4in]{8827.fig2.eps}
\end{center}
\begin{verbatim}
postscript(file="fig2.eps",width=8,height=8,
horiz=F,onefile=F,pointsize=16,paper="special")
c1 <- c(0.9,2.0,2.6,3.3,4.0)
c2 <- c(-1.6,-1.3,-1.1,-1.0,-0.8)
c1t <- c("0.9","2.0","2.6","3.3","4.0")
c2t <- c("-1.6","-1.3","-1.1","-1.0","-0.8")
par(oma=c(3,0,0,0))
plot(c1,ylab=expression(bold("Predicted preference")),xlab="",
ylim=c(-2,4),xlim=c(.8,5.2),type="b",xaxt="n")
axis(1,at=c(1:5),padj=1,labels=c("1st quintile\nLargest\ndistance-\neffect\nslope",
"2nd quintile","3rd quintile","4th quintile",
"5th quintile\nSmallest\ndistance-\neffect\nslope"))
mtext(side=1,line=5,text=expression(bold("Distance-effect slope quintile split")))
lines(c2,type="b")
text(1:5,c1-.25,labels=c1t)
text(1:5,c2+.25,labels=c2t)
abline(h=0,lty=2)
text(3,mean(c1)+.65,labels="$10/$15 choice")
text(3,mean(c2)-.5,labels="$100/$110 choice")
dev.off()
system("bbox fig2.eps")
\end{verbatim}
This figure illustrates the use of padj and multi-line labels.
\subsection{\tt http://journal.sjdm.org/8814/jdm8814.pdf}
\begin{center}
\includegraphics[width=4in]{8814.fig0.eps}
\end{center}
\begin{verbatim}
Intent <- array(c(7.32,7.60,7.80,5.28,7.44,8.24,7.96,7.40,8.08,
7.50,6.76,7.48,7.52,7.28,6.48,6.80,7.72,7.48),c(3,3,2))
dimnames(Intent) <- list(Arguments=c(2,4,6),
Background=c("Positive","Negative","None"),
Frame=c("Positive","Negative"))
Intention <- c(Intent)
Arguments <- rep(c(2,4,6),6)
Background <- rep(rep(c("Positive","Negative","None"),c(3,3,3)),2)
Frame <- rep(c("Positive","Negative"),c(9,9))
postscript(file="fig0.eps",width=8,height=8,
horiz=F,onefile=F,pointsize=18,paper="special")
plot(c(2,4,6),Intention[1:3],xlim=c(2,18),ylim=c(5,8.5),pch=19,col="maroon",
xlab="Number of arguments",ylab="Behavioral intention",xaxt="n")
y.lm <- fitted(lm(Intention[1:3] ~ c(2,4,6)))
segments(2, y.lm[1], 6, y.lm[3],col="maroon")
points(c(8,10,12),Intention[4:6],col="maroon",pch=19)
y.lm <- fitted(lm(Intention[4:6] ~ c(8,10,12)))
segments(8, y.lm[1], 12, y.lm[3],col="maroon")
points(c(14,16,18),Intention[7:9],col="maroon",pch=19)
y.lm <- fitted(lm(Intention[7:9] ~ c(14,16,18)))
segments(14, y.lm[1], 18, y.lm[3],col="maroon")
points(c(2,4,6),Intention[10:12],col="blue")
y.lm <- fitted(lm(Intention[10:12] ~ c(2,4,6)))
segments(2, y.lm[1], 6, y.lm[3],col="blue",lty=2)
points(c(8,10,12),Intention[13:15],col="blue")
y.lm <- fitted(lm(Intention[13:15] ~ c(8,10,12)))
segments(8, y.lm[1], 12, y.lm[3],col="blue",lty=2)
points(c(14,16,18),Intention[16:18],col="blue")
y.lm <- fitted(lm(Intention[16:18] ~ c(14,16,18)))
segments(14, y.lm[1], 18, y.lm[3],col="blue",lty=2)
mtext(side=1,line=1,at=c(1,2,3,3.5,4,5,6,6.5,7,8,9)*2,
text=c(2,4,6,"|",2,4,6,"|",2,4,6))
abline(v=7)
abline(v=13)
text(c(4,10,16),5.15,labels=c("Positive\nbackground",
"Negative\nbackground","No\nbackground"))
legend(14,6.42,legend=c("Gain","Loss"),title="Frame:",
col=c("maroon","blue"),pch=c(19,1))
dev.off()
system("bbox fig0.eps")
\end{verbatim}
This is a fairly complicated example, which illustrates several
things. One is the use of \texttt{lm()} to get properties of
best-fitting lines to superimpose on a plot. In simple cases, it is
usually sufficient to say something like \texttt{abline(lm(Y~X))}.
But here the origin is different for each part of the plot, so we use
fitted values and \texttt{segments()} instead of \texttt{abline()}.
Also shown here is the use of \texttt{mtext()} to add text around the
margins of a plot, just as \texttt{text()} adds test internally.
Finally, we use \texttt{legend()} to specify more carefully where the
legend should go.
\subsection{\tt http://journal.sjdm.org/8801/jdm8801.pdf}
\begin{center}
\includegraphics[width=4in]{8801.fig4.eps}
\end{center}
\begin{verbatim}
library(gplots)
c1 <- c(66,69,63,78,40,70,53)
e1 <- c(3,4,4,4,3,4,8)
postscript(file="fig4.eps",width=10.8,height=8,
horiz=F,onefile=F,pointsize=16,paper="special")
barplot2(height=c1,plot.ci=T,ci.u=c1+e1,ci.l=c1-e1,xaxt="n",yaxt="n",
ylab="Prediction accuracy",ylim=c(0,100),width=c(.5,.5),space=1)
axis(1,at=(1:7)-.25,padj=.5,lty=0,
labels=c("Total\n","Mouselab\n","Eye\ntracking","Consistent\ntrials",
"Inconsistent\ntrials","Choice\ntrials","Deferral\ntrials"))
axis(2,at=c(0,25,50,75,100),labels=c("0%","25%","50%","75%","100%"))
text((1:7)-.25,10,labels=paste(c1,"%",sep=""))
text(3.15,72,labels="n.s.")
text(3.15,65.6,labels="}",cex=1.75,lwd=.1)
dev.off()
system("bbox fig4.eps")
\end{verbatim}
For adding confidence intervals, the easiest way is to use the
\texttt{gplots} package, as illustraged here. This plot also
illustrates the use of \texttt{axis()}, and the use of \texttt{cex} to
make a large character, in this case a bracket. The \texttt{lwd}
option is necessary to keep the bracket from being too thick. Trial
and error are needed.
\subsection{\tt http://journal.sjdm.org/8319/jdm8319.pdf}
\begin{center}
\includegraphics[width=4in]{8319.fig1.eps}
\end{center}
\begin{verbatim}
postscript(file="fig1.eps",family="NimbusSan",width=8,height=8,
horiz=F,onefile=F,pointsize=16,paper="special")
plot(0:100,0:100,type="n",axes=F,,xlab="",ylab="")
rect(0,30,20,70,col="#EEEEFF",border=NA)
rect(20,30,45,70,col="#FFFFDD",border=NA)
rect(45,30,95,70,col="#EEEEFF",border=NA)
lines(c(0,0),c(25,75))
lines(c(0,100),c(30,30),col="red",lty=2)
lines(c(0,100),c(70,70),col="red",lty=2)
lines(c(0,100),c(50,50))
segments(x0=c(0,20,45),y0=c(50,60,45),x1=c(20,45,95),y1=c(60,45,70),lwd=2)
points(95,70,cex=4)
mtext("r(t) = value(Left - Right)",side=2)
text(4,65,expression(r(t) == r(t-1) + f(v[target],v[non-target]) + u[t]), pos=4)
text(10,25,"barrier right",pos=4,col="red")
text(10,75,"barrier left",pos=4,col="red")
text(95,77,"choose left")
text(c(0,20,45,85),c(33,33,33,47),labels=c("left","right","left","time"),pos=4)
lines(c(20,20),c(30,70),lty=3)
lines(c(45,45),c(30,70),lty=3)
dev.off()
system("bbox fig1.eps")
\end{verbatim}
This plot illustrates the inclusion of a mathematical expression as
text, as well as the use of \texttt{rect} to make shaded rectangles.
\subsection{\tt http://journal.sjdm.org/8221/jdm8221.pdf}
\begin{center}
\includegraphics[width=4in]{8221.fig1.eps}
\end{center}
\begin{verbatim}
Ch=14
postscript(file="fig1.eps",family="NimbusSan",width=8,height=8,
horiz=F,onefile=F,pointsize=16,paper="special")
plot(c(0,110),c(0,100),type="n",axes=F,xlab="",ylab="")
rect(0,80,20,100,col="gray80")
rect(0+Ch,80-Ch,20+Ch,100-Ch,col="gray80")
rect(0+2*Ch,80-2*Ch,20+2*Ch,100-2*Ch,col="gray80")
rect(0+3*Ch,80-3*Ch,20+3*Ch,100-3*Ch,col="gray80")
rect(0+4*Ch,80-4*Ch,20+4*Ch,100-4*Ch,col="gray80")
text(20,98,pos=4,labels="fixation cue appears")
text(20+Ch,98-Ch,pos=4,labels="saccad targets appear")
text(20+2*Ch,98-2*Ch,pos=4,labels="go cue presented")
text(20+3*Ch,98-3*Ch,pos=4,labels="saccade executed")
text(20+4*Ch,98-4*Ch,pos=4,labels="reward delivered")
points(c(10, 10+Ch,10+Ch-6,10+Ch+6, 10+2*Ch-6,10+2*Ch+6, 10+3*Ch-6),
c(90, 90-Ch,90-Ch-6,90-Ch+6, 90-2*Ch-6,90-2*Ch+6, 90-3*Ch-6),
pch=20)
arrows(0,70,3.5*Ch,70-3.5*Ch,lwd=2)
text(25,80-3*Ch,labels="time")
arrows(10+3*Ch,90-3*Ch,10+3*Ch-5,90-3*Ch-5,length=0.1)
par(new=TRUE)
xspline(87+c(0,2,0,-2),33+c(4,0,-2,0),open=F,shape=c(0,1,1,1))
dev.off()
system("bbox fig1.eps")
\end{verbatim}
This shows how rectangles can overlap. The droplet of sweetened water
in the last rectangle was a bit of a problem because no character
seemed to have the shape of a droplet. Thus, we used
\texttt{xspline()} to draw it piece by piece (with much trial and
error). Another feature is the setting of the constant \texttt{Ch},
which helped get the dimensions right.
\subsection{\tt http://journal.sjdm.org/8120/jdm8120.pdf}
\begin{center}
\includegraphics[width=4in]{8210.fig1.eps}
\end{center}
\begin{verbatim}
P <- c(0.0000616649,0.0012931876,0.0014858932,0.0034575074,0.0095432743,0.0112784208,0.0198140078,
0.0260565422,0.0378525090,0.0476971273,0.0665802025,0.1160787054,0.1561110462,0.1741858728,
0.2592136466,0.3849843314,0.3970805883,0.4387950690,0.5686058809,0.5880746208,0.6367807765,
0.7164637107,0.7548314071,0.8594174096,0.8637551603,0.8852179374,0.8854362373,0.8904200780,
0.9319782385,0.9411071229,0.9474470330,0.9605232158,0.9621474910,0.9716238220,0.9750371388,
0.9800862502,0.9856935080,0.9923052342,0.9993104279,0.9994746329,0.9997647547,0.9999417310,
0.9999506389,0.9999650462,0.9999825779,0.9999967088,0.9999994243,0.9999999681)
ordinate <- sort(P)
n <- length(ordinate)
plotpos <- seq(0.5/n, (n - 0.5)/n, by = 1/n)
postscript(file="fig1.eps",family="NimbusSan",width=8,height=8,
horiz=F,onefile=F,pointsize=16,paper="special") # was 20, changed for sinica paper
plot(ordinate, plotpos, xlab="Expected probability",
ylab="Observed probability")
abline(0,1,lty=3)
grid()
dev.off()
system("bbox fig1.eps")
\end{verbatim}
This example is of substantive interest. It shows the p-values for a
set of t-tests, one for each subject. The abscissa is the percentile
rank of the p-value. If the data were random, the plot would be on
the diagonal. The 5th percentile would have $p=.05$, because 5\% of
the p-values would be significant at the .05 level. As is apparent,
the curve departs from the expectation at both ends, showing that
subjects show significant effects in both directions. This example is
discussed in Baron (in press).
% Baron, J. (in press). Looking at individual subjects in research on
% judgment and decision making (or anything). \textit{Acta Psychologica
% Sinica}. (Special issue on ``Methodological concerns of the
% experimental behavioral researcher: Questions and answers,'' edited by
% S. Li).
\subsection{Boxes and arrows}
This figure illustrates the use of William Revelle's \texttt{psych}
package for making boxes and arrows. (This plot was for an article
but was omitted in the final version.) It also illustrates how to
include Greek letters.
\begin{center}
\includegraphics[width=4in]{06140.fig1.eps}
\end{center}
\begin{verbatim}
library(psych)
xlim=c(0,6)
ylim=c(0,4)
plot(NA,xlim=xlim,ylim=ylim,axes=FALSE,xlab="",ylab="")
r1 <- dia.rect(1,1,"Vividness\nmanipulation",xlim=xlim,ylim=ylim)
r2 <- dia.rect(3,2.3,"Anticipated\nregret",xlim=xlim,ylim=ylim)
r3 <- dia.rect(5,1,"Exchange\ntickets",xlim=xlim,ylim=ylim)
dia.arrow(from=r1$top,to=r2$left)
dia.arrow(from=r1,to=r3)
dia.arrow(from=r2$right,to=r3$top)
text(1.6,1.85,labels=expression(paste(beta == 1.004,", p < .005")),srt=46)
text(4.4,1.85,labels=expression(paste(chi^2 == 6.2,", p < .05")),srt=-46)
text(3,1.2,labels=expression(chi^2 == 3.07))
dev.copy2eps(file="fig1.eps")
dev.off()
system("bbox fig1.eps")
\end{verbatim}
% \subsection{\tt http://journal.sjdm.org//jdm.pdf}
% \begin{center}
% \includegraphics[width=4in]{.eps}
% \end{center}
% \begin{verbatim}
% \end{verbatim}
\section{Statistics}
This section is not a summary of all statistics but, rather, a
set of notes on procedures that are more useful to the kind of
studies that psychology researchers do.
\subsection{Very basic statistics}
Here is a list of some of the commands psychologists use all the
time:
{\tt t.test(a1,b1)} --- Test the difference between the means of
vectors {\tt a1} and {\tt b1}.
{\tt t.test(a1,b1,paired=TRUE)} or {\tt t.test(a1,b1,p=T)}, or, even
simpler, {\tt t.test(b1-a1)} --- Test the difference between means
of vectors {\tt a1} and {\tt b1} when these represent paired
measures, e.g., variables for the same subject. The vectors can
be parts of matrices or data.frames. A good plot to look at is
\begin{verbatim}
plot(a1,b1)
abline(0,1)
\end{verbatim}
This plots {\tt b1} as a function of {\tt a1} and then draws a
diagonal line with an intercept of 0 and a slope of 1. Another
plot is {\tt matplot(t(cbind(a1,b1)),type="l")}, which shows one
line for each pair.
Sometimes you want to do a t-test comparing two groups
represented in the same vector, such as males and females. For
example, you have a vector called {\tt age1} and a vector called
{\tt sex1}, both of the same length. Subject {\tt i1}'s age and
sex are {\tt age1[i1]} and {\tt sex1[i1]}. Then a test to see if
the sexes differ in age is {\tt
t.test(age1[sex1==0],age1[sex1==1])} (or perhaps {\tt
t.test(age1[sex1==0],age1[sex1==1],var.equal=T)} for the
assumption of equal variance). A good plot to do with this sort
of test is\\
{\tt stripchart(age1~sex1,method='jitter')} (or {\tt
stripchart(age1~sex1,method='stack')} if there are only a few
ages represented).
The binomial test (sign test) for asking whether heads are more
likely than tails (for example) uses {\tt prop.test(h1,n1)},
where {\tt h1} is the number of heads and {\tt n1} is the number
of coin tosses. Suppose you have two vectors {\tt a1} and {\tt
b1} of the same length, representing pair of observations on
the same subjects, and you want to find whether {\tt a1} is
higher than {\tt b1} more often than the reverse. Then you can
say {\tt prop.test(sum(a1>b1), sum(a1>b1)+sum(a1+ character), and there
was an error. Then we can type the command
\verb+kappaFor2 <- emacs(kappaFor2)+ to edit its contents with the
emacs editor; or use \verb+vi()+ for the visual editor.
In the above function there
is a \verb+browser()+ command, commented out by a \verb+#+ character.
The \verb+browser()+ command is used to debug a function. Each time R
runs the \verb+kappaFor2+ function, it stops at the point where a
working \verb+browser()+ command was set, and we can debug the
function by examining the variables inside the function.
We run the function \verb+kappaFor2+ and the results show that the
agreement on Item 2 is not reliably greater than 0, with a two-tail
$p$-value of about $0.20$.
\begin{verbatim}
> kappaFor2(r1 = c(0, 0, 0, 0, 1, 1, 1, 1, 1, 1),
r2 = c(0, 1, 1, 1, 1, 1, 1, 1, 1, 1))
kappa S.E. z.stat p.value
0.2857143 0.2213133 1.2909944 0.1967056
\end{verbatim}
\subsection{Generating random data for testing}
Suppose you want to test that the last two formulas yield the
same result, but you don't have any data handy. You can generate
a sample of 10 Likert-type, 5-point scale responses from 100 Ss
as follows:
\begin{verbatim}
v1 <- matrix(sample(c(1, 2, 3, 4, 5), size=1000, replace=T), ncol=10)
\end{verbatim}
\subsection{Within-subject correlations and regressions}
Suppose you have a set of 8 items with 2 measures for each, and
each subject answers both questions about each item. You want to
find out how the answers to the two questions correlate
\textit{within} each subject, across the 8 items. Then you want
to find the average, across subjects, of these within-subject
correlations. The matrices {\tt m1} and {\tt m2} contain the
data for the questions. Each matrix has 8 columns, corresponding
to the 8 items, and one row per subject. The following will do
it:
\begin{verbatim}
nsub1 <- nrow(m1)
cors1 <- rep(NA,nsub1)
for (i in 1:nsub1) cors1[i] <- cor(m1[i,],m2[i,])
\end{verbatim}
The first line finds the number of subjects, {\tt nsub1}. The
second line creates an vector of {\tt nsub1 NA}'s to hold the
correlations, one for each subject. The third line fills this
vector with the correlations using a loop, which uses {\tt i} to
index the subject. Now, if we want to do a {\tt t} test to see
if the correlations are positive, for example, we can say {\tt
t.test(cors1)}.
Similarly, you can store within-subject regressions in a matrix,
as in the following example.
\begin{verbatim}
# first set up a matrix to hold the results, including the intercept
reg1 <- matrix(NA,4,nsub1) # nsub1 is the number of subjects
for (x in 1:ns) reg1[x] <- lm(y1[x,]~x1[x,]+x2[x,]+x3[x,])$coefficients
t.test(reg1[,1]) # is the mean intercept positive?
t.test(reg1[,2]) # is the mean first coefficient positive?
\end{verbatim}%$
This works because {\tt lm()} produces a list, and element {\tt
coefficients} of that list is the coefficients. This element itself may
be decomposed into the intercept, the first coefficient, and so
on.
\subsection{Linear regression and analysis of variance (anova)}
If you want to find whether {\tt y1} depends on {\tt x1} and {\tt
x2}, the basic thing you need is
\begin{verbatim}
lm(y1 ~ x1 + x2)
\end{verbatim}
If these variables are part of a data frame called df1, then you
can sayl \verb"lm(y1 ~ x1 + x2, data=df1)", or you can say {\tt
attach(df1)} before you run the analysis.
Note that {\tt lm()} by itself doesn't do much that is useful.
If you want a summary table, one way to get it is to say
\begin{verbatim}
summary(lm(y ~ x1 + x2))
\end{verbatim}
The coefficients are unstandardized. If you want standardized
coefficients, use
{\tt summary(lm(scale(y) ~ scale(x1) + scale(x2)))}. The {\tt
scale()} function standardizes vectors by default (and it does
many other things, which you can see from {\tt help(scale)}).
Another way to get a summary is with {\tt anova()}. The {\tt
anova()} command is most useful when you want to compare two
models. For example, suppose that you want to ask whether {\tt
x3} and {\tt x4} together account for additional variance after
x1 and x2 are already included in the regression. You cannot
tell this from the summary table that you get from
\begin{verbatim}
summary(lm(y1 ~ x1 + x2 + x3 + x4))
\end{verbatim}
That is because you get a test for each coefficient, but not the
two together. So, you can do the following sequence:
\begin{verbatim}
model1 <- lm(y1 ~ x1 + x2)
model2 <- lm(y1 ~ x1 + x2 + x3 + x4)
anova(model1,model2)
\end{verbatim}
As you might imagine, this is an extremely flexible mechanism,
which allows you to compare two nested models, one with many
predictors not contained in the other. Note that {\tt anova}
reports sums of squares sequentially, building up by adding the
models successively. It is thus different from the usual report
of a multiple regression, {\tt summary(lm(\dots))}. Note also
that you can add and drop variables from a model without retyping
the model, using the functions {\tt add()} and {\tt drop()}.
\subsection{Advanced analysis of variance examples}
We now turn to repeated-measure analysis of variance using the
{\tt aov()} function in R.\footnote{We drop our convention here
of using numbers in made-up variable names, in order to be
consistent with the original names in the examples we cite.}
The {\tt aov()} function is used to produce a Univariate ANOVA
table similar to the one produced by SAS, SPSS, and Systat. The
SAS syntax of an identical analysis is also listed in example 2
for comparison.\footnote{The following examples are inspired by
the examples of Gabriel Baud-Bovy,
\texttt{} contributed to S-News
\texttt{}, entitled ``ANOVAs
and MANOVAs for repeated measures (solved examples),'' dated
2/17/1998.}\footnote{The statistics theory behind the syntax
can be found in the references so detailed explanations are not
provided here. The examples are: \\
1) Hays, Table 13.21.2, p. 518
(1 dependent variable, 2 independent variables: 0 between, 2 within)\\
2) Maxwell and Delaney, p. 497
(1 dependent variable, 2 independent variables: 0 between, 2 within)\\
3) Stevens, Ch. 13.2, p. 442
(1 dependent variable, 1 independent variable: 0 between, 1 within)\\
4) Stevens, Ch. 13.12, p. 468 (1 dependent variable, 3
independent variables: 1 between, 2 within) }
Sometimes the output of {\tt aov()} is different from the output of a
GLM procedure in SAS and SPSS. This may seem confusing for
someone who is more familiar with SAS and SPSS than with R.
We will use example 1 below to show how this can happen and what we
can do to find a suitable alternative solution. We may need to use
a linear mixed-effects models approach using {\tt lme()} instead of
repeated-measures ANOVA using {\tt aov()}.
\subsubsection{Example 1: Mixed Models Approach to Within-Subject Factors (Hays, 1988, Table
13.21.2, p. 518)}\label{EXAMPONE}
We use the data in Hays (1988) to describe how to carry out
repeated-measures ANOVA with R. The example is simplified as follows.
Suppose a psychologist is asked to conduct a study to
help design the control
panel of a machine that delivers medicine by intravenous infusion.
The main purpose of the study is to find the
best shape and color of
the buttons on the control panel to improve efficiency and prevent
potential errors. The psychologist wants to know how quickly users
(physicians) respond to buttons of different colors and shapes.
Suppose that the psychologist
hypothesizes that bright
colors are easier to see than dark colors so the users
respond to them faster. In addition, she thinks that users
can spot circles faster than squares.
These hypotheses may seem trivial, but they involve the same statistical
techniques that can be extended to more complicated experiments.
The psychologists knows that she will not be able to recruit many
physicians to run the test apparatus. Thus she wants to collect as
many test results as possible from each participating physician. Each
participant works with the test apparatus four times, one with round
red buttons, one with square red buttons, one with round gray buttons,
and one with square gray buttons. Here the users only try each
arrangement once, but the psychologist could ask them to repeat the
tests several times in random order to get a more stable response
time. In that case she would have another effect
she may choose to test (repetition, or the effect of learning).
In psychology, an experimental design like this is
often called a ``repeated measures''
design because each respondent gives several responses. It is also
referred to as a within-subject design because
the measurements are made repeatedly within individual subjects. The
variables {\tt shape} and {\tt color} are therefore called
within-subject variables. It is possible to do the experiment between
subjects. Each data point comes from a
different physician. A completely between-subject experiment is also
called a randomized design. However, the experimenter
would need to recruit four times as many physicians, which
is not efficient.
This example has 2 within-subject variables and
no between subject variables:
\begin{itemize}
\item one dependent variable: time required to solve the puzzles
\item one random effect: subject (see Hays for reasons why)
\item 2 within-subject fixed effects: shape (2 levels), color (2
levels)
\end{itemize}
We first enter the reaction time data into a vector {\tt data1}. Then
we transform {\tt data1} into appropriate format for the repeated
analysis of variance using {\tt aov()}.
\begin{verbatim}
data1<-c(
49,47,46,47,48,47,41,46,43,47,46,45,
48,46,47,45,49,44,44,45,42,45,45,40,
49,46,47,45,49,45,41,43,44,46,45,40,
45,43,44,45,48,46,40,45,40,45,47,40) # across subjects then conditions
\end{verbatim}
We can take a look at the data in a layout that is easier to read.
Each subject takes up a row in the data matrix.
\begin{verbatim}
> matrix(data1, ncol= 4, dimnames =
+ list(paste("subj", 1:12), c("Shape1.Color1", "Shape2.Color1",
+ "Shape1.Color2", "Shape2.Color2")))
Shape1.Color1 Shape2.Color1 Shape1.Color2 Shape2.Color2
subj 1 49 48 49 45
subj 2 47 46 46 43
subj 3 46 47 47 44
subj 4 47 45 45 45
subj 5 48 49 49 48
subj 6 47 44 45 46
subj 7 41 44 41 40
subj 8 46 45 43 45
subj 9 43 42 44 40
subj 10 47 45 46 45
subj 11 46 45 45 47
subj 12 45 40 40 40
\end{verbatim}
Next we use the {\tt data.frame()} function to create a data frame {\tt
Hays.df} that is appropriate for the {\tt aov()} function.
\begin{verbatim}
Hays.df <- data.frame(rt = data1,
subj = factor(rep(paste("subj", 1:12, sep=""), 4)),
shape = factor(rep(rep(c("shape1", "shape2"), c(12, 12)), 2)),
color = factor(rep(c("color1", "color2"), c(24, 24))))
\end{verbatim}
The experimenter is interested in knowing if the shape {\tt (shape)}
and the color {\tt (color)} of the buttons affect the reaction time
{\tt (rt)}. The syntax is:
\begin{verbatim}
aov(rt ~ shape * color + Error(subj/(shape * color)), data=Hays.df)
\end{verbatim}
The model formula,
{\tt rt \verb+~+ shape * color + Error(subj/(shape * color))},
can be divided into two parts. The first part,
{\tt rt \verb+~+ shape * color}, states that reaction time
is affected by the shapes and colors of the buttons. The asterisk is
a shorthand for {\tt shape + color +
shape:color}, representing the shape and color main effects and the
shape by color interaction.
The second part, {\tt Error(subj/(shape * color))}, is very important
for getting the appropriate statistical tests.
We will first explain what the syntax means, then we will explain why
we do it this way.
The {\tt Error(subj/(shape * color))} statement is used to separate
the residual sums of squares into several components
(called error strata).
The statement is equivalent to
{\tt Error(subj + subj:shape + subj:color + subj:shape:color)},
meaning that we want to separate the following error terms: one for
subject, one for subject by shape interaction, one for subject by
color interaction, and one for subject by shape by color interaction.
This syntax generates the appropriate tests for the within-subject
variables {\tt shape} and {\tt color}. You get
\begin{verbatim}
> summary(aov(rt ~ shape * color + Error(subj/(shape*color)), data=Hays.df))
Error: subj
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11 226.500 20.591
Error: subj:shape
Df Sum Sq Mean Sq F value Pr(>F)
shape 1 12.0000 12.0000 7.5429 0.01901 *
Residuals 11 17.5000 1.5909
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Error: subj:color
Df Sum Sq Mean Sq F value Pr(>F)
color 1 12.0000 12.0000 13.895 0.003338 **
Residuals 11 9.5000 0.8636
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Error: subj:shape:color
Df Sum Sq Mean Sq F value Pr(>F)
shape:color 1 1.200e-27 1.200e-27 4.327e-28 1
Residuals 11 30.5000 2.7727
\end{verbatim}
Note that the shape effect is tested against the residuals
of the subject by shape interaction,
shown in the {\tt subj:shape} error stratum.
Similarly, the color effect is tested against subject by color stratum.
The last error stratum with {\tt Error:~subj:shape:color}
offers a test of the {\tt shape:color} interaction.
{\bf Random vs. Fixed Effects} \quad In this simplified example we
have one ``random'' subject effect ({\tt subj}) and two ``fixed''
experimental manipulations ({\tt color} and {\tt shape}). In a
repeated-measures design where the within-subject factors are
considered fixed effects and the only random
effect comes from subjects, the
within-subject factors are tested against
their interaction with the random subject effect.
According to Hays (1988), the appropriate {\sf F} tests are the following:
\begin{itemize}
\item[] F(shape in subj) = MS(shape) / MS(shape : subj) = 7.543
\item[] F(color in subj) = MS(color) / MS(color : subj) = 13.895
\end{itemize}
Without going into the details, we can still understand
the distinction between a random and a fixed effect.
Take the fixed {\tt shape} effect as an example, the psychologist only wants
to compare the reaction time differences between round and square
buttons. She is not concerned about generalizing the effect to buttons
of other shapes. In this case the number of possible
shapes is fixed to two---round and square. The reaction time
differences between the two conditions do
not generalize beyond these two shapes. Similarly, the variable
{\tt color} is also considered fixed (again the effect not
generalizable to colors other than red and gray).
When {\tt color} and {\tt shape} are both considered fixed, they
are tested against the {\tt subj:color} and {\tt subj:shape}
mean squares, respectively. The {\tt Error()} function allows you to
do these comparisons.
In this example the only random effect is the {\tt subj}
effect. The 12 subjects reported here belong
to a random sample of
numerous other potential users of the device.
The study would not be very useful without this
generalizability because the
results of the experiments would only apply to these particular 12 test
subjects.
Without the {\tt Error(subj/(shape * color))} formula, you get the
wrong statistical tests. Note that both {\tt color} and {\tt shape}
are tested against a common entry called ``Residuals'', which is the
sum of all the pieces of residuals in the previous output of {\tt
Error(subj/(shape * color))}, with 11 degrees of freedom in each of
the four error strata.
\begin{verbatim}
> summary(aov(rt ~ shape*color, data=hays.df))
Df Sum Sq Mean Sq F value Pr(>F)
shape 1 12.000 12.000 1.8592 0.1797
color 1 12.000 12.000 1.8592 0.1797
shape:color 1 4.399e-29 4.399e-29 6.816e-30 1.0000
Residuals 44 284.000 6.455
\end{verbatim}
{\bf Note about {\tt Error()}} \quad What goes inside the
{\tt Error()} statement, and the order in which
they are arranged, are important in ensuring correct statistical
tests. Suppose you replaced the asterisk inside
{\tt Error(subj/(shape * color))} with a plus sign,
you got {\tt Error(subj/(shape + color))} instead:
\begin{verbatim}
> summary(aov(rt ~ shape * color + Error(subj/(shape + color)),
data=Hays.df))
[identical output as before ... snipped ]
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
shape:color 1 1.185e-27 1.185e-27 4.272e-28 1
Residuals 11 30.5000 2.7727
\end{verbatim}
Note that {\tt Error()} lumps the {\tt shape:color} and {\tt subj:shape:color}
sums of squares into a ``Within'' error
stratum. The ``Residuals'' in the Within stratum is actually the
last piece of sum of square in the previous output
({\tt subj:shape:color}).
This {\tt Error(subj/(shape + color))} syntax gives you the wrong
statistics when you have more than two
within-subject variables. We will return to this point later.
There is a subtle detail in the {\tt Error(subj/(shape * color))}
syntax that is worth noticing. For reasons that will become
clearer in the next section, the forward slash (the {\tt /} in
{\tt subj/(shape * color)}) tells {\tt aov()} that the shape and
color of the buttons are actually {\it nested} within
individual subjects. That is, the changes in response time due to
shape and color should be considered within the subjects.
{\bf Unbalanced (Nonorthogonal) Designs } \quad
We will conclude the
first example by introducing a common complication in
repeated-measures ANOVA that involves uneven cell size.
In the example above, there are four
possible combinations of shape and color, each containing exactly 12
observations. We call this design ``balanced'' because an equal
number of observations is found in every
combination (or a ``cell'') of the design.
A balanced design is also called an ``orthogonal'' design.
In a balanced design, the {\tt aov()}
syntax above produces
exactly the same univariate test statistics as SPSS and SAS.
However, often the cell size is not even, making the design
``unbalanced'' or ``nonorthogonal''. This complication may be due to
participant attrition or a post-hoc grouping factor that was not originally
allotted the same number of participants. For example,
the experimenter may
stratify the design by gender, that is, by recruiting an equal number
of female and male physicians in the study. But later on the
experimental may decide to try analyzing the reaction time by
years of professional experience as the between-subject factor.
When a design is ``unbalanced'',
its {\tt aov()} test statistics may look
very different from those obtained in SPSS and SAS.
The reason is because statistical packages like the
GLM procedure in SPSS adjusts the default Type III Sums of Squares by the
harmonic mean of the unbalanced cell sizes. The adjustment is
discussed in Maxwell and Delaney (1990, pp. 271-297).
SPSS produces the same output as R if the user tells SPSS to
calculate the Type I SS ({\tt SSTYPE(1)}) or Type II SS
({\tt SSTYPE(2)}) instead of the default {\tt SSTYPE(3)}.
As shown in Maxwell and Delaney (1990),
the calculations of SS1 and SS2 do not involve the harmonic mean.
Maxwell and Delaney discuss the pros and cons of
each type of Sums of Squares.
Apparently SPSS and SAS think that the
harmonic mean SS3 is the right analysis.
Afterall, the SS3 is in general what a
behaivoral scientist seeks.
Readers who are interested in the distinctions
between the different types of SS can find a discussion
in Maxwell and Delaney (1990).
The example {\tt aov()}
analysis below can be compared with the results of SPSS using
{\tt SSTYPE(2)}.
\begin{verbatim}
# add one unbalanced between-subject variable
# n=8 in grp 1; 4 in grp 2
Hays.df$grp <- factor(rep(c(1,1,1,1,1,1,1,1,2,2,2,2), 4))
summary(aov(rt ~ grp*color*shape + Error(subj/shape+color), data=Hays.df))
\end{verbatim}
\begin{verbatim}
GLM
Sh1Col1 Sh2Col1 Sh1Col2 Sh2Col2 BY grp
/WSFACTOR = color 2 Polynomial shape 2 Polynomial
/METHOD = SSTYPE(2)
/CRITERIA = ALPHA(.05)
/WSDESIGN = color shape color*shape
/DESIGN = grp .
\end{verbatim}
\subsubsection{Example 2: Maxwell and Delaney, p. 497}
It is the same design as in the previous example, with two within and
a subject effect. We repeat the same R syntax, then we include the
SAS GLM syntax for the same analysis. Here we have:
\begin{itemize}
\item[] one dependent variable: reaction time
\item[] two independent variables: visual stimuli are tilted at 0, 4, and
8 degrees; with noise absent or present. Each subject
responded to 3 tilt x 2 noise = 6 trials.
\end{itemize}
The data are entered slightly differently; their format is like what
you would usually do with SAS, SPSS, and Systat.
\begin{verbatim}
MD497.dat <- matrix(c(
420, 420, 480, 480, 600, 780,
420, 480, 480, 360, 480, 600,
480, 480, 540, 660, 780, 780,
420, 540, 540, 480, 780, 900,
540, 660, 540, 480, 660, 720,
360, 420, 360, 360, 480, 540,
480, 480, 600, 540, 720, 840,
480, 600, 660, 540, 720, 900,
540, 600, 540, 480, 720, 780,
480, 420, 540, 540, 660, 780),
ncol = 6, byrow = T) # byrow=T so the matrix's layout is exactly like this
\end{verbatim}
Next we transform the data matrix into a data frame.
\begin{verbatim}
MD497.df <- data.frame(
rt = as.vector(MD497.dat),
subj = factor(rep(paste("s", 1:10, sep=""), 6)),
deg = factor(rep(rep(c(0,4,8), c(10, 10, 10)), 2)),
noise = factor(rep(c("no.noise", "noise"), c(30, 30))))
\end{verbatim}
Then we test the main effects and the interaction in one aov() model.
The syntax is the same as in the Hays example:
\begin{verbatim}
taov <- aov(rt ~ deg * noise + Error(subj / (deg * noise)), data=MD497.df)
summary(taov)
\end{verbatim}
The results obtained with R can be compared with those obtained using
SAS or SPSS. The following SAS GLM syntax carries out
the same univariate analysis, plus some
multivariate tests. Maxwell and Delaney summarized the conditions
under which one wants to trust the multivariate results more than the
univariate results.
In SAS, each row contains the data from one subject, across 3 degrees
of tilt and two levels of noise. The GLM syntax has a {\tt class}
option where the between-subject factors are listed (if any).
\begin{verbatim}
data rt1;
input deg0NA deg4NA deg8NA deg0NP deg4NP deg8NP;
cards;
420 420 480 480 600 780
420 480 480 360 480 600
480 480 540 660 780 780
420 540 540 480 780 900
540 660 540 480 660 720
360 420 360 360 480 540
480 480 600 540 720 840
480 600 660 540 720 900
540 600 540 480 720 780
480 420 540 540 660 780
;
proc glm data=rt1;
model deg0NA deg4NA deg8NA deg0NP deg4NP deg8NP = ;
repeated noise 2 (0 1), degree 3 (0 4 8) / summary ;
run;
\end{verbatim}
\subsubsection{Example 3: More Than Two Within-Subject Variables}
Earlier we noted that {\tt Error(subj/(shape * color))}, which uses an
asterisk to connect {\tt shape} and {\tt color},
produces detailed breakdown of the variance components. The
{\tt Error(subj/(shape + color))} statement prints out what you
specifically ask for and lumps the remainder into a ``Within''
stratum. If you have more than two within-subject fixed effects,
the latter will produce some undesirable side effects.
The next hypothetical example
\footnote{contributed by Christophe~Pallier} shows
that {\tt aov(a * b * c + Error(subj/(a+b+c)))}
does not give you the appropriate statistical tests for the two-way
interactions among factors {\tt a}, {\tt b}, and {\tt c}.
The problem is because {\tt Error()} lumps everything
other than {\tt Error:~subj:a}, {\tt Error:~subj:b},
and {\tt Error:~subj:c} into a common entry of residuals.
The S Model book by Chambers and Hastie contains some technical
details that explains this behavior.
\begin{verbatim}
subj <- gl(10, 32, 320) # 10 subjects, each tested 32 times, total length 320
a <- gl(2, 16, 320) # first 16 trials with a1 then next 16 with a2
b <- gl(2, 8, 320) # first 8 triasl with b1, then next 8 with b2, etc.
c <- gl(2, 4, 320)
x <- rnorm(320)
d1 <- data.frame(subj, a, b, c, x)
d2 <- aggregate(x, list(a = a, b = b, c = c, subj = subj), mean)
summary(a1 <- aov(x ~ a * b * c + Error(subj/(a*b*c)), d2))
summary(a2 <- aov(x ~ a * b * c + Error(subj/(a+b+c)), d2))
summary(a3 <- aov(x ~ a * b * c + Error(subj/(a*b*c)), d1))
\end{verbatim}
\subsubsection{Example 4: Stevens, 13.2, p.442; a simpler design with only one
within variable}
\begin{verbatim}
data <- c(
30,14,24,38,26,
28,18,20,34,28,
16,10,18,20,14,
34,22,30,44,30)
Stv.df <- data.frame(rt=data,
subj = factor(rep(paste("subj", 1:5, sep=""), 4)),
drug = factor(rep(paste("drug", 1:4, sep=""), c(5,5,5,5))))
\end{verbatim}
We only have one within-subject variable ({\tt drug}) so that the
syntax is simply drug nested within subject.
\begin{verbatim}
summary(aov(rt ~ drug + Error(subj/drug), data = Stv.df))
\end{verbatim}
\subsubsection{Example 5: Stevens pp. 468 -- 474 (one between, two within)}
The original data came from Elashoff (1981).\footnote{``Data for the panel
session in software for repeated measures analysis of variance.''
Proceedings of the Statistical Computing Section of the American
Statistical Association.} It is a test of drug treatment effect by
one between-subject factor: {\tt group} (two groups of 8 subjects
each) and two within-subject factors: {\tt drug} (2 drugs) and {\tt
dose} (3 doses).
\begin{verbatim}
Ela.mat <-matrix(c(
19,22,28,16,26,22,
11,19,30,12,18,28,
20,24,24,24,22,29,
21,25,25,15,10,26,
18,24,29,19,26,28,
17,23,28,15,23,22,
20,23,23,26,21,28,
14,20,29,25,29,29,
16,20,24,30,34,36,
26,26,26,24,30,32,
22,27,23,33,36,45,
16,18,29,27,26,34,
19,21,20,22,22,21,
20,25,25,29,29,33,
21,22,23,27,26,35,
17,20,22,23,26,28), nrow = 16, byrow = T)
\end{verbatim}
We first put them in a multivariate format, using the {\tt
cbind.data.frame()} function.
\begin{verbatim}
Ela.mul <- cbind.data.frame(subj=1:16, gp=factor(rep(1:2,rep(8,2))), Ela.mat)
dimnames(Ela.mul)[[2]] <-
c("subj","gp","d11","d12","d13","d21","d22","d23") # d12 = drug 1, dose 2
\end{verbatim}
Here is the command for transferring it to the univariate format.
\begin{verbatim}
Ela.uni <- data.frame(effect = as.vector(Ela.mat),
subj = factor(paste("s", rep(1:16, 6), sep="")),
gp = factor(paste("gp", rep(rep(c(1, 2), c(8,8)), 6), sep="")),
drug = factor(paste("dr", rep(c(1, 2), c(48, 48)), sep="")),
dose=factor(paste("do", rep(rep(c(1,2,3), rep(16, 3)), 2), sep="")),
row.names = NULL)
\end{verbatim}
The next command performs the correct tests of all effects.
We use {\tt Error(subj/(dose*drug))} to test {\tt gp}, {\tt drug}, {\tt
dose} and their interactions. Worth noting is that R knows that the {\tt gp}
effect goes with the subject error stratum, although we did not
mention {\tt gp} in the {\tt Error()} statement.
\begin{verbatim}
summary(aov(effect ~ gp * drug * dose + Error(subj/(dose*drug)), data=Ela.uni))
\end{verbatim}
% \subsection{Linear Mixed Effects Models}
% Linear mixed effects models using the {\tt nlme()} library is another
% way to analyze repeated-measures ANOVA. In section~\ref{EXAMPONE},
% we note the complication of unbalanced designs that cause
% {\tt aov()} to produce unexpected results. The
% {\tt lme()} function in the {\tt nlme()} library is more suitable for
% unbalanced designs. Using the {\tt Hays.df} data above, the
% {\tt lme()} syntax is listed below
% (note the ``helmert'' contrast instead of the
% default ``treatment'' contrast). Note that the
% p-values are similar to
% those obtained from R {\tt aov()} but are closer to those obtained
% in SPSS GLM.
% \begin{verbatim}
% > options(contrasts = c("contr.helmert", "contr.poly"))
% > lme1 <- lme(rt ~ grp*color*shape, random = ~ 1 | subj, data=Hays.df)
% > summary(lme1)
% Linear mixed-effects model fit by REML
% Data: Hays.df
% AIC BIC logLik
% 210.0505 226.9392 -95.02523
% Random effects:
% Formula: ~1 | subj
% (Intercept) Residual
% StdDev: 2.071131 1.319722
% Fixed effects: rt ~ grp * color * shape
% Value Std.Error DF t-value p-value
% (Intercept) 44.68750 0.6655590 30 67.14281 0.0000
% grp1 -0.93750 0.6655590 10 -1.40859 0.1893
% color1 -0.46875 0.2020404 30 -2.32008 0.0273
% shape1 -0.56250 0.2020404 30 -2.78410 0.0092
% grp1:color1 0.09375 0.2020404 30 0.46402 0.6460
% grp1:shape1 -0.18750 0.2020404 30 -0.92803 0.3608
% color1:shape1 0.09375 0.2020404 30 0.46402 0.6460
% grp1:color1:shape1 0.28125 0.2020404 30 1.39205 0.1741
% Correlation:
% (Intr) grp1 color1 shape1 grp1:c1 grp1:s1 clr1:1
% grp1 0.333
% color1 0.000 0.000
% shape1 0.000 0.000 0.000
% grp1:color1 0.000 0.000 0.333 0.000
% grp1:shape1 0.000 0.000 0.000 0.333 0.000
% color1:shape1 0.000 0.000 0.000 0.000 0.000 0.000
% grp1:color1:shape1 0.000 0.000 0.000 0.000 0.000 0.000 0.333
% Standardized Within-Group Residuals:
% Min Q1 Med Q3 Max
% -1.4294819 -0.6443506 0.1128209 0.5595731 1.6551236
% Number of Observations: 48
% Number of Groups: 12
% \end{verbatim}
% Here is another example that uses {\tt lme()} to
% analyze an ANOVA with a
% nested random effect (p. 439, Maxwell \& Delaney, 1990).
% A director of a psychiatry clinic wants to know if there is a gender
% difference in how psychologist trainees make intake assessments. The
% director randomly selects six psychologists
% trainees, 3 males and 3 females. Each trainee administers intake
% interviews to n = 4 psychiatric patients, and rates the patients'
% general severity (higher rating is more severe). The effects to be
% tested are {\tt trainee} (random) and {\tt sex} (fixed) effects.
% Here {\tt trainee} is nested within {\tt sex}.
% \begin{verbatim}
% > library(nlme)
% > rating <- c(49, 40, 31, 40, 42, 48, 52, 58, 42, 46, 50, 54, 54,
% 60, 64, 70, 44, 54, 54, 64, 57, 62, 66, 71)
% > sex <- factor(c(rep("M", 12), rep("F", 12)))
% > trainee <- factor(paste("tr", rep(1:6, each = 4), sep = ""))
% > nested1.df <- groupedData(rating ~ sex | sex/trainee)
% >
% > lme1 <- lme(rating ~ sex, random = ~ 1 | trainee)
% > summary(lme1)
% Linear mixed-effects model fit by REML
% Data: NULL
% AIC BIC logLik
% 163.0179 167.3821 -77.50895
% Random effects:
% Formula: ~1 | trainee
% (Intercept) Residual
% StdDev: 4.075667 6.749487
% Fixed effects: rating ~ sex
% Value Std.Error DF t-value p-value
% (Intercept) 60 3.055048 18 19.639627 0.0000
% sexM -14 4.320490 4 -3.240373 0.0317
% Correlation:
% (Intr)
% sexM -0.707
% Standardized Within-Group Residuals:
% Min Q1 Med Q3 Max
% -1.8431744 -0.4632926 -0.1155287 0.6459929 1.4263288
% Number of Observations: 24
% Number of Groups: 6
% \end{verbatim}
% Women psychologists trainees gave significantly higher ratings to
% patients than men. Because this is a balanced design, the square of
% the t-statistic is equal to the {\sf F} = 10.5 statistic derived
% by Maxwell \& Delaney (1990, p. 440). The test for the {\tt trainee}
% effect, however, is not directly available in {\tt lme()}. But it can
% be derived from the {\tt Random effects:} output. The {\sf F}
% ratio for the {\tt trainee} effect is
% MS(sex/trainee)/MS(within). According to Maxwell \& Delaney (p. 437,
% Table 10.4), MS(sex/trainee) = Error(within) + n * Error(trainee)
% \verb@6.75^2 + 4 * 4.076^2 = 112@; and MS(within) is just the {\tt Residual}
% squared = \verb+6.75^2 = 45.56+. Thus {\sf F} = 2.45, which is not
% statistically significant.
% The same analysis can be carried out with the {\tt lmer()} function in
% the {\tt lme4} library.
% \begin{verbatim}
% > library(lme4) # generalized linear mixed-effects models by Bates
% Loading required package: Matrix
% Loading required package: lattice
% > lmer(rating ~ sex + (1 | trainee), data = nested1.df)
% Linear mixed-effects model fit by REML
% Formula: rating ~ sex + (1 | trainee)
% Data: nested1.df
% AIC BIC logLik MLdeviance REMLdeviance
% 161.0 164.6 -77.51 163.1 155.0
% Random effects:
% Groups Name Variance Std.Dev.
% trainee (Intercept) 16.611 4.0756
% Residual 45.556 6.7495
% number of obs: 24, groups: trainee, 6
% Fixed effects:
% Estimate Std. Error t value
% (Intercept) 60.000 3.055 19.64
% sexM -14.000 4.320 -3.24
% Correlation of Fixed Effects:
% (Intr)
% sexM -0.707
% \end{verbatim}
% One of the coauthors of the {\tt nlme} library, Douglas Bates, gave
% reasons why {\tt lmer()} and {\tt lme()} do not routinely print out F
% tests for random effects
% (https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html, last
% accessed December, 2006). They can nevertheless be derived if
% necessary (and appropriate).
% This kind of nested design is common. Patients can be treated by
% different therapists (random) who are nested within different
% psychotherapy methods (fixed or random, depending on the research
% question). Therapists can be nested within sites (almost always
% random). {\tt aov()} and {\tt Error()} do not always generate the
% most appropriate statistical tests.
% Often these random effects are treated as fixed. This problem was
% first identified by Clark (1973) (if not the first, one of the first).
% In his ``language-as-fixed-effect fallacy'' paper, Clark showed that
% many studies in linguistics had a mistake in treating the effect
% associated with words as a fixed effect. When the stimuli you present
% to your subjects are a random sample of numerous other possible
% stimuli, the word effect should be treated as a random effect. Many
% statistically significant findings disappeared when Clark treated them
% appropriately as random effects.
% Finally, linear mixed effects
% models are the preferred methods to analyze longitudinal data,
% especially when the repeated measurements are made over staggered
% assessments (at uneven time intervals). In health psychology, often
% assessments are made at baseline (e.g., before treatment), then again
% at 3 and 12 months post treatment completion. Assessments made at
% baseline and 3 months post treatment completion are likely to be
% highly correlated as compared to assessments made between 3 and 12
% months.
\subsubsection{Other Useful Functions for ANOVA}
As we discussed earlier, we can use the {\tt tapply()} function to
calculate the means across various conditions. We can think of it as
using one statement to run the {\tt mean()} function 12 times! The
output matrix is very useful for plotting.
\begin{verbatim}
tapply(Ela.uni$effect, IND = list(Ela.uni$gp, Ela.uni$drug, Ela.uni$dose),
FUN = mean)
\end{verbatim}
We can also easily custom design a function {\tt se()} to calculate
the standard error for the means. R does not have a built-in function
for that purpose, but there is really no need because the standard
error is just the square root (R has the {\tt sqrt()} function) of the
variance ({\tt var()}), divided by the number of observations ({\tt
length()}). We can use one line of {\tt
tapply()} to get all standard errors. The {\tt se()} makes it
easy to find the confidence intervals for those means. Later we will
demonstrate how to use the means and standard errors we got from {\tt
tapply()} to plot the data.
\begin{verbatim}
se <- function(x)
{
y <- x[!is.na(x)] # remove the missing values
sqrt(var(as.vector(y))/length(y))
}
\end{verbatim}
In R, we not only can use the built-in functions such as {\tt
aov()} to do the analyses, we can also take advantage of R's flexibility
and do many analyses by hand. The following examples demonstrate that
some of the ANOVA tests we did earlier with the {\tt aov()} function
can also be done manually with contrasts.
\begin{enumerate}
\item We can use the following contrast to test the group effect. On
the left hand side of the {\tt aov()} model, we use matrix
multiplication ({\tt \%*\%}) to apply the contrast ({\tt contr}) to
each person's 6
data points. As a result, each person's 6 data points become one
number that is actually the person's total score summed across all
conditions. The matrix multiplication is the same as doing {\tt 1 *
d11 + 1 * d12 +
1 * d13 + 1 * d21 + 1 * d22 + 1 * d23} for each person.
Then we use the {\tt aov()} function to compare the total scores
across the two
groups. We can verify that in this output the F statistic for the
{\tt gp} marginal effect is exactly the same as the one in the the
previous {\tt aov(... Error())} output, although the sums of squares
are different because the contrast is not scaled to length
1.
\begin{verbatim}
contr <- matrix(c(
1,
1,
1,
1,
1,
1), ncol = 1)
taov <- aov(cbind(d11,d12,d13,d21,d22,d23) %*% contr ~ gp, data = Ela.mul)
summary(taov, intercept = T)
\end{verbatim}
\item The following contrast, when combined with the {\tt aov()}
function, will test the drug main effect and drug:group interaction.
The contrast {\tt c(1, 1, 1, -1, -1, -1)} applies positive 1's to
columns 1:3 and negative 1's to columns 4:6. Columns 1:3 contain the
data for drug 1 and 4:6 for drug 2, respectively. So the contrast and
the matrix multiplication generates a difference score between drugs 1
and 2. When we use {\tt aov()} to compare this difference among two
groups, we actually test the {\tt drug:gp} interaction.
\begin{verbatim}
contr <- matrix(c(
1,
1,
1,
-1,
-1,
-1), ncol = 1)
tmp<-aov(cbind(d11,d12,d13,d21,d22,d23) %*% contr ~ gp, Ela.mul)
summary(tmp,intercept= T)
\end{verbatim}
\item The next contrast, when combined with the {\tt manova()}
function in R-1.2.0 or later, tests the dose main effect and the
dose:group interaction. The first contrast {\tt c(1, 0, -1, 1, 0, -1)}
tests if the difference between dose 1 and dose 3 are statistically
significant across
groups; and the second contrast {\tt c(0, 1, -1, 0, 1, -1)} tests
the difference between dose 2 and dose 3 across two groups. When
tested simultaneously with {\tt manova()}, we get
\begin{verbatim}
contr <- matrix(c(
1, 0,
0, 1,
-1,-1,
1, 0,
0, 1,
-1,-1), nrow = 6, byrow = T)
\end{verbatim}
\begin{verbatim}
tmp <- manova(cbind(d11,d12,d13,d21,d22,d23) %*% contr ~ gp, Ela.mul)
summary(tmp, test="Wilks", intercept = T)
\end{verbatim}
\item Another {\tt manova()} contrast, which tests {\tt drug:dose}
interaction and {\tt drug:dose:group} interaction.
\begin{verbatim}
contr <- matrix(c(
1,-1,
0, 2,
-1,-1,
-1, 1,
0,-2,
1, 1), nrow = 6, byrow = T)
tmp<-manova(cbind(d11,d12,d13,d21,d22,d23) %*% contr ~ gp, Ela.mul)
summary(tmp, test="Wilks", intercept = T)
\end{verbatim}
\end{enumerate}
\subsubsection{Graphics with error bars}
Next we will demonstrate how to use R's powerful graphics functions to
add error bars to a plot. The example uses the Elashoff example
discussed earlier. In this example we will briefly show how visual
representations compliment the statistical tests. We use R's {\tt
jpg()} graphics driver to generate a graph that can be viewed by a web
browser. The command syntax may appear intimidating for beginners,
but it is worth the increased efficiency in the long run.
You can also use the {\tt postscript()} graphics driver to generate
presentation-quality plots. PostScript files can be transformed into
PDF format so that nowadays the graphs generated by R can be viewed
and printed by virtually anyone.\footnote{One converter is {\tt
ps2pdf}, or try GhostScript.}
Typically the graphs are first generated interactively with drivers
like {\tt X11()}, then the commands are saved and edited into a
script file. A command syntax script eliminates the need to save
bulky graphic files.
First we start the graphics driver {\tt jpg()} and name the
file where the graph(s) will be saved.
\begin{verbatim}
attach(Ela.uni)
jpeg(file = "ElasBar.jpg")
\end{verbatim}
Then we find the means, the standard errors, and the 95\% confidence
bounds of the means.
\begin{verbatim}
tmean <- tapply(effect, IND = list(gp, drug, dose), mean)
tse <- tapply(effect, IND = list(gp, drug, dose), se)
tbarHeight <- matrix(tmean, ncol=3)
dimnames(tbarHeight) <- list(c("gp1dr1","gp2dr1","gp1dr2","gp2dr2"),
c("dose1","dose2" ,"dose3"))
tn <- tapply(effect, IND = list(gp, drug, dose), length)
tu <- tmean + qt(.975, df=tn-1) * tse # upper bound of 95% confidence int.
tl <- tmean + qt(.025, df=tn-1) * tse # lower bound
tcol <- c("blue", "darkblue", "yellow", "orange") # color of the bars
\end{verbatim}
After all the numbers are computed, we start building the barplot.
First we plot the bars without the confidence intervals, axes, labels,
and tick marks. Note that the {\tt barplot()} function returns the
x-axis values at where the center of the bars are plotted. Later we
will use the values in {\tt tbars} to add additional pieces.
\begin{verbatim}
tbars <- barplot(height=tbarHeight, beside=T, space=c(0, 0.5),
ylab="", xlab="", axes=F, names.arg=NULL, ylim=c(-15, 40),
col=tcol)
\end{verbatim}
Then we add the 95\% confidence intervals of the means to the bars.
\begin{verbatim}
segments(x0=tbars, x1=tbars, y0=tl, y1=tu)
segments(x0=tbars-.1, x1=tbars+0.1, y0=tl, y1=tl) # lower whiskers
segments(x0=tbars-.1, x1=tbars+0.1, y0=tu, y1=tu) # upper whiskers
\end{verbatim}
The axes labels are added.
\begin{verbatim}
axis(2, at = seq(0, 40, by=10), labels = rep("", 5), las=1)
tx <- apply(tbars, 2, mean) # center positions for 3 clusters of bars
\end{verbatim}
We plot the horizontal axis manually so that we can ask R to put
things at exactly where we want them.
\begin{verbatim}
segments(x0=0, x1=max(tbars)+1.0, y0=0, y1=0, lty=1, lwd = 2)
text(c("Dose 1", "Dose 2", "Dose 3"), x = tx, y = -1.5, cex =1.5)
mtext(text=seq(0,40,by=10), side = 2, at = seq(0,40,by=10),
line = 1.5, cex =1.5, las=1)
mtext(text="Drug Effectiveness", side = 2, line = 2.5, at = 20, cex =1.8)
\end{verbatim}
Finally we want to plot the legend of the graph manually. R also has
a {\tt legend()} function, although less flexible.
\begin{verbatim}
tx1 <- c(0, 1, 1, 0)
ty1 <- c(-15, -15, -13, -13)
polygon(x=tx1, y=ty1, col=tcol[1])
polygon(x=tx1, y=ty1 + 2.5, col=tcol[2]) # 2nd, moved 2.5 points up
polygon(x=tx1, y=ty1 + 5, col=tcol[3]) # 3rd
polygon(x=tx1, y=ty1 + 7.5, col=tcol[4]) # last
\end{verbatim}
Finally, we add the legend labels for the filled rectangles.
\begin{verbatim}
text(x = 2.0, y = -14, labels="group 1, drug 1", cex = 1.5, adj = 0)
text(x = 2.0, y = -11.5, labels="group 2, drug 1", cex = 1.5, adj = 0)
text(x = 2.0, y = -9, labels="group 1, drug 2", cex = 1.5, adj = 0)
text(x = 2.0, y = -6.5, labels="group 2, drug 2", cex = 1.5, adj = 0)
\end{verbatim}
The greatest effectiveness is attained by subjects in group 2 when
drug 2 is given. This suggests a group by drug interaction, which is
confirmed by the \verb+aov()+ results outlined earlier. It also
indicates an increasing effectiveness from dose 1 to 3, which is
also confirmed by the statistics.
\subsubsection{Another way to do error bars using plotCI()}
The \texttt{gplots} library has a function \texttt{plotCI()},
which does confidence intervals for plots. Here is an example,
which is Figure 1 in
\texttt{http://journal.sjdm.org/06137/jdm06137.htm}. Note the use of a
small horizontal offset ($-.01$) so that the error bars do not
overlap. The font ``NimbusSan'' is supposed to fit well with
Times Roman.
\begin{verbatim}
library("gplots")
m.pg <- c(-2.6400, 3.6000, 6.0000, 3.6800, 5.4400)
se.pg <- c(1.71938, 1.86548, 1.74738, 1.94484, 1.83492)
m.pl <- c(-4.9600, -3.7600, -2.3200, -.1600, 6.5600)
se.pl <- c(1.47024, 1.72170, 1.79139, 1.36587, 1.56852)
postscript(file="fig1.eps",family="NimbusSan",
width=8,height=8,horiz=F,pointsize=18,paper="special")
plotCI(y=c(m.pg,m.pl),x=c(c(1:5)-.01,c(1:5)+.01),uiw=c(se.pg,se.pl),
ylim=c(-6,8),ylab="Net IGT score",xlab="Block",lty=rep(c(1,2),c(5,5)))
lines(y=m.pg,x=c(1:5)-.01,lty=1)
lines(y=m.pl,x=c(1:5)+.01,lty=2)
legend(3.6,-3.7,legend=c("Prior gain","Prior loss"),lty=c(1,2))
dev.off()
\end{verbatim}
\subsection{Use {\tt Error()} for repeated-measure ANOVA}
In this section we give an intuitive explanation to the use of the
{\tt Error()} statement for repeated-measure analysis of variance.
These explanations are different than what are typically covered in
advanced textbooks. The conventional method focuses on
deriving the appropriate error terms for specific statistical
tests. We use an intuitive method, which will show that using {\tt
Error()}
inside an {\tt aov()} function is actually the same as performing
t-tests using contrasts.
The conventional explanation is
computationally and theoretically equivalent to what we are about to
summarize. Detailed theoretical explanations can be found in most
advanced textbooks, including the book by Hoaglin, Mosteller, and
Tukey (1991). Explanations of the technical details can be found in
the book by Chambers and Hastie (1992).
We first review Analysis of Variance using a method commonly seen in
most introductory textbooks. This method uses an ANOVA table to
describes how much of the total variability is accounted for by all
the related variables. An ANOVA table is exactly what {\tt aov()}
does for you. We first apply this method to the {\tt Hays.df} data
described earlier (but repeated here), then we use the ANOVA table to
explain why we must add the {\tt Error()} statement in an {\tt aov()}
command in order to get the appropriate significance tests. Finally
we draw a connection between {\tt Error()} and specific t-tests
tailored for repeated-measure data.
\subsubsection{Basic ANOVA table with {\tt aov()}}
The {\tt aov()} function generates a basic ANOVA table if {\tt
Error()} is not inserted. Applying a simple {\tt aov()} to the {\tt
Hays.df} data, you get an ANOVA table like the following:
\begin{verbatim}
summary(aov(rt ~ subj * color * shape, data = Hays.df))
Df Sum Sq Mean Sq
subj 11 226.500 20.591
color 1 12.000 12.000
shape 1 12.000 12.000
subj:color 11 9.500 0.864
subj:shape 11 17.500 1.591
color:shape 1 1.493e-27 1.493e-27
subj:color:shape 11 30.500 2.773
\end{verbatim}
{\sf R} analyzes how reaction time differs depending on the subjects,
color and the shape of the stimuli. Also, you can have {\sf R} tell
you how they interact with one another. A simple plot of the data may
suggest an interaction between color and shape. A {\tt color:shape}
interaction occurs if, for example, the color yellow is easier to
recognize than red when it comes in a particular shape. The subjects
may recognize yellow squares much faster than any other color and
shape combinations. Therefore the effect of color on reaction time is
not the same for all shapes. We call this an interaction.
The above {\tt aov()} statement divides the total sum of squares in
the reaction time into pieces. By looking at the size of the sums of
squares ({\tt Sum Sq} in the table), you can get a rough idea that
there is a lot of variability among subjects and negligible in the
{\tt color:shape} interaction.
So we are pretty sure that the effect of color does not depend on what
shape it is. The sum of square for {\tt color:shape} is negligible.
Additionally, the {\tt subj} variable has very high variability,
although this is not very interesting because this happens all the
time. We always know for sure that some subjects respond faster
than others.
Obviously we want to know if different colors or shapes make a
difference in the response time. One might naturally think that we
do not need the {\tt subj} variable in the {\tt aov()} statement.
Unfortunately doing so in a repeated design can cause misleading
results:
\begin{verbatim}
summary(aov(rt ~ color * shape, data = Hays.df))
Df Sum Sq Mean Sq F value Pr(>F)
color 1 12.000 12.000 1.8592 0.1797
shape 1 12.000 12.000 1.8592 0.1797
color:shape 1 1.246e-27 1.246e-27 1.931e-28 1.0000
Residuals 44 284.000 6.455
\end{verbatim}
This output can easily deceive you into thinking that there is
nothing statistically significant! This is where {\tt Error()}
is needed to give you the appropriate test statistics.
\subsubsection{Using {\tt Error()} within {\tt aov()}}
It is important to remember that {\tt summary()} generates incorrect
results if you give it the wrong model. Note that in the statement
above the {\tt summary()} function automatically compares each sum
of square with the residual sum of square and prints out the {\sf F}
statistics accordingly. In addition, because the {\tt aov()}
function does not contain the {\tt subj} variable, {\tt aov()} lumps
every sum of squares related to the {\tt subj} variable into this
big {\tt Residuals} sum of squares. You can verify this by adding
up those entries in our basic ANOVA table ($226.5 + 9.5 + 17.5 +
1.49E-27 + 30 = 284$).
{\sf R} does not complain about the above syntax, which assumes
that you want to test each effect against the sum of residual errors
related to the subjects. This leads to incorrect {\sf F} statistics.
The residual error related to the subjects is
not the correct error term for all. Next we will explain how to
find the correct error terms using the {\tt Error()} statement. We
will then use a simple t-test to show you why we want to do that.
\subsubsection{The Appropriate Error Terms}
In a repeated-measure design like that in Hays, the appropriate
error term for the {\tt color} effect is the {\tt
subj:color} sum of squares. Also the error term for the other
within-subject, {\tt shape} effect is the {\tt subj:shape} sum of
squares. The error term for the {\tt color:shape} interaction is
then the {\tt subj:color:shape} sum of squares. A general
discussion can be found in Hoaglin's book. In the next section we
will examine in some detail the test of the {\tt color} effect.
For now we will focus on the appropriate analyses using {\tt
Error()}. We must add an {\tt Error(subj/(shape + color))}
statement within {\tt aov()}. This repeats an earlier analysis.
\begin{verbatim}
summary(aov(rt ~ color * shape + Error(subj/(color + shape)), data = Hays.df))
Error: subj
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11 226.500 20.591
Error: subj:color
Df Sum Sq Mean Sq F value Pr(>F)
color 1 12.0000 12.0000 13.895 0.003338 **
Residuals 11 9.5000 0.8636
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Error: subj:shape
Df Sum Sq Mean Sq F value Pr(>F)
shape 1 12.0000 12.0000 7.5429 0.01901 *
Residuals 11 17.5000 1.5909
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
color:shape 1 1.139e-27 1.139e-27 4.108e-28 1
Residuals 11 30.5000 2.7727
\end{verbatim}
As we mentioned before, the {\tt Error(subj/(color * shape))}
statement is the short hand for dividing all the residual sums of
squares---in this case all subject-related sums of squares---into
three error strata.
The \verb+Error()+ statement says that we want three
error terms separated in the ANOVA table: one for {\tt subj}, {\tt
subj:color}, and {\tt subj:shape}, respectively. The {\tt
summary()} and {\tt aov()} functions are smart enough to do the rest
for you. The effects are arranged according to where they belong.
In the output the {\tt color} effect is tested against the correct
error term {\tt subj:color}, etc. If you add up all the {\tt
Residuals} entries in the table, you will find that it is exactly
$284$, the sum of all subject-related sums of squares.
\subsubsection{Sources of the Appropriate Error Terms}
In this section we use simple examples of t-tests to demonstrate the
need of the appropriate error terms. Rigorous explanations can be
found in Edwards (1985) and Hoaglin, Mosteller, and Tukey (1991).
We will demonstrate that the appropriate error term for an effect in
a repeated ANOVA is exactly identical to the standard error in
a t statistic for testing the same effect.
Let's use the data in Hays (1988), which we show here again as
{\tt hays.mat} (See earlier example for how to read in the
data).
\begin{verbatim}
hays.mat
Shape1.Color1 Shape2.Color1 Shape1.Color2 Shape2.Color2
subj 1 49 48 49 45
subj 2 47 46 46 43
subj 3 46 47 47 44
subj 4 47 45 45 45
subj 5 48 49 49 48
subj 6 47 44 45 46
subj 7 41 44 41 40
subj 8 46 45 43 45
subj 9 43 42 44 40
subj 10 47 45 46 45
subj 11 46 45 45 47
subj 12 45 40 40 40
\end{verbatim}
In a repeated-measure experiment the four measurements of reaction
time are correlated by design because they are from the same subject.
A subject who responds quickly in one condition is likely to respond
quickly in other conditions as well.
To take into consideration these differences, the comparisons of
reaction time should be tested with differences across
conditions. When we take the differences, we use each subject as
his/her own control. So the difference in reaction time has the
subject's baseline speed subtracted out. In the {\tt hays.mat}
data we test the {\tt color} effect by a simple t-test comparing
the differences between the columns of ``Color1'' and ``Color2.''
Using the \verb+t.test()+ function, this is done by
\begin{verbatim}
t.test(x = hays.mat[, 1] + hays.mat[, 2], y = hays.mat[, 3] + hays.mat[, 4],
+ paired = T)
Paired t-test
data: hays.mat[, 1] + hays.mat[, 2] and hays.mat[, 3] + hays.mat[, 4]
t = 3.7276, df = 11, p-value = 0.003338
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.819076 3.180924
sample estimates:
mean of the differences
2
\end{verbatim}
An alternative is to test if a contrast is equal to zero, we
talked about this in earlier sections:
\begin{verbatim}
t.test(hays.mat %*% c(1, 1, -1, -1))
One Sample t-test
data: hays.mat %*% c(1, 1, -1, -1)
t = 3.7276, df = 11, p-value = 0.003338
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
0.819076 3.180924
sample estimates:
mean of x
2
\end{verbatim}
This {\tt c(1, 1, -1, -1)} contrast is identical to the first t-test.
The matrix multiplication (the {\tt \%*\%} operand) takes care of the
arithmetic. It multiplies the first column by a constant
1, add column 2, then subtract from that columns 3 and 4. This tests
the \verb+color+ effect. Note
that the p-value of this t test is the same as the p-values for the
first t test and the earlier {\sf F} test.
It can be proven algebraically that the square of a t-statistic is
identical to the {\sf F} test for the same effect. So this fact can
be used to double check the results. The square of our t-statistic
for {\tt color} is $3.7276^2 = 13.895$, which is identical to the {\sf
F} statistic for {\tt color}.
Now we are ready to draw the connection between a t-statistic for the
contrast and the {\sf F}-statistic in an ANOVA table for
repeated-measure \verb+aov()+. The t statistic
is a ratio between the effect size to be tested and the standard
error of that effect. The larger the ratio, the stronger the effect
size. The formula can be described as follows:
\begin{equation}
t = \frac{\bar x_1 - \bar x_2}{s/\sqrt{n}},
\end{equation}
where the numerator is the observed differences and the denominator
can be interpreted as the expected differences due to chance. If the
actual difference is substantially larger than what you would expect,
then you tend to think that the difference is not due to random
chance.
Similarly, an {\sf F} test contrasts the observed variability with the
expected variability. In a repeated design we must find an
appropriate denominator by adding the the {\tt
Error()} statement inside an {\tt aov()} function.
The next two commands show
that the error sum of squares of the contrast is exactly identical to
the {\tt Residual} sum of squares for the {\tt subj:color} error
stratum.
\begin{verbatim}
tvec <- hays.mat %*% c(1, 1, -1, -1)/2
sum((tvec - mean(tvec))^2)
[1] 9.5
\end{verbatim}
The sum of squares of the contrast is exactly $9.5$, identical to the
residual sum of squares for the correct {\sf F} test. The scaling
factor $1/2$ is critical because it provides correct scaling for the
numbers. By definition a statistical contrast should have a vector
length of 1. This is done by dividing each element of the contrast
vector by 2, turning it to {\tt c(1/2, 1/2, -1/2, -1/2)}. The scaling
does not affect the t-statistics. But it becomes important when
we draw a connection between a t-test and an {\sf F} test.
You get the standard error of the t-statistic if you do the
following:
\begin{verbatim}
sqrt(sum((tvec - mean(tvec))^2 / 11) / 12)
[1] 0.2682717
\end{verbatim}
The first division of 11 is for calculating the variance; then you
divide the variance by the sample size of 12, take the square root,
you have the standard error for the t-test. You can verify it
by running {\tt se(hays.mat \%*\% c(1, 1, -1, -1)/2)}.
\subsubsection{Verify the Calculations Manually}
All the above calculations by \verb+aov()+ can be verified manually.
This section
summarizes some of the technical details. This also gives you a
flavor of how Analysis Of Variance can be done by matrix algebra.
First we re-arrange the raw data into a three-dimensional array.
Each element of the array is one data point, and the three
dimensions are for the subject, the shape, and the color,
respectively.
\begin{verbatim}
hays.A <- array(hays.dat, dim=c(12, 2, 2))
dimnames(hays.A) <- list(paste("subj", 1:12, sep = ""),
+ c("Shape1", "Shape2"), c("Color1", "Color2"))
\end{verbatim}
Because at this point we want to solve for the effect of {\tt
color}, we use the {\tt apply()} function to average the reaction
time over the two shapes.
\begin{verbatim}
Ss.color <- apply(hays.A, c(1, 3), mean) # Ss x color: average across shape
\end{verbatim}
Next we test a t-test contrast for the color effect, which is the same as
\verb"t.test(Ss.color %*% c(1, -1))". Also note that the square of
the t statistic is exactly the same as the {\sf F} test.
\begin{verbatim}
Contr <- c(1, -1)
Ss.color.Contr <- Ss.color %*% Contr
mean(Ss.color.Contr) / (sqrt(var(Ss.color.Contr) / length(Ss.color.Contr)))
[,1]
[1,] 3.727564
\end{verbatim}
The above t-test compares the mean of the contrast against
the standard error of the contrast, which is
{\tt sqrt(var(Ss.color.Contr) / length(Ss.color.Contr))}
Now we can verify that the sum of square of the contrast is exactly the
same as the error term when we use {\tt aov()} with the
{\tt Error(subj:color)} stratum.
\begin{verbatim}
sum((Ss.color.Contr - mean(Ss.color.Contr))^2)
[1] 9.5
\end{verbatim}
\subsection{Sphericity}
Sphericity is an assumption of repeated measure ANOVA. It means that
the variance-covariance structure of the repeated measure ANOVA
follows a certain pattern. In a repeated measure design, several
measurements are made to the same groups of subjects under different
conditions (and/or at different time). Sphericity is, in a nutshell,
that the variances of the differences between the repeated
measurements should be about the same. This is a complicated concept
best explained with examples. Simple examples tend to oversimplify
things, because they often rely on idealized hypothetical data that
rarely occur in real life. Although their usefulness in relieving the
burden of learning overcompensates for the over-simplification. In
what follows the example in section (cite the previous section???) is
used again to explain what sphericity is, how to calculate and test
it, and what alternative statistics can be used if the sphericity
assumption is not met.
\subsubsection{Why is sphericity important?}
Violations of the sphericity assumption lead to biased p-values. The
alpha error of a test may be set at 5\%, but the test may be actually
rejecting the null hypothesis 10\% of the time. This raises doubts
of the conclusions of the repeated measure ANOVA.
Imagine a study about weight gains
of new born babies at the Intensive Care Unit. The weight of the
babies is measured every day for a critical period of three days. On
the average the babies gain 100 grams between days 1 and 2, and they
gain 150 grams between days 2 and 3. Sphericity says that the
variance of the weight gain between days 1 and 2 should be about the
same as the variance of the weight gain between days 2 and 3 (and also
between days 1 and 3). If not, the variance observed between
different time periods are confounded with the correlation of the
measurements made at different time.
Suppose the variance of the first weight gain is 20 and the
variance of the second weight gain is 100, then the measurements made
at times 1 and 2 are likely to be correlated more closely than
measurements made at times 2 and 3. As a result
the variance over the whole 3-day period (what is often called the
variance of the ``time'' effect in ANOVA jargon) fluctuates over time
and is not reliable in describing the overall growth pattern in
critically ill babies.
In repeated measure experiments the same subjects are tested multiple
times under different conditions. It is a good idea to check if the
responses made under some conditions are correlated more closely than
responses made under other conditions.
There is a statistic, The Greenhouse-Geisser epsilon $\epsilon$, which
measures
by how much the
sphericity assumption is violated. Epsilon is then used to
adjust for the
potential bias in the
$F$ statistic.
Epsilon can be 1, which
means that the sphericity assumption is met perfectly.
An epsilon
smaller than 1 means that the sphericity assumption is violated. The
further it deviates from 1, the worse the violation. In real life
epsilon is rarely exactly 1. If it is not much smaller than 1, then we feel
comfortable with the results of repeated measure
ANOVA. Thus the question is how small is too small. We will get to
that below. Additionally, we will talk about two remedies
when sphericity is violated: 1) correct for the p-value, and 2) use
procedures that do not depend on sphericity, such as MANOVA.
\subsubsection{How to estimate the Greenhouse-Geisser epsilon?}
The Greenhouse-Geisser epsilon is derived from the variance-covariance
matrix of the data. The \verb+MD497.dat+ example above involves a
study where study subjects judged stimuli under 3 different angles of
rotation, at $0$, $2$, and $4$ degrees angle from the horizontal. In
this subsection we estimate the Greenhouse-Geisser epsilon associated
with the rotation of the stimuli. The 3 measurements of reaction time
are stored in \verb+x0+, \verb+x2+, and \verb+x4+, respectively.
\begin{verbatim}
x0 <- apply(MD497.dat[, c(1, 4)], 1, mean) # avg across noise present/absent
x2 <- apply(MD497.dat[, c(2, 5)], 1, mean)
x4 <- apply(MD497.dat[, c(3, 6)], 1, mean)
\end{verbatim}
We need to first calculate the variance-covariance matrix of the 3
variables. The \verb+var()+ function calculates the
variance-covariance matrix if the data are arranged in a matrix, like
\verb+S <- var(cbind(x0, x2, x4))+:
\begin{verbatim}
> var(cbind(x0, x2, x4))
x0 x2 x4
x0 4090 3950 4350
x2 3950 6850 6150
x4 4350 6150 8850
\end{verbatim}
The diagonal entries are the variances and the off diagonal entries are
the covariances. From this variance-covariance matrix the $\epsilon$
statistic can be estimated:
$$
\hat \epsilon = \frac{k^2 (\bar s_{ii} - \bar s)^2}
{(k - 1) (\sum\!\sum s_{ij}^2 - 2k \sum \bar s_{i.}^2 + k^2 \bar s^2) },
$$
\noindent where $\bar s_{ii}$ is the mean of the entries on the main
diagonal of
\verb+S+, which can be shown by \verb+mean(diag(S))+ to equal
6596.667;
$\bar s$ is the mean of all entries, $\bar s_{i.}$ is the mean of all
entries in row $i$ of $S$, and $s_{ij}$ is the 9 individual
entries in the variance-covariance matrix.
\begin{verbatim}
S <- var(cbind(x0, x2, x4))
k <- 3; D < - k^2 * ( mean ( diag ( S ) ) - mean ( S ) ) ^2
N1 <- sum(S^2)
N2 <- 2 * k * sum(apply(S, 1, mean)^2)
N3 <- k^2 * mean(S)^2
D / ((k - 1) * (N1 - N2 + N3))
\end{verbatim}
\noindent Which returns \verb+[1] 0.9616365+ for
the value of $\hat \epsilon$. This value
rounds to the 0.9616 value calculated by SAS (the SAS GLM
syntax is provided above).
There are three important values of $\epsilon$.
It can be $1$ when
the sphericity is met perfectly. It can be as low as
$\epsilon = 1/(k - 1)$,
which produces the lower bound of epsilon (the worst case
scenario). The worst case scenario depends on $k$, the number of
levels in the repeated measure factor. In this example $k = 3$.
Each subject is tested under 3 different levels of stimuli
rotation. Epsilon can be as low as $0.50$ when $k = 3$. Note that
the sphericity assumption does not apply when $k = 2$. Another way to
view it is that even the lowest epsilon is 1. Thus a repeated measure
design with only 2 levels does not involve violations of the
sphericity assumption.
Adjustment of the $F$ statistic can be made against either the
estimated epsilon, in this case $\hat \epsilon = 0.962$; or against
the worst case epsilon of $0.50$. It depends on how conservative one
wants to be. If the cost of falsely rejecting the null hypothesis is
high, then one may want to adjust against the worst possible (and very
conservative) epsilon. Both SPSS and SAS use the estimated
value to make the Greenhouse-Geisser adjustment. The
Greenhouse-Geisser adjustment made by SPSS and SAS is different from
the adjustment originally proposed by Greenhouse and Geisser
(1959). Although the adjustment made by SPSS and SAS is
considered more reasonable (Stevens, 1990).
The estimated epsilon is used to adjust
the degrees of freedom associated with the $F$ statistic from the
unadjusted $(k - 1)$ and $(k - 1)(n - 1)$ to $\epsilon (k - 1)$ and
$\epsilon (k - 1)(n - 1)$. Severe violations of the sphericity
assumption (as $\hat \epsilon \rightarrow 0$) may decrease the degrees
of freedom so much that the $F$ statistic is no longer statistically
significant. The p-value associated with the adjusted
$F$ can be obtained from the
\verb+pf()+ function.
>From the previous \verb+aov()+ output we get an $F$ statistic of
40.719 for the variable \verb+deg+. The numerator
degree of freedom is $(k - 1) = 2$ and the denominator degrees of
freedom is $(k - 1)(n - 1) = (3 - 1)(10 - 1) = 18$. These can be
verified with the output of the previous analysis.
Suppose the value of epsilon is assigned to
{\tt epsi <- D / ((k - 1) * (N1 - N2 + N3))}.
We can then use \verb+epsi+ to weigh down
the degrees of freedom.
\begin{verbatim}
> 2 * (1 - pf(40.719, df1=2*epsi, df2=18*epsi))
[1] 6.80353e-07
\end{verbatim}
The $F$ statistic is still statistically
significant below the 0.0001 level. The negligible violation of the
sphericity assumption does not affect the conclusion we
make.
\paragraph{Huynh-Feldt correction.}
The Greenhouse-Geisser epsilon tends to underestimate epsilon when
epsilon is greater than 0.70 (Stevens, 1990). An estimated
$\epsilon = 0.96$ may be actually 1.
Huynh-Feldt correction is less conservative. The Huynh-Feldt epsilon
is calculated from the Greenhouse-Geisser epsilon,
$$
\bar \epsilon = \frac{n (k - 1) \, \hat \epsilon - 2}%
{(k - 1)[\, (n - 1) - (k - 1) \, \hat \epsilon\, ]}.
$$
The Huynh-Feldt epsilon can be easily calculated:
\begin{verbatim}
> epsiHF <- (10 * (k-1) * epsi - 2) / ((k-1) * ((10-1) - (k-1)*epsi))
[1] 1.217564
> 2 * (1 - pf(40.719, df1=2*epsiHF, df2=18*epsiHF))
[1] 2.633106e-08
\end{verbatim}
The Huynh-Feldt epsilon is 1.2176 with an adjusted p-value lower than
0.0001. Again, the Huynh-Feldt correction does not change the
conclusion. In this example, the univariate results are preferred
because an $\epsilon = 0.96$ is very close to 1.0. The
sphericity corrections do not change the conclusions.
MANOVA (and multivariate tests) may be better if the
Greenhouse-Geisser and the Huynh-Feldt corrections do not agree, which
may happen when epsilon drops below $0.70$. When epsilon drops below
$0.40$, both the G-G and H-F corrections may indicate that the
violation of sphericity is affecting the adjusted p-values.
MANOVA is not always appropriate, though. MANOVA usually requires a
larger sample size. Maxwell and Delaney
(1990, p.~602) suggest a rough rule of thumb that the sample size
$n$ should be greater than $k + 10$. In
the present example the sample size $n$ is 10, which is smaller than
$k + 10 = 13$. In such cases the MANOVA procedures are not
recommended. Fortunately we already have strong univariate results.
\subsection{Logistic regression}
Multiple regression is usually not appropriate when we regress a
dichotomous (yes-no) variable on continuous predictors. The
assumptions of normally distributed error are typically violated.
So we usually use logistic regression instead. That is, we
assume that the probability of a ``yes'' is certain function of
a weighted sum of the predictors, the inverse logit. In other
words, if $Y$ is the probability of a ``yes'' for a given set of
predictor values $X_1$, $X_2$, \dots, the model says that
\[ \log \frac{Y}{1-Y} = b_0 + b_1X_1 + b_2X_2 + \dots + \mbox{error} \]
The function $\log \frac{Y}{1-Y}$ is the logit function. This is
the ``link function'' in logistic regression. Other link
functions are possible in R. If we represent the right side of
this equation as $X$, then the inverse function is
\[ Y = \frac{e^X}{1+e^X} \]
In R, when using such transformations as this one, we use {\tt
glm} (the generalized linear model) instead of {\tt lm}. We
specify the ``family'' of the model to get the right
distribution. Here the family is called {\tt binomial}. Suppose
the variable {\tt y} has a value of 0 or 1 for each subject, and
the predictors are {\tt x1, x2}, and {\tt x3}. We can thus say
\begin{verbatim}
summary(glm(y ~ x1 + x2 + x3, family=binomial))
\end{verbatim}
\noindent
to get the basic analysis, including {\tt p} values for each
predictor. Psychologists often like to ask whether the overall
regression is significant before looking at the individual
predictors. Unfortunately, R does not report the overall
significance as part of the {\tt summary} command. To get a test
of overall significance, you must compare two models. One way to
do this is:
\begin{verbatim}
glm1 <- glm(y ~ x1 + x2 + x3, family=binomial)
glm0 <- glm(y ~ 1, family=binomial)
anova(glm0,glm1,test="Chisq")
\end{verbatim}
\subsection{Log-linear models}
Another use of {\tt glm} is log-linear analysis, where the family
is {\tt poisson} rather than {\tt binomial}. Suppose we have a
table called {\tt t1.data} like the following (which you could
generate with the help of \texttt{expand.grid()}). Each row
represents the levels of the variables of interest. The last
column represents the number of subjects with that combination of
levels. The dependent measure is actually expens vs.\ notexpens.
The classification of subjects into these categories depended on
whether the subject chose the expensive treatment or not. The
variable ``cancer'' has three values (cervic, colon, breast)
corresponding to the three scenarios, so R makes two dummy
variables, ``cancercervic'' and ``cancercolon''. The variable
``cost'' has the levels ``expens'' and ``notexp''. The variable
``real'' is ''real'' vs. ``hyp'' (hypothetical).
\begin{verbatim}
cancer cost real count
colon notexp real 37
colon expens real 20
colon notexp hyp 31
colon expens hyp 15
cervic notexp real 27
cervic expens real 28
cervic notexp hyp 52
cervic expens hyp 6
breast notexp real 22
breast expens real 32
breast notexp hyp 25
breast expens hyp 27
\end{verbatim}
The following sequence of commands does one analysis:
\begin{verbatim}
t1 <- read.table("t1.data",header=T)
summary(glm(count ~ cancer + cost + real + cost*real,
family=poisson(), data=t1)
\end{verbatim}
This analysis asks whether ``cost'' and ``real'' interact in
determining ``count,'' that is, whether the response is affected
by ``real.'' See the chapter on Generalized Linear Models in
Venables and Ripley (1999) for more discussion of how this
works.
\subsection{Rasch models}
R can be used to fit a Rasch model in Item Response Theory (IRT).
The Rasch model was proposed in the 1960's by a Danish statistician
Georg Rasch.
The basic Rasch model is used to separate the ability of
test takers and the quality of
the test. In this section we will be going over
how R can be used to estimate IRT models.
Take the following table as an example. The table shows the
results of 35 students taking an 18-item ability test. Each item
is coded as correct (1) or incorrect (0). The data come from the
first example of a popular free software called MiniStep.
MiniStep is the evaluation version of a larger, more complete
software called WinStep. Ministep can be downloaded from
http://www.winstep.org. It works only under the Windows
environment.
At the bottom of the table is a row showing the percentages of
students who answer each question correctly. Questions 1, 2, and 3
are the easiest because all students get it right. Question 18 is the
hardest because none of the students knows the answer to it.
Questions like these are not useful in distinguishing students who
possess high ability and students who possess low ability. Questions
like these are therefore omitted from the calculation.
\begin{verbatim}
name q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12 q13 q14 q15 q16 q17 q18
Richard 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
Tracie 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
Walter 1 1 1 1 1 1 1 1 1 0 0 1 0 0 0 0 0 0
Blaise 1 1 1 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0
Ron 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
William 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
Susan 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 0
Linda 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
Kim 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
Carol 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
Pete 1 1 1 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0
Brenda 1 1 1 1 1 0 1 0 1 1 0 0 0 0 0 0 0 0
Mike 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 0 0 0
Zula 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
Frank 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
Dorothy 1 1 1 1 1 1 1 1 1 0 1 0 0 0 0 0 0 0
Rod 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 0 0 0
Britton 1 1 1 1 1 1 1 1 1 1 0 0 1 0 0 0 0 0
Janet 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0
David 1 1 1 1 1 1 1 1 1 1 0 0 1 0 0 0 0 0
Thomas 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 0 0 0
Betty 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
Bert 1 1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 0 0
Rick 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 1 1 0
Don 1 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0
Barbara 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
Adam 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
Audrey 1 1 1 1 1 1 1 1 1 0 1 0 0 0 0 0 0 0
Anne 1 1 1 1 1 1 0 0 1 1 1 0 0 1 0 0 0 0
Lisa 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0
James 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
Joe 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
Martha 1 1 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0
Elsie 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 0 0 0
Helen 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
----------------------------------------------------------------------------
1 1 1 .91 .88 .86 .88 .77 .86 .69 .34 .17 .20 .09 .03 .03 .03 0.0
\end{verbatim}
We also need to examine if there are students who get all or none of
the items right. If we omit items 1, 2, 3, and 18, and we take a look
at the the 0's and 1's for each student, then we see that Helen fails
on all remaining items. Helen therefore tells us nothing
about which item is harder and which is easier. All items are beyond
Helen's ability.
Students who can get all items
right also do not give us useful information.
But in this example none of the students can answer all
items correctly.
This is how many specific-purpose statistical packages prepare the raw
data for IRT analysis. These packages set aside any items or persons
that provide no useful information for the analysis. They analyze the
IRT model with the remaining data. From the solutions derived from
the remaining data, these packages extrapolate to come up with
estimates for items and persons first set aside.
In R we can use a conditional logit model to perform IRT analyses.
The conditinal logit model is in the survival analysis library
maintained by Thomas~Lumley at University of
Washington. A description of the survival analysis library and the
additional packages it depends on can be
obtained by the R command \verb+library(help=survival)+. The R codes
offered in this section have been tested with version 2.11-4 of the
survival library in R-1.9.0.
The original derivation by Rasch was different but equivalent to a
fixed-effect conditional logit model. Note that the following
procedures in R are suitable only for binary items (e.g., pass/fail on
test items).
The Rasch model states that the log odds of a person $s$ answering an
item $i$ correctly is a function of the person's ability ($\theta_s$)
and the item's difficulty ($\beta_i$). Difficult items are hard to
get right even for people with high ability. The odds of getting an
item right decreases with item difficulty (and thus the minus sign
before $\beta_i$). The notations we use here are from Embretson and
Reise (2000).
$$
{\rm logit}(p_{i,s}) = \log({\Pr(p_{i,s}) \over \Pr(1-p_{i,s})}) =
\theta_s - \beta_i,
$$
where
$$
s = 1, \ldots,
{\rm number\ of\ persons}; \quad i = 1, \ldots, {\rm number\ of\ items}; \quad
{\rm and\ }
\theta_s \approx {\rm Normal}(0, \tau).
$$
The example data, after removing Helen and items 1, 2, 3, and 18,
contain 34 students ($s = 1, 2, \ldots 34$) and
17 items ($i = i, 2, \ldots 17$). Richard is
person 1, Tracie is person 2, ..., and so on. To evaluate Richard's
chance of answering item 3 correctly, we take Richard's abiliy (say,
it is 3 in logit unit), substract it out with item 3's difficulty
(say, it is 2.2 in logit unit), and we get a 0.8 logit. We can convert the
logit scale back to probability by taking \verb=exp(0.8)/(1 + exp(0.8))=
% note in \verb I changed the '+' to '=' or the plus sign disappears
and we get \verb+0.69+. Richard has a 69\% chance of getting item 3
right.
The question is, how do we calculate the $\theta$'s and the $\beta$'s
in the equation?
We need to first load the survival analysis library by typing
\verb+library(survival)+. Then we
need to set aside uninformative items and persons. This is done by
the following commands:
\begin{verbatim}
> library(survival)
> exam1 <- read.csv(file="exam1.dat", sep="", header=T, row.names=NULL)
> exam1$id <- factor(paste("s", 1:dim(exam1)[1], sep=""))
> #
> # percentage of Ss who can correctly answer the items
> #
> round(apply(exam1[, 3:20], 2, mean), 2) # items 1, 2, 3 and 18 are no good
q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12 q13 q14 q15 q16
1.00 1.00 1.00 0.91 0.89 0.86 0.89 0.77 0.86 0.69 0.34 0.17 0.20 0.09 0.03 0.03
q17 q18
0.03 0.00
> round(apply(exam1[, 7:19], 1, mean), 2) # Helen got every good item wrong
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0.23 0.46 0.46 0.15 0.46 0.46 0.77 0.46 0.46 0.54 0.38 0.31 0.46 0.54 0.69 0.46
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
0.38 0.54 0.38 0.54 0.62 0.62 0.62 0.77 0.15 0.46 0.23 0.46 0.46 0.38 0.46 0.54
33 34 35
0.15 0.62 0.00
> exam1[35, ]
name sex q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12 q13 q14 q15 q16 q17 q18 id
35 Helen F 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 s35
> # remove items 1, 2, 3, and 18 because they provide no useful information
> names(exam1)
[1] "name" "sex" "q1" "q2" "q3" "q4" "q5" "q6" "q7" "q8"
[11] "q9" "q10" "q11" "q12" "q13" "q14" "q15" "q16" "q17" "q18"
[21] "id"
> exam1.1 <- exam1[, -c(3,4,5,20)]
> # remove Helen, who answered all good items incorrectly
> exam1.1 <- exam1[ -dim(exam1)[1], ]
\end{verbatim} % $
To run the conditional logit model \verb+clog()+, we need to
reshape the data from a wide format to a long format.
The \verb+reshape()+ command takes \verb+q4:q17+ and
organize them into a new variable called \verb+resp+.
Another new varible \verb+item+ is created to identify
questions 4 to 17.
\begin{verbatim}
> exam1.1 <- reshape(exam1.1, varying=list(paste("q", 4:17, sep="")),
+ timevar="item", v.names="resp", direction="long")
\end{verbatim}
In the equation the item difficulty parameters (the $\beta$'s)
have a minus sign before them. In \verb+clog()+ a parameter
with a minus sign is represented by a dummy variable of -1.
For example, in the long form the items are identified by the
variable \verb+item+. An excerpt of the full data shows that
Martha answered item 1 correctly, Elsie also answered item 1
correctly. The answers of Martha and Elsie for item 1 are followed
by answers of Richard, Tracie, and Walter to the next item (item 2).
\begin{verbatim}
> exam1.1[34:38, ]
name sex id item resp
s33.1 Martha F s33 1 1
s34.1 Elsie F s34 1 1
s1.2 Richard M s1 2 1
s2.2 Tracie F s2 2 1
s3.2 Walter M s3 2 1
\end{verbatim}
Using the above data excerpt as an example, we can create dummy
variables for items 1 and 2. The R command is
\verb+tt <- factor(exam1.1[34:38, "item"]); dummy <- diag(nlevels(tt))[tt, ]+
and you have:
\begin{verbatim}
> cbind(exam1.1[34:38, ], 0 - dummy)
name sex id item resp 1 2
s34.1 Elsie F s34 1 1 -1 0
s35.1 Helen F s35 1 1 -1 0
s1.2 Richard M s1 2 1 0 -1
s2.2 Tracie F s2 2 1 0 -1
s3.2 Walter M s3 2 1 0 -1
\end{verbatim}
Note that the last two columns are added with \verb+cbind()+ and they
can be used to code the $\beta$ parameters for items 1 and 2, respectively.
When this is applied to the 476 rows in the data set, it produces a
dummy variable matrix of 476 rows and 14 columns. Each column of the
matrix will be used to estimate the item difficulty associated
with each of questions 4 to 17.
\begin{verbatim}
> ex.i.dummy <- diag(nlevels(factor(exam1.1$item)))[factor(exam1.1$item), ]
> ex.i.dummy <- 0 - ex.i.dummy # turns (0, 1) into (0, -1)
> ex.i.dummy <- data.frame(ex.i.dummy, row.names=NULL)
> dimnames(ex.i.dummy)[[2]] <- paste("i", 4:17, sep="")
> dim(ex.i.dummy)
[1] 476 14
\end{verbatim}
Finally, a call to \verb+clogit()+ completes the calculations.
Worth noting is how the syntax is written. The \verb+strata()+
function tells \verb+clogit()+ to estimate item
difficulty associated with each dummy variable by holding
constant the influence across different persons.
The syntax also specifies
that the item difficulty paramaters
be estimated for all items except \verb+i4+. This is because
\verb+i4+ is the reference level when all dummy variables are
zero. The \verb+clogit()+ function is not able to define a
reference level if \verb+i4+ is added. The \verb+clogit()+
function with \verb+i4+ does not know what to do when all
dummy variables are zero.
\begin{verbatim}
> attach(ex.i.dummy) # attach the dummy variables so I can call i4:i17
> # item 4 is the reference
> exam2.clog<- clogit(resp ~ i5+i6+i7+i8+i9+i10+i11+i12+i13+i14+
+ i15+i16+i17+ strata(id), data=exam1.1)
> summary(exam2.clog)
Call:
clogit(resp ~ i5 + i6 + i7 + i8 + i9 + i10 + i11 + i12 + i13 +
i14 + i15 + i16 + i17 + strata(id), data = exam1.1)
n= 476
coef exp(coef) se(coef) z p
i5 0.514 1.67 1.025 0.501 6.2e-01
i6 0.923 2.52 0.991 0.931 3.5e-01
i7 0.514 1.67 1.025 0.501 6.2e-01
i8 1.875 6.52 0.955 1.964 5.0e-02
i9 0.923 2.52 0.991 0.931 3.5e-01
i10 2.598 13.44 0.944 2.752 5.9e-03
i11 4.517 91.54 0.954 4.736 2.2e-06
i12 5.792 327.53 1.028 5.632 1.8e-08
i13 5.537 253.95 1.010 5.484 4.1e-08
i14 6.825 920.74 1.136 6.007 1.9e-09
i15 8.095 3278.95 1.403 5.772 7.8e-09
i16 8.095 3278.95 1.403 5.772 7.8e-09
i17 8.095 3278.95 1.403 5.772 7.8e-09
exp(coef) exp(-coef) lower .95 upper .95
i5 1.67 0.597954 0.224 12.5
i6 2.52 0.397477 0.361 17.6
i7 1.67 0.597954 0.224 12.5
i8 6.52 0.153373 1.004 42.4
i9 2.52 0.397477 0.361 17.6
i10 13.44 0.074425 2.112 85.5
i11 91.54 0.010924 14.120 593.5
i12 327.53 0.003053 43.650 2457.7
i13 253.95 0.003938 35.104 1837.0
i14 920.74 0.001086 99.293 8538.0
i15 3278.95 0.000305 209.826 51240.2
i16 3278.95 0.000305 209.826 51240.2
i17 3278.95 0.000305 209.826 51240.2
Rsquare= 0.525 (max possible= 0.659 )
Likelihood ratio test= 355 on 13 df, p=0
Wald test = 108 on 13 df, p=0
Score (logrank) test = 277 on 13 df, p=0
\end{verbatim}
With the \verb+i4+ being the reference
level, the coeffiencts associated with
the dummy variables should be interpreted as the differences
in item difficulty between \verb+i4+ and each of
\verb+i5+ to \verb+i17+. For example, the \verb+0.514+ associated
with \verb+i5+ means that \verb+i5+ is more difficult (remember,
higher $\beta$ means higher difficulty) than \verb+i4+ by
\verb+0.514+. In addition, items 15, 16, and 17 are equally
difficult. They are more difficult than \verb+i4+ by
\verb+8.095+ logits. This translates to a reduction in the
probability of answering items 15, 16, and 17 correctly. On a
probability scale, this 8.095 logit is equivalent to 0.999
reduction in probability (try
\verb=exp(8.095) / ( 1 + exp(8.095))=).
The \verb+clogit()+ function in R is general-purpose routine designed
for fitting a conditional logit model. Thus it does not automatically
print out many other statistics specific to Rasch model. For example,
the \verb+summary()+ function does not automatically generate the
estimated ability level for each person. Other packages like MiniStep
does that for you. But \verb+clogit()+ in R provides estimates of
item difficulty, which are the most important
information for developing surveys or tests.
\subsection{Conjoint analysis}
In true conjoint analysis, we present a subject with stimuli made
by crossing at least two different variables. For example, in
one study, Baron and Andrea Gurmankin presented each subject with
56 items in a random order. The items consisted of each of each
of eight medical conditions (ranging from a wart to immediate
death) at each of seven probability levels, and the subjects
provided a badness rating for each of the 56 items. The
assumption of conjoint analysis is that the subject behaves as if
he represents the disutility of each condition as a number, and
the probability as another number, adds the two numbers, and
transforms the result monotonically into the response scale
provided. The representation of probability is a monotonic
function of the given probability.
When we analyze the data, we try to recover the three
transformations so as to get the best fit assuming that this
model is true. The three transformations are the assignment of
numbers to the probabilities, the assignment of numbers to the
conditions, and the function relating the result to the response
scale. (In this case, if the subject followed expected-utility
theory, the probability would be transformed logarithmically, so
that the additive representation corresponded to multiplication.)
R does not have a conjoint analysis package as such. But the
Acepack package contains a function called \texttt{ace()}, for
``alternating conditional expectations,'' which does essentially
what we want. It maximizes the variance in the dependent
variable (the response) that is explained by the predictors
(probability and condition), by using an iterative process. Here
is an example in which the response is called \texttt{bad}, which
is a matrix in which the rows are subjects, and within each row
the probabilities are in groups of 8, with the conditions
repeated in each group.
\begin{verbatim}
probs <- rep(1:7,rep(8,7))
conds <- gl(8,1,56)
cnames <- c("wart","toe","deaf1","leg1","leg2","blind","bbdd","death")
pnames <- c(".001",".0032",".01",".032",".1",".32","1")
c.wt.ace <- matrix(0,ns,8) # resulting numbers for conditions
p.wt.ace <- matrix(0,ns,7) # resulting transformed probabilities
bad.ace <- matrix(0,ns,56) # transformed responses
for (i in 1:ns) # fit the model for each subject
{av1 <- ace(cbind(probs,conds),bad[i,],lin=1)
bad.ace[i,] <- av1$ty
p.wt.ace[i,] <- av1$tx[8*0:6+1,1]
c.wt.ace[i,] <- av1$tx[1:8,2]
}
\end{verbatim} %$
In the end, the matrices \texttt{p.wt.ace, c.wt.ace} and
\texttt{bad.ace} should have the transformed numbers, one row per
subject.
\subsection{Imputation of missing data}
Schafer and Graham (2002) provide a good review of methods for
dealing with missing data. R provides many of the methods that
they discuss. One method is multiple imputation, which is found
in the Hmisc package. Each missing datum is inferred repeated
from different samples of other variables, and the repeated
inferences are used to determine the error. It turns out that
this method works best with the \texttt{ols()} function from the
Design package rather than with (the otherwise equivalent)
\texttt{lm()} function. Here is an example, using the data set
\texttt{t1}.
\begin{verbatim}
library(Hmisc)
f <- aregImpute(~v1+v2+v3+v4, n.impute=20,
fweighted=.2, defaultLinear=T, data=t1)
library(Design)
fmp <- fit.mult.impute(v1~v2+v3, ols, f, data=t1)
summary(fmp)
\end{verbatim}
The first command (\texttt{f}) imputes missing values in all four
variables, using the other three for each variable. The second
command (\texttt{fmp}) estimates a regression model in which
\texttt{v1} is predicted from two of the remaining variables. A
variable can be used (and should be used, if it is useful) for
imputation even when it is not used in the model.
\section{References}
{\everypar = {\hangafter=1 \parindent 0pt \hangindent 1.5em}
\noindent
Chambers, J. M., \& Hastie, T. J. (1992). \textit{Statistical
models in S.} Pacific Grove, CA: Wadsworth \& Brooks Cole
Advanced Books and Software.
Clark, H. H. (1973). The language-as-fixed-effect fallacy: A
critique of language statistics in psychological research.
\textit{Journal of Verbal Learning \& Verbal Behavior, 12}, 335--359.
Hays, W. L. (1988, 4th ed.) \textit{Statistics.} New York: Holt,
Rinehart and Winston.
Hoaglin, D. C., Mosteller, F., \& Tukey, J. W. (Eds.)
(1983). \textit{Understanding robust and exploratory data
analysis.} New York: Wiley.
Lord, F. M., \& Novick, M. R. (1968). {\it Statistical theories
of mental test scores}. Reading, MA: Addison-Wesley.
Maxwell, S. E. \& Delaney, H. D. (1990) \textit{Designing Experiments and
Analyzing Data: A model comparison perspective.} Pacific
Grove, CA: Brooks/Cole.
Schafer, J. L., \& Graham, J. W. (2002). Missing data: Our view
of the state of the art. \texttt{Psychological Methods, 7},
147--177.
Stevens, J. (1992, 2nd ed) \textit{Applied Multivariate Statistics
for the Social Sciences.} Hillsdale, NJ: Erlbaum.
Venables, W. N., \& Ripley, B. D. (1999). \textit{Modern applied
statistics with S--PLUS} (3rd Ed.). New York: Springer.
}
\end{document}