% \VignetteIndexEntry{Overview of the psych package}
% \VignettePackage{psych}
% \VignetteKeywords{multivariate}
% \VignetteKeyword{models}
% \VignetteKeyword{Hplot}
%\VignetteDepends{psych}
%\documentclass[doc]{apa}
\documentclass[11pt]{article}
%\documentclass[11pt]{amsart}
\usepackage{geometry} % See geometry.pdf to learn the layout options. There are lots.
\geometry{letterpaper} % ... or a4paper or a5paper or ...
%\geometry{landscape} % Activate for for rotated page geometry
\usepackage[parfill]{parskip} % Activate to begin paragraphs with an empty line rather than an indent
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{epstopdf}
\usepackage{mathptmx}
\usepackage{helvet}
\usepackage{courier}
\usepackage{epstopdf}
\usepackage{makeidx} % allows index generation
\usepackage[authoryear,round]{natbib}
\usepackage{gensymb}
\usepackage{longtable}
%\usepackage{geometry}
\usepackage{amssymb}
\usepackage{amsmath}
%\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}
\usepackage{Sweave}
%\usepackage{/Volumes/'Macintosh HD'/Library/Frameworks/R.framework/Versions/2.13/Resources/share/texmf/tex/latex/Sweave}
%\usepackage[ae]{Rd}
%\usepackage[usenames]{color}
%\usepackage{setspace}
\bibstyle{apacite}
\bibliographystyle{apa} %this one plus author year seems to work?
%\usepackage{hyperref}
\usepackage[colorlinks=true,citecolor=blue]{hyperref} %this makes reference links hyperlinks in pdf!
\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}
\usepackage{multicol} % used for the two-column index
\usepackage[bottom]{footmisc}% places footnotes at page bottom
\let\proglang=\textsf
\newcommand{\R}{\proglang{R}}
%\newcommand{\pkg}[1]{{\normalfont\fontseries{b}\selectfont #1}}
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\fun}[1]{{\texttt{#1}\index{#1}\index{R function!#1}}}
\newcommand{\pfun}[1]{{\texttt{#1}\index{#1}\index{R function!#1}\index{R function!psych package!#1}}}\newcommand{\Rc}[1]{{\texttt{#1}}} %R command same as Robject
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpkg}[1]{{\textit{#1}\index{#1}\index{R package!#1}}} %different from pkg - which is better?
\newcommand{\iemph}[1]{{\emph{#1}\index{#1}}}
\newcommand{\wrc}[1]{\marginpar{\textcolor{blue}{#1}}} %bill's comments
\newcommand{\wra}[1]{\textcolor{blue}{#1}} %bill's comments
\newcommand{\ve}[1]{{\textbf{#1}}} %trying to get a vector command
\makeindex % used for the subject index
\title{How to do a factor analysis with the psych package}
\author{William Revelle\\Department of Psychology\\Northwestern University}
%\affiliation{Northwestern University}
%\acknowledgements{Written to accompany the psych package. Comments should be directed to William Revelle \\ \url{revelle@northwestern.edu}}
%\date{} % Activate to display a given date or no date
\begin{document}
\input{factor-concordance}
\maketitle
\tableofcontents
\newpage
\subsection{Jump starting the \Rpkg{psych} package to do factor analysis --a guide for the impatient}
You have installed \Rpkg{psych} (section \ref{sect:starting}) and you want to use it for a factor analysis without reading much more. What should you do?
If you have already read the \href{http://personality-project.org/r/psych/HowTo/getting_started.pdf}{getting started} documentation and have \R{} \citep{R} and the \Rpkg{psych} package \citep{psych} installed then just
\begin{enumerate}
\item Activate the \Rpkg{psych} package:
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
library(psych)
\end{Sinput}
\end{Schunk}
\end{scriptsize}
\item Input your data (section \ref{sect:read}). Go to your friendly text editor or data manipulation program (e.g., Excel) and copy the data to the clipboard. Include a first line that has the variable labels. Paste it into \Rpkg{psych} using the \pfun{read.clipboard.tab} command:
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
myData <- read.clipboard.tab()
\end{Sinput}
\end{Schunk}
\end{scriptsize}
\item Make sure that what you just read is right. Describe it (section~\ref{sect:describe}) and perhaps look at the first and last few lines:
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
describe(myData)
headTail(myData)
\end{Sinput}
\end{Schunk}
\end{scriptsize}
\item Test for the number of factors in your data using parallel analysis (\pfun{fa.parallel}, section \ref{sect:fa.parallel}) or Very Simple Structure (\pfun{vss}, \ref{sect:vss}) .
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
fa.parallel(myData)
vss(myData)
\end{Sinput}
\end{Schunk}
\end{scriptsize}
\item Factor analyze (see section \ref{sect:fa}) the data with a specified number of factors (the default is 1), the default method is minimum residual, the default rotation for more than one factor is oblimin. There are many more possibilities (see sections \ref{sect:minres}-\ref{sect:wls}). Compare the solution to a hierarchical cluster analysis using the ICLUST algorithm \citep{revelle:iclust} (see section \ref{sect:iclust}). Also consider a hierarchical factor solution to find coefficient $\omega$ (see \ref{sect:omega}).
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
fa(myData)
iclust(myData)
omega(myData)
\end{Sinput}
\end{Schunk}
\end{scriptsize}
\end{enumerate}
You might find reading this entire \href{http://personality-project.org/r/psych/vignettes/overview.pdf} {over view vignette} helpful to get a broader understanding of what can be done in \R{} using the \Rpkg{psych}. Remember that the help command (?) is available for every function. Try running the examples for each help page.
\newpage
\section{Overview of this and related documents}
This is one of a number of ``How To use the \Rpkg{psych} package to do X'' documents. The entire set may be found at \href{personality-project.org/r/psych/index.html}{personality-project.org/r/psych}. Some of these are extracted from the \href{http://personality-project.org/r/psych/vignettes/overview.pdf}{overview of psych} vignette.
\section{Input the data}
A number of \Rpkg{psych} functions facilitate the entry of data and finding basic descriptive statistics.
Remember, to run any of the \Rpkg{psych} functions, it is necessary to make the package active by using the \fun{library} command:
\begin{Schunk}
\begin{Sinput}
> library(psych)
\end{Sinput}
\end{Schunk}
The other packages, once installed, will be called automatically by \Rpkg{psych}.
It is possible to automatically load \Rpkg{psych} and other functions by creating and then saving a ``.First" function: e.g.,
\begin{Schunk}
\begin{Sinput}
.First <- function(x) {library(psych)}
\end{Sinput}
\end{Schunk}
\subsection{Data input from the clipboard}
\label{sect:read}
There are of course many ways to enter data into \R. Reading from a local file using \fun{read.table} is perhaps the most preferred. However, many users will enter their data in a text editor or spreadsheet program and then want to copy and paste into \R{}. This may be done by using \fun{read.table} and specifying the input file as ``clipboard" (PCs) or ``pipe(pbpaste)" (Macs). Alternatively, the \pfun{read.clipboard} set of functions are perhaps more user friendly:
\begin{description}
\item [\pfun{read.clipboard}] is the base function for reading data from the clipboard.
\item [\pfun{read.clipboard.csv}] for reading text that is comma delimited.
\item [\pfun{read.clipboard.tab}] for reading text that is tab delimited (e.g., copied directly from an Excel file).
\item [\pfun{read.clipboard.lower}] for reading input of a lower triangular matrix with or without a diagonal. The resulting object is a square matrix.
\item [\pfun{read.clipboard.upper}] for reading input of an upper triangular matrix.
\item[\pfun{read.clipboard.fwf}] for reading in fixed width fields (some very old data sets)
\end{description}
For example, given a data set copied to the clipboard from a spreadsheet, just enter the command
\begin{Schunk}
\begin{Sinput}
> my.data <- read.clipboard()
\end{Sinput}
\end{Schunk}
This will work if every data field has a value and even missing data are given some values (e.g., NA or -999). If the data were entered in a spreadsheet and the missing values were just empty cells, then the data should be read in as a tab delimited or by using the \pfun{read.clipboard.tab} function.
\begin{Schunk}
\begin{Sinput}
> my.data <- read.clipboard(sep="\t") #define the tab option, or
> my.tab.data <- read.clipboard.tab() #just use the alternative function
\end{Sinput}
\end{Schunk}
For the case of data in fixed width fields (some old data sets tend to have this format), copy to the clipboard and then specify the width of each field (in the example below, the first variable is 5 columns, the second is 2 columns, the next 5 are 1 column the last 4 are 3 columns).
\begin{Schunk}
\begin{Sinput}
> my.data <- read.clipboard.fwf(widths=c(5,2,rep(1,5),rep(3,4))
\end{Sinput}
\end{Schunk}
\subsection{Basic descriptive statistics}
\label{sect:describe}
Before doing any analysis, it is important to make sure that your data are what you think they are. You should examine the basic descriptive statistics using the \pfun{describe} function.
Consider the data set \pfun{sat.act} which includes data from 700 web based participants on 3 demographic variables and 3 ability measures.
\begin{description}
\item[\pfun{describe}] reports means, standard deviations, medians, min, max, range, skew, kurtosis and standard errors for integer or real data. Non-numeric data, although the statistics are meaningless, will be treated as if numeric (based upon the categorical coding of the data), and will be flagged with an *.
\item[\pfun{describeBy}] reports descriptive statistics broken down by some categorizing variable (e.g., gender, age, etc.)
\end{description}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> library(psych)
> data(sat.act)
> describe(sat.act) #basic descriptive statistics
\end{Sinput}
\begin{Soutput}
var n mean sd median trimmed mad min max range skew kurtosis se
gender 1 700 1.65 0.48 2 1.68 0.00 1 2 1 -0.61 -1.62 0.02
education 2 700 3.16 1.43 3 3.31 1.48 0 5 5 -0.68 -0.07 0.05
age 3 700 25.59 9.50 22 23.86 5.93 13 65 52 1.64 2.42 0.36
ACT 4 700 28.55 4.82 29 28.84 4.45 3 36 33 -0.66 0.53 0.18
SATV 5 700 612.23 112.90 620 619.45 118.61 200 800 600 -0.64 0.33 4.27
SATQ 6 687 610.22 115.64 620 617.25 118.61 200 800 600 -0.59 -0.02 4.41
\end{Soutput}
\end{Schunk}
\end{scriptsize}
%These data may then be analyzed by groups defined in a logical statement or by some other variable. E.g., break down the descriptive data for males or females. These descriptive data can also be seen graphically using the \pfun{error.bars.by} function (Figure~\ref{fig:error.bars}). By setting skew=FALSE and ranges=FALSE, the output is limited to the most basic statistics.
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> #basic descriptive statistics by a grouping variable.
%> describeBy(sat.act,sat.act$gender,skew=FALSE,ranges=FALSE)
%\end{Sinput}
%\begin{Soutput}
%group: 1
% var n mean sd se
%gender 1 247 1.00 0.00 0.00
%education 2 247 3.00 1.54 0.10
%age 3 247 25.86 9.74 0.62
%ACT 4 247 28.79 5.06 0.32
%SATV 5 247 615.11 114.16 7.26
%SATQ 6 245 635.87 116.02 7.41
%--------------------------------------------------------------------------
%group: 2
% var n mean sd se
%gender 1 453 2.00 0.00 0.00
%education 2 453 3.26 1.35 0.06
%age 3 453 25.45 9.37 0.44
%ACT 4 453 28.42 4.69 0.22
%SATV 5 453 610.66 112.31 5.28
%SATQ 6 442 596.00 113.07 5.38
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%
%The output from the \pfun{describeBy} function can be forced into a matrix form for easy analysis by other programs. In addition, describeBy can group by several grouping variables at the same time.
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> sa.mat <- describeBy(sat.act,list(sat.act$gender,sat.act$education),
%+ skew=FALSE,ranges=FALSE,mat=TRUE)
%> headTail(sa.mat)
%\end{Sinput}
%\begin{Soutput}
% item group1 group2 var n mean sd se
%gender1 1 1 0 1 27 1 0 0
%gender2 2 2 0 1 30 2 0 0
%gender3 3 1 1 1 20 1 0 0
%gender4 4 2 1 1 25 2 0 0
%... ... ... ... ... ...
%SATQ9 69 1 4 6 51 635.9 104.12 14.58
%SATQ10 70 2 4 6 86 597.59 106.24 11.46
%SATQ11 71 1 5 6 46 657.83 89.61 13.21
%SATQ12 72 2 5 6 93 606.72 105.55 10.95
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%\subsubsection{Basic data cleaning using \pfun{scrub}}
%\label{sect:scrub}
%If, after describing the data it is apparent that there were data entry errors that need to be globally replaced with NA, or only certain ranges of data will be analyzed, the data can be ``cleaned" using the \pfun{scrub} function.
%
%Consider a data set of 10 rows of 12 columns with values from 1 - 120. All values of columns 3 - 5 that are less than 30, 40, or 50 respectively, or greater than 70 in any of the three columns will be replaced with NA. In addition, any value exactly equal to 45 will be set to NA. (max and isvalue are set to one value here, but they could be a different value for every column).
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> x <- matrix(1:120,ncol=10,byrow=TRUE)
%> colnames(x) <- paste('V',1:10,sep='')
%> new.x <- scrub(x,3:5,min=c(30,40,50),max=70,isvalue=45,newvalue=NA)
%> new.x
%\end{Sinput}
%\begin{Soutput}
% V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
% [1,] 1 2 NA NA NA 6 7 8 9 10
% [2,] 11 12 NA NA NA 16 17 18 19 20
% [3,] 21 22 NA NA NA 26 27 28 29 30
% [4,] 31 32 33 NA NA 36 37 38 39 40
% [5,] 41 42 43 44 NA 46 47 48 49 50
% [6,] 51 52 53 54 55 56 57 58 59 60
% [7,] 61 62 63 64 65 66 67 68 69 70
% [8,] 71 72 NA NA NA 76 77 78 79 80
% [9,] 81 82 NA NA NA 86 87 88 89 90
%[10,] 91 92 NA NA NA 96 97 98 99 100
%[11,] 101 102 NA NA NA 106 107 108 109 110
%[12,] 111 112 NA NA NA 116 117 118 119 120
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%Note that the number of subjects for those columns has decreased, and the minimums have gone up but the maximums down. Data cleaning and examination for outliers should be a routine part of any data analysis.
%
%\subsubsection{Recoding categorical variables into dummy coded variables}
%Sometimes categorical variables (e.g., college major, occupation, ethnicity) are to be analyzed using correlation or regression. To do this, one can form ``dummy codes'' which are merely binary variables for each category. This may be done using \pfun{dummy.code}. Subsequent analyses using these dummy coded variables may be using \pfun{biserial} or point biserial (regular Pearson r) to show effect sizes and may be plotted in e.g., \pfun{spider} plots.
%
%\subsection{Simple descriptive graphics}
%Graphic descriptions of data are very helpful both for understanding the data as well as communicating important results. Scatter Plot Matrices (SPLOMS) using the \pfun{pairs.panels} function are useful ways to look for strange effects involving outliers and non-linearities. \pfun{error.bars.by} will show group means with 95\% confidence boundaries.
%
%\subsubsection{Scatter Plot Matrices}
%Scatter Plot Matrices (SPLOMS) are very useful for describing the data. The \pfun{pairs.panels} function, adapted from the help menu for the \fun{pairs} function produces xy scatter plots of each pair of variables below the diagonal, shows the histogram of each variable on the diagonal, and shows the \iemph{lowess} locally fit regression line as well. An ellipse around the mean with the axis length reflecting one standard deviation of the x and y variables is also drawn. The x axis in each scatter plot represents the column variable, the y axis the row variable (Figure~\ref{fig:pairs.panels}). When plotting many subjects, it is both faster and cleaner to set the plot character (pch) to be '.'. (See Figure~\ref{fig:pairs.panels} for an example.)
%
%\begin{description}
%\label{sect:pairs}
%\item[\pfun{pairs.panels} ] will show the pairwise scatter plots of all the variables as well as histograms, locally smoothed regressions, and the Pearson correlation. When plotting many data points (as in the case of the sat.act data, it is possible to specify that the plot character is a period to get a somewhat cleaner graphic.
%\end{description}
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> png( 'pairspanels.png' )
%> pairs.panels(sat.act,pch='.')
%> dev.off()
%\end{Sinput}
%\begin{Soutput}
%null device
% 1
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%\includegraphics{pairspanels}
%\caption{Using the \pfun{pairs.panels} function to graphically show relationships. The x axis in each scatter plot represents the column variable, the y axis the row variable. Note the extreme outlier for the ACT. The plot character was set to a period (pch='.') in order to make a cleaner graph. }
%\label{fig:pairs.panels}
%\end{center}
%\end{figure}
%
%Another example of \pfun{pairs.panels} is to show differences between experimental groups. Consider the data in the \pfun{affect} data set. The scores reflect post test scores on positive and negative affect and energetic and tense arousal. The colors show the results for four movie conditions: depressing, frightening movie, neutral, and a comedy.
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> png('affect.png')
%> pairs.panels(affect[14:17],bg=c("red","black","white","blue")[affect$Film],pch=21,
%+ main="Affect varies by movies ")
%> dev.off()
%\end{Sinput}
%\begin{Soutput}
%null device
% 1
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%\includegraphics{affect}
%\caption{Using the \pfun{pairs.panels} function to graphically show relationships. The x axis in each scatter plot represents the column variable, the y axis the row variable. The coloring represent four different movie conditions. }
%\label{fig:pairs.panels2}
%\end{center}
%\end{figure}
%
%\subsubsection{Means and error bars}
%\label{sect:errorbars}
%Additional descriptive graphics include the ability to draw \iemph{error bars} on sets of data, as well as to draw error bars in both the x and y directions for paired data. These are the functions
%
%\begin{description}
%\item [\pfun{error.bars}] show the 95 \% confidence intervals for each variable in a data frame or matrix. These errors are based upon normal theory and the standard errors of the mean. Alternative options include +/- one standard deviation or 1 standard error. If the data are repeated measures, the error bars will be reflect the between variable correlations.
%\item [\pfun{error.bars.by}] does the same, but grouping the data by some condition.
%\item [\pfun{error.crosses}] draw the confidence intervals for an x set and a y set of the same size.
%\end{description}
%
%The use of the \pfun{error.bars.by} function allows for graphic comparisons of different groups (see Figure~\ref{fig:error.bars}). Five personality measures are shown as a function of high versus low scores on a ``lie" scale. People with higher lie scores tend to report being more agreeable, conscientious and less neurotic than people with lower lie scores. The error bars are based upon normal theory and thus are symmetric rather than reflect any skewing in the data.
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> data(epi.bfi)
%> error.bars.by(epi.bfi[,6:10],epi.bfi$epilie<4)
%\end{Sinput}
%\end{Schunk}
%\includegraphics{factor-007}
%\end{scriptsize}
%\caption{Using the \pfun{error.bars.by} function shows that self reported personality scales on the Big Five Inventory vary as a function of the Lie scale on the EPI. }
%\label{fig:error.bars}
%\end{center}
%\end{figure}
%
%Although not recommended, it is possible to use the \pfun{error.bars} function to draw bar graphs with associated error bars. (This kind of`\iemph{dynamite plot} (Figure~\ref{fig:dynamite}) can be very misleading in that the scale is arbitrary. Go to a discussion of the problems in presenting data this way at \url{http://emdbolker.wikidot.com/blog:dynamite}. In the example shown, note that the graph starts at 0, although is out of the range. This is a function of using bars, which always are assumed to start at zero. Consider other ways of showing your data.
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> error.bars.by(sat.act[5:6],sat.act$gender,bars=TRUE,
%+ labels=c("Male","Female"),ylab="SAT score",xlab="")
%\end{Sinput}
%\end{Schunk}
%\includegraphics{factor-008}
%\end{scriptsize}
%\caption{A ``Dynamite plot" of SAT scores as a function of gender is one way of misleading the reader. By using a bar graph, the range of scores is ignored. Bar graphs start from 0. }
%\label{fig:dynamite}
%\end{center}
%\end{figure}
%
%
%\subsubsection{Two dimensional displays of means and errors}
%Yet another way to display data for different conditions is to use the \pfun{errorCrosses} function. For instance, the effect of various movies on both ``Energetic Arousal'' and ``Tense Arousal'' can be seen in one graph and compared to the same movie manipulations on ``Positive Affect'' and ``Negative Affect''. Note how Energetic Arousal is increased by three of the movie manipulations, but that Positive Affect increases following the Happy movie only.
%
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> op <- par(mfrow=c(1,2))
%> data(affect)
%> colors <- c("black","red","white","blue")
%> films <- c("Sad","Horror","Neutral","Happy")
%> affect.stats <- errorCircles("EA2","TA2",data=affect,group="Film",labels=films,xlab="Energetic Arousal",ylab="Tense Arousal",ylim=c(10,22),xlim=c(8,20),pch=16,cex=2,col=colors,
%+ main =' Movies effect on arousal')
%> errorCircles("PA2","NA2",data=affect.stats,labels=films,xlab="Positive Affect",ylab="Negative Affect",pch=16,cex=2,col=colors,
%+ main ="Movies effect on affect")
%> op <- par(mfrow=c(1,1))
%\end{Sinput}
%\end{Schunk}
%\includegraphics{factor-009}
%\end{scriptsize}
%\caption{The use of the \pfun{errorCircles} function allows for two dimensional displays of means and error bars. The first call to \pfun{errorCircles} finds descriptive statistics for the \iemph{affect} data.frame based upon the grouping variable of Film. These data are returned and then used by the second call which examines the effect of the same grouping variable upon different measures. The size of the circles represent the relative sample sizes for each group. The data are from the PMC lab and reported in \cite{smillie:jpsp}.}
%\label{fig:errorCircles}
%\end{center}
%\end{figure}
%
%\clearpage
%\subsubsection{Back to back histograms}
%The \pfun{bi.bars} function summarize the characteristics of two groups (e.g., males and females) on a second variable (e.g., age) by drawing back to back histograms (see Figure~\ref{fig:bibars}).
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> data(bfi)
%> with(bfi,{bi.bars(age,gender,ylab="Age",main="Age by males and females")})
%\end{Sinput}
%\end{Schunk}
%\includegraphics{factor-010}
%\end{scriptsize}
%\caption{A bar plot of the age distribution for males and females shows the use of \pfun{bi.bars}. The data are males and females from 2800 cases collected using the \iemph{SAPA} procedure and are available as part of the \pfun{bfi} data set. }
%\label{fig:bibars}
%\end{center}
%\end{figure}
%
%\clearpage
%\subsubsection{Correlational structure}
%\label{sect:lowerCor}
%There are many ways to display correlations. Tabular displays are probably the most common. The output from the \fun{cor} function in core R is a rectangular matrix. \pfun{lowerMat} will round this to (2) digits and then display as a lower off diagonal matrix. \pfun{lowerCor} calls \fun{cor} with \emph{use=`pairwise', method=`pearson'} as default values and returns (invisibly) the full correlation matrix and displays the lower off diagonal matrix.
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> lowerCor(sat.act)
%\end{Sinput}
%\begin{Soutput}
% gendr edctn age ACT SATV SATQ
%gender 1.00
%education 0.09 1.00
%age -0.02 0.55 1.00
%ACT -0.04 0.15 0.11 1.00
%SATV -0.02 0.05 -0.04 0.56 1.00
%SATQ -0.17 0.03 -0.03 0.59 0.64 1.00
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%When comparing results from two different groups, it is convenient to display them as one matrix, with the results from one group below the diagonal, and the other group above the diagonal. Use \pfun{lowerUpper} to do this:
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> female <- subset(sat.act,sat.act$gender==2)
%> male <- subset(sat.act,sat.act$gender==1)
%> lower <- lowerCor(male[-1])
%\end{Sinput}
%\begin{Soutput}
% edctn age ACT SATV SATQ
%education 1.00
%age 0.61 1.00
%ACT 0.16 0.15 1.00
%SATV 0.02 -0.06 0.61 1.00
%SATQ 0.08 0.04 0.60 0.68 1.00
%\end{Soutput}
%\begin{Sinput}
%> upper <- lowerCor(female[-1])
%\end{Sinput}
%\begin{Soutput}
% edctn age ACT SATV SATQ
%education 1.00
%age 0.52 1.00
%ACT 0.16 0.08 1.00
%SATV 0.07 -0.03 0.53 1.00
%SATQ 0.03 -0.09 0.58 0.63 1.00
%\end{Soutput}
%\begin{Sinput}
%> both <- lowerUpper(lower,upper)
%> round(both,2)
%\end{Sinput}
%\begin{Soutput}
% education age ACT SATV SATQ
%education NA 0.52 0.16 0.07 0.03
%age 0.61 NA 0.08 -0.03 -0.09
%ACT 0.16 0.15 NA 0.53 0.58
%SATV 0.02 -0.06 0.61 NA 0.63
%SATQ 0.08 0.04 0.60 0.68 NA
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%It is also possible to compare two matrices by taking their differences and displaying one (below the diagonal) and the difference of the second from the first above the diagonal:
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> diffs <- lowerUpper(lower,upper,diff=TRUE)
%> round(diffs,2)
%\end{Sinput}
%\begin{Soutput}
% education age ACT SATV SATQ
%education NA 0.09 0.00 -0.05 0.05
%age 0.61 NA 0.07 -0.03 0.13
%ACT 0.16 0.15 NA 0.08 0.02
%SATV 0.02 -0.06 0.61 NA 0.05
%SATQ 0.08 0.04 0.60 0.68 NA
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%\subsubsection{Heatmap displays of correlational structure}
%\label{sect:corplot}
%Perhaps a better way to see the structure in a correlation matrix is to display a \emph{heat map} of the correlations. This is just a matrix color coded to represent the magnitude of the correlation. This is useful when considering the number of factors in a data set. Consider the \pfun{Thurstone} data set which has a clear 3 factor solution (Figure~\ref{fig:cor.plot}) or a simulated data set of 24 variables with a circumplex structure (Figure~\ref{fig:cor.plot.circ}). The color coding represents a ``heat map'' of the correlations, with darker shades of red representing stronger negative and darker shades of blue stronger positive correlations. As an option, the value of the correlation can be shown.
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> png('corplot.png')
%> cor.plot(Thurstone,numbers=TRUE,main="9 cognitive variables from Thurstone")
%> dev.off()
%\end{Sinput}
%\begin{Soutput}
%null device
% 1
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%\includegraphics{corplot.png}
%\caption{The structure of correlation matrix can be seen more clearly if the variables are grouped by factor and then the correlations are shown by color. By using the 'numbers' option, the values are displayed as well. }
%\label{fig:cor.plot}
%\end{center}
%\end{figure}
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> png('circplot.png')
%> circ <- sim.circ(24)
%> r.circ <- cor(circ)
%> cor.plot(r.circ,main='24 variables in a circumplex')
%> dev.off()
%\end{Sinput}
%\begin{Soutput}
%null device
% 1
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%\includegraphics{circplot.png}
%\caption{Using the cor.plot function to show the correlations in a circumplex. Correlations are highest near the diagonal, diminish to zero further from the diagonal, and the increase again towards the corners of the matrix. Circumplex structures are common in the study of affect.}
%\label{fig:cor.plot.circ}
%\end{center}
%\end{figure}
%
%Yet another way to show structure is to use ``spider'' plots. Particularly if variables are ordered in some meaningful way (e.g., in a circumplex), a spider plot will show this structure easily. This is just a plot of the magnitude of the correlation as a radial line, with length ranging from 0 (for a correlation of -1) to 1 (for a correlation of 1). (See Figure~\ref{fig:cor.plot.spider}).
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> png('spider.png')
%> op<- par(mfrow=c(2,2))
%> spider(y=c(1,6,12,18),x=1:24,data=r.circ,fill=TRUE,main="Spider plot of 24 circumplex variables")
%> op <- par(mfrow=c(1,1))
%> dev.off()
%\end{Sinput}
%\begin{Soutput}
%null device
% 1
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%\includegraphics{spider.png}
%\caption{A spider plot can show circumplex structure very clearly. Circumplex structures are common in the study of affect.}
%\label{fig:cor.plot.spider}
%\end{center}
%\end{figure}
%
%
%\subsection{Testing correlations}
%\label{sect:corr.test}
%Correlations are wonderful descriptive statistics of the data but some people like to test whether these correlations differ from zero, or differ from each other. The \fun{cor.test} function (in the \Rpkg{stats} package) will test the significance of a single correlation, and the \fun{rcorr} function in the \Rpkg{Hmisc} package will do this for many correlations. In the \Rpkg{psych} package, the \pfun{corr.test} function reports the correlation (Pearson, Spearman, or Kendall) between all variables in either one or two data frames or matrices, as well as the number of observations for each case, and the (two-tailed) probability for each correlation. Unfortunately, these probability values have not been corrected for multiple comparisons and so should be taken with a great deal of salt. Thus, in \pfun{corr.test} and \pfun{corr.p} the raw probabilities are reported below the diagonal and the probabilities adjusted for multiple comparisons using (by default) the Holm correction are reported above the diagonal (Table~\ref{tab:corr.test}). (See the \fun{p.adjust} function for a discussion of \cite{holm:79} and other corrections.)
%
%\begin{table}[htdp]
%\caption{The \pfun{corr.test} function reports correlations, cell sizes, and raw and adjusted probability values. \pfun{corr.p} reports the probability values for a correlation matrix. By default, the adjustment used is that of \cite{holm:79}.}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> corr.test(sat.act)
%\end{Sinput}
%\begin{Soutput}
%Call:corr.test(x = sat.act)
%Correlation matrix
% gender education age ACT SATV SATQ
%gender 1.00 0.09 -0.02 -0.04 -0.02 -0.17
%education 0.09 1.00 0.55 0.15 0.05 0.03
%age -0.02 0.55 1.00 0.11 -0.04 -0.03
%ACT -0.04 0.15 0.11 1.00 0.56 0.59
%SATV -0.02 0.05 -0.04 0.56 1.00 0.64
%SATQ -0.17 0.03 -0.03 0.59 0.64 1.00
%Sample Size
% gender education age ACT SATV SATQ
%gender 700 700 700 700 700 687
%education 700 700 700 700 700 687
%age 700 700 700 700 700 687
%ACT 700 700 700 700 700 687
%SATV 700 700 700 700 700 687
%SATQ 687 687 687 687 687 687
%Probability values (Entries above the diagonal are adjusted for multiple tests.)
% gender education age ACT SATV SATQ
%gender 0.00 0.17 1.00 1.00 1 0
%education 0.02 0.00 0.00 0.00 1 1
%age 0.58 0.00 0.00 0.03 1 1
%ACT 0.33 0.00 0.00 0.00 0 0
%SATV 0.62 0.22 0.26 0.00 0 0
%SATQ 0.00 0.36 0.37 0.00 0 0
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%\label{tab:corr.test}
%\end{table}%
%
%
%Testing the difference between any two correlations can be done using the \pfun{r.test} function. The function actually does four different tests (based upon an article by \cite{steiger:80b}, depending upon the input:
%
%1) For a sample size n, find the t and p value for a single correlation as well as the confidence interval.
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> r.test(50,.3)
%\end{Sinput}
%\begin{Soutput}
%Correlation tests
%Call:r.test(n = 50, r12 = 0.3)
%Test of significance of a correlation
% t value 2.18 with probability < 0.034
% and confidence interval 0.02 0.53
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%2) For sample sizes of n and n2 (n2 = n if not specified) find the z of the difference between the z transformed correlations divided by the standard error of the difference of two z scores.
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> r.test(30,.4,.6)
%\end{Sinput}
%\begin{Soutput}
%Correlation tests
%Call:r.test(n = 30, r12 = 0.4, r34 = 0.6)
%Test of difference between two independent correlations
% z value 0.99 with probability 0.32
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%
%3) For sample size n, and correlations ra= r12, rb= r23 and r13 specified, test for the difference of two dependent correlations (Steiger case A).
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> r.test(103,.4,.5,.1)
%\end{Sinput}
%\begin{Soutput}
%Correlation tests
%Call:r.test(n = 103, r12 = 0.4, r34 = 0.5, r23 = 0.1)
%Test of difference between two correlated correlations
% t value -0.89 with probability < 0.37
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%4) For sample size n, test for the difference between two dependent correlations involving different variables. (Steiger case B).
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> r.test(103,.5,.6,.7,.5,.5,.8) #steiger Case B
%\end{Sinput}
%\begin{Soutput}
%Correlation tests
%Call:r.test(n = 103, r12 = 0.5, r34 = 0.6, r23 = 0.7, r13 = 0.5, r14 = 0.5,
% r24 = 0.8)
%Test of difference between two dependent correlations
% z value -1.2 with probability 0.23
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%
%To test whether a matrix of correlations differs from what would be expected if the population correlations were all zero, the function \pfun{cortest} follows \cite{steiger:80b} who pointed out that the sum of the squared elements of a correlation matrix, or the Fisher z score equivalents, is distributed as chi square under the null hypothesis that the values are zero (i.e., elements of the identity matrix). This is particularly useful for examining whether correlations in a single matrix differ from zero or for comparing two matrices. Although obvious, \pfun{cortest} can be used to test whether the \pfun{sat.act} data matrix produces non-zero correlations (it does). This is a much more appropriate test when testing whether a residual matrix differs from zero.
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> cortest(sat.act)
%\end{Sinput}
%\begin{Soutput}
%Tests of correlation matrices
%Call:cortest(R1 = sat.act)
% Chi Square value 1325.42 with df = 15 with probability < 1.8e-273
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%\subsection{Polychoric, tetrachoric, polyserial, and biserial correlations}
%
%The Pearson correlation of dichotomous data is also known as the $\phi$ coefficient. If the data, e.g., ability items, are thought to represent an underlying continuous although latent variable, the $\phi$ will underestimate the value of the Pearson applied to these latent variables. One solution to this problem is to use the \pfun{tetrachoric} correlation which is based upon the assumption of a bivariate normal distribution that has been cut at certain points. The \pfun{draw.tetra} function demonstrates the process (Figure~\ref{fig:tetra}). A simple generalization of this to the case of the multiple cuts is the \pfun{polychoric} correlation.
%
%\begin{figure}[htbp]
%\begin{center}
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> draw.tetra()
%\end{Sinput}
%\end{Schunk}
%\includegraphics{factor-023}
%\end{scriptsize}
%\caption{The tetrachoric correlation estimates what a Pearson correlation would be given a two by two table of observed values assumed to be sampled from a bivariate normal distribution. The $\phi$ correlation is just a Pearson r performed on the observed values.}
%\label{fig:tetra}
%\end{center}
%\end{figure}
%
%Other estimated correlations based upon the assumption of bivariate normality with cut points include the \pfun{biserial} and \pfun{polyserial} correlation.
%
%If the data are a mix of continuous, polytomous and dichotomous variables, the \pfun{mixed.cor} function will calculate the appropriate mixture of Pearson, polychoric, tetrachoric, biserial, and polyserial correlations.
%
%The correlation matrix resulting from a number of tetrachoric or polychoric correlation matrix sometimes will not be positive semi-definite. This will also happen if the correlation matrix is formed by using pair-wise deletion of cases. The \pfun{cor.smooth} function will adjust the smallest eigen values of the correlation matrix to make them positive, rescale all of them to sum to the number of variables, and produce a ``smoothed'' correlation matrix. An example of this problem is a data set of \pfun{burt} which probably had a typo in the original correlation matrix. Smoothing the matrix corrects this problem.
%
%\subsection{Multiple regression from data or correlation matrices}
%
%The typical application of the \fun{lm} function is to do a linear model of one Y variable as a function of multiple X variables. Because \fun{lm} is designed to analyze complex interactions, it requires raw data as input. It is, however, sometimes convenient to do \iemph{multiple regression} from a correlation or covariance matrix. The \pfun{set.cor} function will do this, taking a set of y variables predicted from a set of x variables, perhaps with a set of z covariates removed from both x and y. Consider the \iemph{Thurstone} correlation matrix and find the multiple correlation of the last five variables as a function of the first 4.
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> set.cor(y = 5:9,x=1:4,data=Thurstone)
%\end{Sinput}
%\begin{Soutput}
%Call: set.cor(y = 5:9, x = 1:4, data = Thurstone)
%
%Multiple Regression from matrix input
%
%Beta weights
% 4.Letter.Words Suffixes Letter.Series Pedigrees Letter.Group
%Sentences 0.09 0.07 0.25 0.21 0.20
%Vocabulary 0.09 0.17 0.09 0.16 -0.02
%Sent.Completion 0.02 0.05 0.04 0.21 0.08
%First.Letters 0.58 0.45 0.21 0.08 0.31
%
%Multiple R
%4.Letter.Words Suffixes Letter.Series Pedigrees Letter.Group
% 0.69 0.63 0.50 0.58 0.48
%
%Multiple R2
%4.Letter.Words Suffixes Letter.Series Pedigrees Letter.Group
% 0.48 0.40 0.25 0.34 0.23
%
%Various estimates of between set correlations
%Squared Canonical Correlations
%[1] 0.6280 0.1478 0.0076 0.0049
%Chisq of canonical correlations
%NULL
%
% Average squared canonical correlation = 0.2
% Cohen's Set Correlation R2 = 0.69
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
%By specifying the number of subjects in correlation matrix, appropriate estimates of standard errors, t-values, and probabilities are also found. The next example finds the regressions with variables 1 and 2 used as covariates. The $\hat{\beta}$ weights for variables 3 and 4 do not change, but the multiple correlation is much less. It also shows how to find the residual correlations between variables 5-9 with variables 1-4 removed.
%
%\begin{scriptsize}
%\begin{Schunk}
%\begin{Sinput}
%> sc <- set.cor(y = 5:9,x=3:4,data=Thurstone,z=1:2)
%\end{Sinput}
%\begin{Soutput}
%Call: set.cor(y = 5:9, x = 3:4, data = Thurstone, z = 1:2)
%
%Multiple Regression from matrix input
%
%Beta weights
% 4.Letter.Words Suffixes Letter.Series Pedigrees Letter.Group
%Sent.Completion 0.02 0.05 0.04 0.21 0.08
%First.Letters 0.58 0.45 0.21 0.08 0.31
%
%Multiple R
%4.Letter.Words Suffixes Letter.Series Pedigrees Letter.Group
% 0.51 0.40 0.19 0.15 0.28
%
%Multiple R2
%4.Letter.Words Suffixes Letter.Series Pedigrees Letter.Group
% 0.258 0.158 0.034 0.022 0.078
%
%Various estimates of between set correlations
%Squared Canonical Correlations
%[1] 0.405 0.023
%Chisq of canonical correlations
%NULL
%
% Average squared canonical correlation = 0.21
% Cohen's Set Correlation R2 = 0.42
%\end{Soutput}
%\begin{Sinput}
%> round(sc$residual,2)
%\end{Sinput}
%\begin{Soutput}
% 4.Letter.Words Suffixes Letter.Series Pedigrees Letter.Group
%4.Letter.Words 0.52 0.11 0.09 0.06 0.13
%Suffixes 0.11 0.60 -0.01 0.01 0.03
%Letter.Series 0.09 -0.01 0.75 0.28 0.37
%Pedigrees 0.06 0.01 0.28 0.66 0.20
%Letter.Group 0.13 0.03 0.37 0.20 0.77
%\end{Soutput}
%\end{Schunk}
%\end{scriptsize}
%
\section{Item and scale analysis}
The main functions in the \Rpkg{psych} package are for analyzing the structure of items and of scales and for finding various estimates of scale reliability. These may be considered as problems of dimension reduction (e.g., factor analysis, cluster analysis, principal components analysis) and of forming and estimating the reliability of the resulting composite scales.
\subsection{Dimension reduction through factor analysis and cluster analysis}
\label{sect:fa}
Parsimony of description has been a goal of science since at least the famous dictum commonly attributed to William of Ockham to not multiply entities beyond necessity\footnote{Although probably neither original with Ockham nor directly stated by him \citep{thornburn:1918}, Ockham's razor remains a fundamental principal of science.}. The goal for parsimony is seen in psychometrics as an attempt either to describe (components) or to explain (factors) the relationships between many observed variables in terms of a more limited set of components or latent factors.
The typical data matrix represents multiple items or scales usually thought to reflect fewer underlying constructs\footnote{\cite{cattell:fa78} as well as \cite{maccallum:07} argue that the data are the result of many more factors than observed variables, but are willing to estimate the major underlying factors.}. At the most simple, a set of items can be be thought to represent a random sample from one underlying domain or perhaps a small set of domains. The question for the psychometrician is how many domains are represented and how well does each item represent the domains. Solutions to this problem are examples of \iemph{factor analysis} (\iemph{FA}), \iemph{principal components analysis} (\iemph{PCA}), and \iemph{cluster analysis} (\emph{CA}). All of these procedures aim to reduce the complexity of the observed data. In the case of FA, the goal is to identify fewer underlying constructs to explain the observed data. In the case of PCA, the goal can be mere data reduction, but the interpretation of components is frequently done in terms similar to those used when describing the latent variables estimated by FA. Cluster analytic techniques, although usually used to partition the subject space rather than the variable space, can also be used to group variables to reduce the complexity of the data by forming fewer and more homogeneous sets of tests or items.
At the data level the data reduction problem may be solved as a \iemph{Singular Value Decomposition} of the original matrix, although the more typical solution is to find either the \iemph{principal components} or \iemph{factors} of the covariance or correlation matrices. Given the pattern of regression weights from the variables to the components or from the factors to the variables, it is then possible to find (for components) individual \index{component scores} \emph{component} or \iemph{cluster scores} or estimate (for factors) \iemph{factor scores}.
Several of the functions in \Rpkg{psych} address the problem of data reduction.
\begin{description}
\item[\pfun{fa}] incorporates five alternative algorithms: \iemph{minres factor analysis}, \iemph{principal axis factor analysis}, \iemph{weighted least squares factor analysis}, \iemph{generalized least squares factor analysis} and \iemph{maximum likelihood factor analysis}. That is, it includes the functionality of three other functions that will be eventually phased out.
\item [\pfun{factor.minres}] Minimum residual factor analysis is a least squares, iterative solution to the factor problem. minres attempts to minimize the residual (off-diagonal) correlation matrix. It produces solutions similar to maximum likelihood solutions, but will work even if the matrix is singular.
\item [\pfun{factor.pa}] Principal Axis factor analysis is a least squares, iterative solution to the factor problem. PA will work for cases where maximum likelihood techniques (\fun{factanal}) will not work. The original communality estimates are either the squared multiple correlations (\pfun{smc}) for each item or 1.
\item [\pfun{factor.wls}] Weighted least squares factor analysis is a least squares, iterative solution to the factor problem. It minimizes the (weighted) squared residual matrix. The weights are based upon the independent contribution of each variable.
\item [\pfun{principal}] Principal Components Analysis reports the largest n eigen vectors rescaled by the square root of their eigen values.
\item [\pfun{factor.congruence}] The congruence between two factors is the cosine of the angle between them. This is just the cross products of the loadings divided by the sum of the squared loadings. This differs from the correlation coefficient in that the mean loading is not subtracted before taking the products. \pfun{factor.congruence} will find the cosines between two (or more) sets of factor loadings.
\item [\pfun{vss}] Very Simple Structure \cite{revelle:vss} applies a goodness of fit test to determine the optimal number of factors to extract. It can be thought of as a quasi-confirmatory model, in that it fits the very simple structure (all except the biggest c loadings per item are set to zero where c is the level of complexity of the item) of a factor pattern matrix to the original correlation matrix. For items where the model is usually of complexity one, this is equivalent to making all except the largest loading for each item 0. This is typically the solution that the user wants to interpret. The analysis includes the \pfun{MAP} criterion of \cite{velicer:76} and a $\chi^2$ estimate.
\item [\pfun{fa.parallel}] The parallel factors technique compares the observed eigen values of a correlation matrix with those from random data.
\item [\pfun{fa.plot}] will plot the loadings from a factor, principal components, or cluster analysis (just a call to plot will suffice). If there are more than two factors, then a SPLOM of the loadings is generated.
\item[\pfun{fa.diagram}] replaces \pfun{fa.graph} and will draw a path diagram representing the factor structure. It does not require Rgraphviz and thus is probably preferred.
\item[\pfun{fa.graph}] requires \fun{Rgraphviz} and will draw a graphic representation of the factor structure. If factors are correlated, this will be represented as well.
\item[\pfun{iclust} ] is meant to do item cluster analysis using a hierarchical clustering algorithm specifically asking questions about the reliability of the clusters \citep{revelle:iclust}. Clusters are formed until either coefficient $\alpha$ \cite{cronbach:51} or $\beta$ \cite{revelle:iclust} fail to increase.
\end{description}
\subsubsection{Minimum Residual Factor Analysis}
\label{sect:minres}
The factor model is an approximation of a correlation matrix by a matrix of lower rank. That is, can the correlation matrix, $\vec{_nR_n}$ be approximated by the product of a factor matrix, $\vec{_nF_k}$ and its transpose plus a diagonal matrix of uniqueness.
\begin{equation}
R = FF' + U^2
\end{equation}
The maximum likelihood solution to this equation is found by \fun{factanal} in the \Rpkg{stats} package. Five alternatives are provided in \Rpkg{psych}, all of them are included in the \pfun{fa} function and are called by specifying the factor method (e.g., fm=``minres", fm=``pa", fm=``"wls", fm="gls" and fm="ml"). In the discussion of the other algorithms, the calls shown are to the \pfun{fa} function specifying the appropriate method.
\pfun{factor.minres} attempts to minimize the off diagonal residual correlation matrix by adjusting the eigen values of the original correlation matrix. This is similar to what is done in \fun{factanal}, but uses an ordinary least squares instead of a maximum likelihood fit function. The solutions tend to be more similar to the MLE solutions than are the \pfun{factor.pa} solutions. \iemph{min.res} is the default for the \pfun{fa} function.
A classic data set, collected by \cite{thurstone:41} and then reanalyzed by \cite{bechtoldt:61} and discussed by \cite{mcdonald:tt}, is a set of 9 cognitive variables with a clear bi-factor structure \cite{holzinger:37}. The minimum residual solution was transformed into an oblique solution using the default option on rotate which uses an oblimin transformation (Table~\ref{tab:factor.minres}). Alternative rotations and transformations include ``none", ``varimax", ``quartimax", ``bentlerT", and ``geominT" (all of which are orthogonal rotations). as well as ``promax", ``oblimin", ``simplimax", ``bentlerQ, and``geominQ" and ``cluster" which are possible oblique transformations of the solution. The default is to do a oblimin transformation, although prior versions defaulted to varimax. The measures of factor adequacy reflect the multiple correlations of the factors with the best fitting linear regression estimates of the factor scores \citep{grice:01}.
\begin{table}[htdp]
\caption{Three correlated factors from the Thurstone 9 variable problem. By default, the solution is transformed obliquely using oblimin. The extraction method is (by default) minimum residual.}
\begin{scriptsize}
\begin{center}
\begin{Schunk}
\begin{Sinput}
> f3t <- fa(Thurstone,3,n.obs=213)
> f3t
\end{Sinput}
\begin{Soutput}
Factor Analysis using method = minres
Call: fa(r = Thurstone, nfactors = 3, n.obs = 213)
Standardized loadings (pattern matrix) based upon correlation matrix
MR1 MR2 MR3 h2 u2
Sentences 0.91 -0.04 0.04 0.82 0.18
Vocabulary 0.89 0.06 -0.03 0.84 0.16
Sent.Completion 0.83 0.04 0.00 0.73 0.27
First.Letters 0.00 0.86 0.00 0.73 0.27
4.Letter.Words -0.01 0.74 0.10 0.63 0.37
Suffixes 0.18 0.63 -0.08 0.50 0.50
Letter.Series 0.03 -0.01 0.84 0.72 0.28
Pedigrees 0.37 -0.05 0.47 0.50 0.50
Letter.Group -0.06 0.21 0.64 0.53 0.47
MR1 MR2 MR3
SS loadings 2.64 1.86 1.50
Proportion Var 0.29 0.21 0.17
Cumulative Var 0.29 0.50 0.67
Proportion Explained 0.44 0.31 0.25
Cumulative Proportion 0.44 0.75 1.00
With factor correlations of
MR1 MR2 MR3
MR1 1.00 0.59 0.54
MR2 0.59 1.00 0.52
MR3 0.54 0.52 1.00
Test of the hypothesis that 3 factors are sufficient.
The degrees of freedom for the null model are 36 and the objective function was 5.2 with Chi Square of 1081.97
The degrees of freedom for the model are 12 and the objective function was 0.01
The root mean square of the residuals (RMSR) is 0.01
The df corrected root mean square of the residuals is 0.02
The harmonic number of observations is 213 with the empirical chi square 0.58 with prob < 1
The total number of observations was 213 with MLE Chi Square = 2.82 with prob < 1
Tucker Lewis Index of factoring reliability = 1.027
RMSEA index = 0 and the 90 % confidence intervals are NA NA
BIC = -61.51
Fit based upon off diagonal values = 1
Measures of factor score adequacy
MR1 MR2 MR3
Correlation of scores with factors 0.96 0.92 0.90
Multiple R square of scores with factors 0.93 0.85 0.81
Minimum correlation of possible factor scores 0.86 0.71 0.63
\end{Soutput}
\end{Schunk}
\end{center}
\end{scriptsize}
\label{tab:factor.minres}
\end{table}%
\subsubsection{Principal Axis Factor Analysis}
An alternative, least squares algorithm, \pfun{factor.pa}, does a Principal Axis factor analysis by iteratively doing an eigen value decomposition of the correlation matrix with the diagonal replaced by the values estimated by the factors of the previous iteration. This OLS solution is not as sensitive to improper matrices as is the maximum likelihood method, and will sometimes produce more interpretable results. It seems as if the SAS example for PA uses only one iteration. Setting the max.iter parameter to 1 produces the SAS solution.
The solutions from the \pfun{fa}, the \pfun{factor.minres} and \pfun{factor.pa} as well as the \pfun{principal} functions can be rotated or transformed with a number of options. Some of these call the \Rpkg{GPArotation} package. Orthogonal rotations are \fun{varimax} and \fun{quartimax}. Oblique transformations include \fun{oblimin}, \fun{quartimin} and then two targeted rotation functions \pfun{Promax} and \pfun{target.rot}. The latter of these will transform a loadings matrix towards an arbitrary target matrix. The default is to transform towards an independent cluster solution.
Using the Thurstone data set, three factors were requested and then transformed into an independent clusters solution using \pfun{target.rot} (Table~\ref{tab:Thurstone}).
\begin{table}[htdp]
\caption{The 9 variable problem from Thurstone is a classic example of factoring where there is a higher order factor, g, that accounts for the correlation between the factors. The extraction method was principal axis. The transformation was a targeted transformation to a simple cluster solution.}
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> f3 <- fa(Thurstone,3,n.obs = 213,fm="pa")
> f3o <- target.rot(f3)
> f3o
\end{Sinput}
\begin{Soutput}
Call: NULL
Standardized loadings (pattern matrix) based upon correlation matrix
PA1 PA2 PA3 h2 u2
Sentences 0.89 -0.03 0.07 0.81 0.19
Vocabulary 0.89 0.07 0.00 0.80 0.20
Sent.Completion 0.83 0.04 0.03 0.70 0.30
First.Letters -0.02 0.85 -0.01 0.73 0.27
4.Letter.Words -0.05 0.74 0.09 0.57 0.43
Suffixes 0.17 0.63 -0.09 0.43 0.57
Letter.Series -0.06 -0.08 0.84 0.69 0.31
Pedigrees 0.33 -0.09 0.48 0.37 0.63
Letter.Group -0.14 0.16 0.64 0.45 0.55
PA1 PA2 PA3
SS loadings 2.45 1.72 1.37
Proportion Var 0.27 0.19 0.15
Cumulative Var 0.27 0.46 0.62
Proportion Explained 0.44 0.31 0.25
Cumulative Proportion 0.44 0.75 1.00
PA1 PA2 PA3
PA1 1.00 0.02 0.08
PA2 0.02 1.00 0.09
PA3 0.08 0.09 1.00
\end{Soutput}
\end{Schunk}
\end{scriptsize}
\end{center}
\label{tab:Thurstone}
\end{table}
\subsubsection{Weighted Least Squares Factor Analysis}
\label{sect:wls}
Similar to the minres approach of minimizing the squared residuals, factor method ``wls" weights the squared residuals by their uniquenesses. This tends to produce slightly smaller overall residuals. In the example of weighted least squares, the output is shown by using the \pfun{print} function with the cut option set to 0. That is, all loadings are shown (Table~\ref{tab:Thurstone.wls}).
\begin{table}[htdp]
\caption{The 9 variable problem from Thurstone is a classic example of factoring where there is a higher order factor, g, that accounts for the correlation between the factors. The factors were extracted using a weighted least squares algorithm. All loadings are shown by using the cut=0 option in the \pfun{print.psych} function.}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> f3w <- fa(Thurstone,3,n.obs = 213,fm="wls")
> print(f3w,cut=0,digits=3)
\end{Sinput}
\begin{Soutput}
Factor Analysis using method = wls
Call: fa(r = Thurstone, nfactors = 3, n.obs = 213, fm = "wls")
Standardized loadings (pattern matrix) based upon correlation matrix
WLS1 WLS2 WLS3 h2 u2
Sentences 0.905 -0.034 0.040 0.822 0.178
Vocabulary 0.890 0.066 -0.031 0.835 0.165
Sent.Completion 0.833 0.034 0.007 0.735 0.265
First.Letters -0.002 0.855 0.003 0.731 0.269
4.Letter.Words -0.016 0.743 0.106 0.629 0.371
Suffixes 0.180 0.626 -0.082 0.496 0.504
Letter.Series 0.033 -0.015 0.838 0.719 0.281
Pedigrees 0.381 -0.051 0.464 0.505 0.495
Letter.Group -0.062 0.209 0.632 0.527 0.473
WLS1 WLS2 WLS3
SS loadings 2.647 1.864 1.488
Proportion Var 0.294 0.207 0.165
Cumulative Var 0.294 0.501 0.667
Proportion Explained 0.441 0.311 0.248
Cumulative Proportion 0.441 0.752 1.000
With factor correlations of
WLS1 WLS2 WLS3
WLS1 1.000 0.591 0.535
WLS2 0.591 1.000 0.516
WLS3 0.535 0.516 1.000
Test of the hypothesis that 3 factors are sufficient.
The degrees of freedom for the null model are 36 and the objective function was 5.198 with Chi Square of 1081.968
The degrees of freedom for the model are 12 and the objective function was 0.014
The root mean square of the residuals (RMSR) is 0.006
The df corrected root mean square of the residuals is 0.014
The harmonic number of observations is 213 with the empirical chi square 0.531 with prob < 1
The total number of observations was 213 with MLE Chi Square = 2.886 with prob < 0.996
Tucker Lewis Index of factoring reliability = 1.0264
RMSEA index = 0 and the 90 % confidence intervals are NA NA
BIC = -61.45
Fit based upon off diagonal values = 1
Measures of factor score adequacy
WLS1 WLS2 WLS3
Correlation of scores with factors 0.964 0.923 0.902
Multiple R square of scores with factors 0.929 0.853 0.814
Minimum correlation of possible factor scores 0.858 0.706 0.627
\end{Soutput}
\end{Schunk}
\end{scriptsize}
\label{tab:Thurstone.wls}
\end{table}
The unweighted least squares solution may be shown graphically using the \pfun{fa.plot} function which is called by the generic \fun{plot} function (Figure~\ref{fig:thurstone}. Factors were transformed obliquely using a oblimin. These solutions may be shown as item by factor plots (Figure~\ref{fig:thurstone} or by a structure diagram (Figure~\ref{fig:thurstone.diagram}.
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> plot(f3t)
\end{Sinput}
\end{Schunk}
\includegraphics{factor-029}
\end{scriptsize}
\caption{A graphic representation of the 3 oblique factors from the Thurstone data using \pfun{plot}. Factors were transformed to an oblique solution using the oblimin function from the GPArotation package.}
\label{fig:thurstone}
\end{center}
\end{figure}
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> fa.diagram(f3t)
\end{Sinput}
\end{Schunk}
\includegraphics{factor-030}
\end{scriptsize}
\caption{A graphic representation of the 3 oblique factors from the Thurstone data using \pfun{fa.diagram}. Factors were transformed to an oblique solution using oblimin.}
\label{fig:thurstone.diagram}
\end{center}
\end{figure}
A comparison of these three approaches suggests that the minres solution is more similar to a maximum likelihood solution and fits slightly better than the pa or wls solutions. Comparisons with SPSS suggest that the pa solution matches the SPSS OLS solution, but that the minres solution is slightly better. At least in one test data set, the weighted least squares solutions, although fitting equally well, had slightly different structure loadings. Note that the rotations used by SPSS will sometimes use the ``Kaiser Normalization''. By default, the rotations used in psych do not normalize, but this can be specified as an option in \pfun{fa}.
\subsubsection{Principal Components analysis (PCA)}
An alternative to factor analysis, which is unfortunately frequently confused with \iemph{factor analysis}, is \iemph{principal components analysis}. Although the goals of \iemph{PCA} and \iemph{FA} are similar, PCA is a descriptive model of the data, while FA is a structural model. Psychologists typically use PCA in a manner similar to factor analysis and thus the \pfun{principal} function produces output that is perhaps more understandable than that produced by \fun{princomp} in the \Rpkg{stats} package. Table~\ref{tab:pca} shows a PCA of the Thurstone 9 variable problem rotated using the \pfun{Promax} function. Note how the loadings from the factor model are similar but smaller than the principal component loadings. This is because the PCA model attempts to account for the entire variance of the correlation matrix, while FA accounts for just the \iemph{common variance}. This distinction becomes most important for small correlation matrices. Also note how the goodness of fit statistics, based upon the residual off diagonal elements, is much worse than the \pfun{fa} solution.
\begin{table}[htdp]
\caption{The Thurstone problem can also be analyzed using Principal Components Analysis. Compare this to Table~\ref{tab:Thurstone}. The loadings are higher for the PCA because the model accounts for the unique as well as the common variance.The fit of the off diagonal elements, however, is much worse than the \pfun{fa} results.}
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> p3p <-principal(Thurstone,3,n.obs = 213,rotate="Promax")
> p3p
\end{Sinput}
\begin{Soutput}
Principal Components Analysis
Call: principal(r = Thurstone, nfactors = 3, rotate = "Promax", n.obs = 213)
Standardized loadings (pattern matrix) based upon correlation matrix
PC1 PC2 PC3 h2 u2
Sentences 0.92 0.01 0.01 0.86 0.14
Vocabulary 0.90 0.10 -0.05 0.86 0.14
Sent.Completion 0.91 0.04 -0.04 0.83 0.17
First.Letters 0.01 0.84 0.07 0.78 0.22
4.Letter.Words -0.05 0.81 0.17 0.75 0.25
Suffixes 0.18 0.79 -0.15 0.70 0.30
Letter.Series 0.03 -0.03 0.88 0.78 0.22
Pedigrees 0.45 -0.16 0.57 0.67 0.33
Letter.Group -0.19 0.19 0.86 0.75 0.25
PC1 PC2 PC3
SS loadings 2.83 2.19 1.96
Proportion Var 0.31 0.24 0.22
Cumulative Var 0.31 0.56 0.78
Proportion Explained 0.41 0.31 0.28
Cumulative Proportion 0.41 0.72 1.00
With component correlations of
PC1 PC2 PC3
PC1 1.00 0.51 0.53
PC2 0.51 1.00 0.44
PC3 0.53 0.44 1.00
Test of the hypothesis that 3 components are sufficient.
The degrees of freedom for the null model are 36 and the objective function was 5.2
The degrees of freedom for the model are 12 and the objective function was 0.62
The total number of observations was 213 with MLE Chi Square = 127.9 with prob < 1.6e-21
Fit based upon off diagonal values = 0.98
\end{Soutput}
\end{Schunk}
\end{scriptsize}
\end{center}
\label{tab:pca}
\end{table}
\subsubsection{Hierarchical and bi-factor solutions}
\label{sect:omega}
For a long time structural analysis of the ability domain have considered the problem of factors that are themselves correlated. These correlations may themselves be factored to produce a higher order, general factor. An alternative \citep{holzinger:37,jensen:weng} is to consider the general factor affecting each item, and then to have group factors account for the residual variance. Exploratory factor solutions to produce a hierarchical or a bifactor solution are found using the \pfun{omega} function. This technique has more recently been applied to the personality domain to consider such things as the structure of neuroticism (treated as a general factor, with lower order factors of anxiety, depression, and aggression).
Consider the 9 Thurstone variables analyzed in the prior factor analyses. The correlations between the factors (as shown in Figure~\ref{fig:thurstone.diagram} can themselves be factored. This results in a higher order factor model (Figure~\ref{fig:omega}). An an alternative solution is to take this higher order model and then solve for the general factor loadings as well as the loadings on the residualized lower order factors using the \iemph{Schmid-Leiman} procedure. (Figure ~\ref{fig:omega.2}). Yet another solution is to use structural equation modeling to directly solve for the general and group factors.
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> om.h <- omega(Thurstone,n.obs=213,sl=FALSE)
> op <- par(mfrow=c(1,1))
\end{Sinput}
\end{Schunk}
\includegraphics{factor-032}
\end{scriptsize}
\caption{A higher order factor solution to the Thurstone 9 variable problem}
\label{fig:omega}
\end{center}
\end{figure}
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> om <- omega(Thurstone,n.obs=213)
\end{Sinput}
\end{Schunk}
\includegraphics{factor-033}
\end{scriptsize}
\caption{A bifactor factor solution to the Thurstone 9 variable problem}
\label{fig:omega.2}
\end{center}
\end{figure}
Yet another approach to the bifactor structure is do use the \pfun{bifactor} rotation function in either \Rpkg{psych} or in \Rpkg{GPArotation}. This does the rotation discussed in \cite{jennrich:11}.
\subsubsection{Item Cluster Analysis: iclust}
\label{sect:iclust}
An alternative to factor or components analysis is \iemph{cluster analysis}. The goal of cluster analysis is the same as factor or components analysis (reduce the complexity of the data and attempt to identify homogeneous subgroupings). Mainly used for clustering people or objects (e.g., projectile points if an anthropologist, DNA if a biologist, galaxies if an astronomer), clustering may be used for clustering items or tests as well. Introduced to psychologists by \cite{tryon:39} in the 1930's, the cluster analytic literature exploded in the 1970s and 1980s \citep{blashfield:80,blashfield:88,everitt:74,hartigan:75}. Much of the research is in taxonmetric applications in biology \citep{sneath:73,sokal:63} and marketing \citep{cooksey:06} where clustering remains very popular. It is also used for taxonomic work in forming clusters of people in family \citep{henry:05} and clinical psychology \citep{martinent:07,mun:08}. Interestingly enough it has has had limited applications to psychometrics. This is unfortunate, for as has been pointed out by e.g. \citep{tryon:35,loevinger:53}, the theory of factors, while mathematically compelling, offers little that the geneticist or behaviorist or perhaps even non-specialist finds compelling. \cite{cooksey:06} reviews why the \pfun{iclust} algorithm is particularly appropriate for scale construction in marketing.
\emph{Hierarchical cluster analysis} \index{hierarchical cluster analysis} forms clusters that are nested within clusters. The resulting \iemph{tree diagram} (also known somewhat pretentiously as a \iemph{rooted dendritic structure}) shows the nesting structure. Although there are many hierarchical clustering algorithms in \R{} (e.g., \fun{agnes}, \fun{hclust}, and \pfun{iclust}), the one most applicable to the problems of scale construction is \pfun{iclust} \citep{revelle:iclust}.
\begin{enumerate}
\item Find the proximity (e.g. correlation) matrix,
\item Identify the most similar pair of items
\item Combine this most similar pair of items to form a new variable (cluster),
\item Find the similarity of this cluster to all other items and clusters,
\item Repeat steps 2 and 3 until some criterion is reached (e.g., typicallly, if only one cluster remains or in \pfun{iclust} if there is a failure to increase reliability coefficients $\alpha$ or $\beta$).
\item Purify the solution by reassigning items to the most similar cluster center.
\end{enumerate}
\pfun{iclust} forms clusters of items using a hierarchical clustering algorithm until one of two measures of internal consistency fails to increase \citep{revelle:iclust}. The number of clusters may be specified a priori, or found empirically. The resulting statistics include the average split half reliability, $\alpha$ \citep{cronbach:51}, as well as the worst split half reliability, $\beta$ \citep{revelle:iclust}, which is an estimate of the general factor saturation of the resulting scale (Figure~\ref{fig:iclust}). Cluster loadings (corresponding to the structure matrix of factor analysis) are reported when printing (Table~\ref{tab:iclust}). The pattern matrix is available as an object in the results.
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> data(bfi)
> ic <- iclust(bfi[1:25])
\end{Sinput}
\end{Schunk}
\includegraphics{factor-034}
\end{scriptsize}
\caption{Using the \pfun{iclust} function to find the cluster structure of 25 personality items (the three demographic variables were excluded from this analysis). When analyzing many variables, the tree structure may be seen more clearly if the graphic output is saved as a pdf and then enlarged using a pdf viewer.}
\label{fig:iclust}
\end{center}
\end{figure}
\begin{table}[htdp]
\caption{The summary statistics from an iclust analysis shows three large clusters and smaller cluster.}
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> summary(ic) #show the results
\end{Sinput}
\begin{Soutput}
ICLUST (Item Cluster Analysis)Call: iclust(r.mat = bfi[1:25])
ICLUST
Purified Alpha:
C20 C16 C15 C21
0.80 0.81 0.73 0.61
Guttman Lambda6*
C20 C16 C15 C21
0.82 0.81 0.72 0.61
Original Beta:
C20 C16 C15 C21
0.63 0.76 0.67 0.27
Cluster size:
C20 C16 C15 C21
10 5 5 5
Purified scale intercorrelations
reliabilities on diagonal
correlations corrected for attenuation above diagonal:
C20 C16 C15 C21
C20 0.80 -0.291 0.40 -0.33
C16 -0.24 0.815 -0.29 0.11
C15 0.30 -0.221 0.73 -0.30
C21 -0.23 0.074 -0.20 0.61
\end{Soutput}
\end{Schunk}
\end{scriptsize}
\end{center}
\label{tab:iclust}
\end{table}%
The previous analysis (Figure~\ref{fig:iclust}) was done using the Pearson correlation. A somewhat cleaner structure is obtained when using the \pfun{polychoric} function to find polychoric correlations (Figure~\ref{fig:iclust.poly}). Note that the first time finding the polychoric correlations some time, but the next three analyses were done using that correlation matrix (r.poly\$rho). When using the console for input, \pfun{polychoric} will report on its progress while working using \pfun{progressBar}.
\begin{table}[htdp]
\caption{The \pfun{polychoric} and the \pfun{tetrachoric} functions can take a long time to finish and report their progress by a series of dots as they work. The dots are suppressed when creating a Sweave document.}
\begin{center}
\begin{tiny}
\begin{Schunk}
\begin{Sinput}
> data(bfi)
> r.poly <- polychoric(bfi[1:25]) #the ... indicate the progress of the function
\end{Sinput}
\begin{Soutput}
\end{Soutput}
\end{Schunk}
\end{tiny}
\end{center}
\label{tab:bad}
\end{table}%
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> ic.poly <- iclust(r.poly$rho,title="ICLUST using polychoric correlations")
> iclust.diagram(ic.poly)
\end{Sinput}
\end{Schunk}
\includegraphics{factor-037}
\end{scriptsize}
\caption{ICLUST of the BFI data set using polychoric correlations. Compare this solution to the previous one (Figure~\ref{fig:iclust}) which was done using Pearson correlations. }
\label{fig:iclust.poly}
\end{center}
\end{figure}
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> ic.poly <- iclust(r.poly$rho,5,title="ICLUST using polychoric correlations for nclusters=5")
> iclust.diagram(ic.poly)
\end{Sinput}
\end{Schunk}
\includegraphics{factor-038}
\end{scriptsize}
\caption{ICLUST of the BFI data set using polychoric correlations with the solution set to 5 clusters. Compare this solution to the previous one (Figure~\ref{fig:iclust.poly}) which was done without specifying the number of clusters and to the next one (Figure~\ref{fig:iclust.3}) which was done by changing the beta criterion. }
\label{fig:iclust.5}
\end{center}
\end{figure}
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> ic.poly <- iclust(r.poly$rho,beta.size=3,title="ICLUST beta.size=3")
\end{Sinput}
\end{Schunk}
\includegraphics{factor-039}
\end{scriptsize}
\caption{ICLUST of the BFI data set using polychoric correlations with the beta criterion set to 3. Compare this solution to the previous three (Figure~\ref{fig:iclust},~\ref{fig:iclust.poly}, \ref{fig:iclust.5}).}
\label{fig:iclust.3}
\end{center}
\end{figure}
\begin{table}[htdp]
\caption{The output from \pfun{iclust}includes the loadings of each item on each cluster. These are equivalent to factor structure loadings. By specifying the value of cut, small loadings are suppressed. The default is for cut=0.su }
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> print(ic,cut=.3)
\end{Sinput}
\begin{Soutput}
ICLUST (Item Cluster Analysis)
Call: iclust(r.mat = bfi[1:25])
Purified Alpha:
C20 C16 C15 C21
0.80 0.81 0.73 0.61
G6* reliability:
C20 C16 C15 C21
0.77 0.81 0.59 0.45
Original Beta:
C20 C16 C15 C21
0.63 0.76 0.67 0.27
Cluster size:
C20 C16 C15 C21
10 5 5 5
Item by Cluster Structure matrix:
O P C20 C16 C15 C21
A1 C20 C20
A2 C20 C20 0.59
A3 C20 C20 0.65
A4 C20 C20 0.43
A5 C20 C20 0.65
C1 C15 C15 0.54
C2 C15 C15 0.62
C3 C15 C15 0.54
C4 C15 C15 0.31 -0.66
C5 C15 C15 -0.30 0.36 -0.59
E1 C20 C20 -0.50
E2 C20 C20 -0.61 0.34
E3 C20 C20 0.59 -0.39
E4 C20 C20 0.66
E5 C20 C20 0.50 0.40 -0.32
N1 C16 C16 0.76
N2 C16 C16 0.75
N3 C16 C16 0.74
N4 C16 C16 -0.34 0.62
N5 C16 C16 0.55
O1 C20 C21 -0.53
O2 C21 C21 0.44
O3 C20 C21 0.39 -0.62
O4 C21 C21 -0.33
O5 C21 C21 0.53
With eigenvalues of:
C20 C16 C15 C21
3.2 2.6 1.9 1.5
Purified scale intercorrelations
reliabilities on diagonal
correlations corrected for attenuation above diagonal:
C20 C16 C15 C21
C20 0.80 -0.29 0.40 -0.33
C16 -0.24 0.81 -0.29 0.11
C15 0.30 -0.22 0.73 -0.30
C21 -0.23 0.07 -0.20 0.61
Cluster fit = 0.68 Pattern fit = 0.96 RMSR = 0.05
NULL
\end{Soutput}
\end{Schunk}
\end{scriptsize}
\end{center}
\label{tab:iclust}
\end{table}%
A comparison of these four cluster solutions suggests both a problem and an advantage of clustering techniques. The problem is that the solutions differ. The advantage is that the structure of the items may be seen more clearly when examining the clusters rather than a simple factor solution.
\subsection{Confidence intervals using bootstrapping techniques}
Exploratory factoring techniques are sometimes criticized because of the lack of statistical information on the solutions. Overall estimates of goodness of fit including $\chi^{2}$ and RMSEA are found in the \pfun{fa} and \pfun{omega} functions. Confidence intervals for the factor loadings may be found by doing multiple bootstrapped iterations of the original analysis. This is done by setting the n.iter parameter to the desired number of iterations. This can be done for factoring of Pearson correlation matrices as well as polychoric/tetrachoric matrices (See Table~\ref{tab:bootstrap}). Although the example value for the number of iterations is set to 20, more conventional analyses might use 1000 bootstraps. This will take much longer.
\begin{table}[htdp]
\caption{An example of bootstrapped confidence intervals on 10 items from the Big 5 inventory. The number of bootstrapped samples was set to 20. More conventional bootstrapping would use 100 or 1000 replications. }
\begin{tiny}
\begin{center}
\begin{Schunk}
\begin{Sinput}
> fa(bfi[1:10],2,n.iter=20)
\end{Sinput}
\begin{Soutput}
Factor Analysis with confidence intervals using method = minres
Call: fa(r = bfi[1:10], nfactors = 2, n.iter = 20)
Factor Analysis using method = minres
Call: fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate,
scores = scores, residuals = residuals, SMC = SMC, covar = covar,
missing = FALSE, impute = impute, min.err = min.err, max.iter = max.iter,
symmetric = symmetric, warnings = warnings, fm = fm, alpha = alpha,
oblique.scores = oblique.scores, np.obs = np.obs)
Standardized loadings (pattern matrix) based upon correlation matrix
MR2 MR1 h2 u2
A1 0.07 -0.40 0.15 0.85
A2 0.02 0.65 0.44 0.56
A3 -0.03 0.77 0.57 0.43
A4 0.15 0.44 0.26 0.74
A5 0.02 0.62 0.39 0.61
C1 0.57 -0.05 0.30 0.70
C2 0.62 -0.01 0.39 0.61
C3 0.54 0.03 0.30 0.70
C4 -0.66 0.01 0.43 0.57
C5 -0.57 -0.05 0.35 0.65
MR2 MR1
SS loadings 1.80 1.77
Proportion Var 0.18 0.18
Cumulative Var 0.18 0.36
Proportion Explained 0.50 0.50
Cumulative Proportion 0.50 1.00
With factor correlations of
MR2 MR1
MR2 1.00 0.32
MR1 0.32 1.00
Test of the hypothesis that 2 factors are sufficient.
The degrees of freedom for the null model are 45 and the objective function was 2.03 with Chi Square of 5664.89
The degrees of freedom for the model are 26 and the objective function was 0.17
The root mean square of the residuals (RMSR) is 0.04
The df corrected root mean square of the residuals is 0.07
The harmonic number of observations is 2762 with the empirical chi square 403.38 with prob < 2.6e-69
The total number of observations was 2800 with MLE Chi Square = 464.04 with prob < 9.2e-82
Tucker Lewis Index of factoring reliability = 0.865
RMSEA index = 0.078 and the 90 % confidence intervals are 0.071 0.084
BIC = 257.67
Fit based upon off diagonal values = 0.98
Measures of factor score adequacy
MR2 MR1
Correlation of scores with factors 0.86 0.88
Multiple R square of scores with factors 0.74 0.77
Minimum correlation of possible factor scores 0.49 0.54
Coefficients and bootstrapped confidence intervals
low MR2 upper low MR1 upper
A1 0.03 0.07 0.10 -0.44 -0.40 -0.36
A2 -0.03 0.02 0.05 0.61 0.65 0.70
A3 -0.05 -0.03 0.00 0.71 0.77 0.81
A4 0.10 0.15 0.19 0.40 0.44 0.49
A5 -0.02 0.02 0.06 0.57 0.62 0.66
C1 0.52 0.57 0.61 -0.09 -0.05 -0.02
C2 0.58 0.62 0.65 -0.05 -0.01 0.04
C3 0.48 0.54 0.59 -0.01 0.03 0.07
C4 -0.71 -0.66 -0.61 -0.03 0.01 0.03
C5 -0.63 -0.57 -0.53 -0.09 -0.05 -0.01
Interfactor correlations and bootstrapped confidence intervals
lower estimate upper
1 0.25 0.32 0.36
\end{Soutput}
\begin{Sinput}
>
\end{Sinput}
\end{Schunk}
\end{center}
\end{tiny}
\label{tab:bootstrap}
\end{table}%
\subsection{Comparing factor/component/cluster solutions}
Cluster analysis, factor analysis, and principal components analysis all produce structure matrices (matrices of correlations between the dimensions and the variables) that can in turn be compared in terms of the \iemph{congruence coefficient} which is just the cosine of the angle between the dimensions $$c_{f_{i}f_{j}} = \frac{\sum_{k=1}^{n}{f_{ik}f_{jk}}} {\sum{f_{ik}^{2}}\sum{f_{jk}^{2}}}.$$ Consider the case of a four factor solution and four cluster solution to the Big Five problem.
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> f4 <- fa(bfi[1:25],4,fm="pa")
> factor.congruence(f4,ic)
\end{Sinput}
\begin{Soutput}
C20 C16 C15 C21
PA1 0.92 -0.32 0.44 -0.40
PA2 -0.26 0.95 -0.33 0.12
PA3 0.35 -0.24 0.88 -0.37
PA4 0.29 -0.12 0.27 -0.90
\end{Soutput}
\end{Schunk}
\end{scriptsize}
A more complete comparison of oblique factor solutions (both minres and principal axis), bifactor and component solutions to the Thurstone data set is done using the \pfun{factor.congruence} function. (See table~\ref{tab:congruence}).
\begin{table}[htdp]
\caption{Congruence coefficients for oblique factor, bifactor and component solutions for the Thurstone problem.}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> factor.congruence(list(f3t,f3o,om,p3p))
\end{Sinput}
\begin{Soutput}
MR1 MR2 MR3 PA1 PA2 PA3 g F1* F2* F3* h2 PC1 PC2 PC3
MR1 1.00 0.06 0.09 1.00 0.06 0.13 0.72 1.00 0.06 0.09 0.74 1.00 0.08 0.04
MR2 0.06 1.00 0.08 0.03 1.00 0.06 0.60 0.06 1.00 0.08 0.57 0.04 0.99 0.12
MR3 0.09 0.08 1.00 0.01 0.01 1.00 0.52 0.09 0.08 1.00 0.51 0.06 0.02 0.99
PA1 1.00 0.03 0.01 1.00 0.04 0.05 0.67 1.00 0.03 0.01 0.69 1.00 0.06 -0.04
PA2 0.06 1.00 0.01 0.04 1.00 0.00 0.57 0.06 1.00 0.01 0.54 0.04 0.99 0.05
PA3 0.13 0.06 1.00 0.05 0.00 1.00 0.54 0.13 0.06 1.00 0.53 0.10 0.01 0.99
g 0.72 0.60 0.52 0.67 0.57 0.54 1.00 0.72 0.60 0.52 0.99 0.69 0.58 0.50
F1* 1.00 0.06 0.09 1.00 0.06 0.13 0.72 1.00 0.06 0.09 0.74 1.00 0.08 0.04
F2* 0.06 1.00 0.08 0.03 1.00 0.06 0.60 0.06 1.00 0.08 0.57 0.04 0.99 0.12
F3* 0.09 0.08 1.00 0.01 0.01 1.00 0.52 0.09 0.08 1.00 0.51 0.06 0.02 0.99
h2 0.74 0.57 0.51 0.69 0.54 0.53 0.99 0.74 0.57 0.51 1.00 0.71 0.56 0.49
PC1 1.00 0.04 0.06 1.00 0.04 0.10 0.69 1.00 0.04 0.06 0.71 1.00 0.06 0.00
PC2 0.08 0.99 0.02 0.06 0.99 0.01 0.58 0.08 0.99 0.02 0.56 0.06 1.00 0.05
PC3 0.04 0.12 0.99 -0.04 0.05 0.99 0.50 0.04 0.12 0.99 0.49 0.00 0.05 1.00
\end{Soutput}
\end{Schunk}
\end{scriptsize}
\label{tab:congruence}
\end{table}%
\subsection{Determining the number of dimensions to extract.}
How many dimensions to use to represent a correlation matrix is an unsolved problem in psychometrics. There are many solutions to this problem, none of which is uniformly the best. Henry Kaiser once said that ``a solution to the number-of factors problem in factor analysis is easy, that he used to make up one every morning before breakfast. But the problem, of course is to find \emph{the} solution, or at least a solution that others will regard quite highly not as the best" \cite{horn:79}.
Techniques most commonly used include
1) Extracting factors until the chi square of the residual matrix is not significant.
2) Extracting factors until the change in chi square from factor n to factor n+1 is not significant.
3) Extracting factors until the eigen values of the real data are less than the corresponding eigen values of a random data set of the same size (parallel analysis) \pfun{fa.parallel} \citep{horn:65}.
4) Plotting the magnitude of the successive eigen values and applying the scree test (a sudden drop in eigen values analogous to the change in slope seen when scrambling up the talus slope of a mountain and approaching the rock face \citep{cattell:scree}.
5) Extracting factors as long as they are interpretable.
6) Using the Very Structure Criterion (\pfun{vss}) \citep{revelle:vss}.
7) Using Wayne Velicer's Minimum Average Partial (\pfun{MAP}) criterion \citep{velicer:76}.
8) Extracting principal components until the eigen value < 1.
Each of the procedures has its advantages and disadvantages. Using either the chi square test or the change in square test is, of course, sensitive to the number of subjects and leads to the nonsensical condition that if one wants to find many factors, one simply runs more subjects. Parallel analysis is partially sensitive to sample size in that for large samples the eigen values of random factors will be very small. The scree test is quite appealing but can lead to differences of interpretation as to when the scree``breaks". Extracting interpretable factors means that the number of factors reflects the investigators creativity more than the data. vss, while very simple to understand, will not work very well if the data are very factorially complex. (Simulations suggests it will work fine if the complexities of some of the items are no more than 2). The eigen value of 1 rule, although the default for many programs, seems to be a rough way of dividing the number of variables by 3 and is probably the worst of all criteria.
An additional problem in determining the number of factors is what is considered a factor. Many treatments of factor analysis assume that the residual correlation matrix after the factors of interest are extracted is composed of just random error. An alternative concept is that the matrix is formed from major factors of interest but that there are also numerous minor factors of no substantive interest but that account for some of the shared covariance between variables. The presence of such minor factors can lead one to extract too many factors and to reject solutions on statistical grounds of misfit that are actually very good fits to the data. This problem is partially addressed later in the discussion of simulating complex structures using \pfun{sim.structure} and of small extraneous factors using the \pfun{sim.minor} function.
\subsubsection{Very Simple Structure}
\label{sect:vss}
The \pfun{vss} function compares the fit of a number of factor analyses with the loading matrix ``simplified" by deleting all except the c greatest loadings per item, where c is a measure of factor complexity \cite{revelle:vss}. Included in \pfun{vss} is the MAP criterion (Minimum Absolute Partial correlation) of \cite{velicer:76}.
Using the Very Simple Structure criterion for the bfi data suggests that 4 factors are optimal (Figure~\ref{fig:vss}). However, the MAP criterion suggests that 5 is optimal.
\begin{figure}[htbp]
\begin{center}
\begin{Schunk}
\begin{Sinput}
> vss <- vss(bfi[1:25],title="Very Simple Structure of a Big 5 inventory")
\end{Sinput}
\end{Schunk}
\includegraphics{factor-044}
\caption{The Very Simple Structure criterion for the number of factors compares solutions for various levels of item complexity and various numbers of factors. For the Big 5 Inventory, the complexity 1 and 2 solutions both achieve their maxima at four factors. This is in contrast to parallel analysis which suggests 6 and the MAP criterion which suggests 5. }
\label{fig:vss}
\end{center}
\end{figure}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> vss
\end{Sinput}
\begin{Soutput}
Very Simple Structure of Very Simple Structure of a Big 5 inventory
Call: VSS(x = x, n = n, rotate = rotate, diagonal = diagonal, fm = fm,
n.obs = n.obs, plot = plot, title = title)
VSS complexity 1 achieves a maximimum of 0.58 with 4 factors
VSS complexity 2 achieves a maximimum of 0.74 with 4 factors
The Velicer MAP criterion achieves a minimum of 0.01 with 5 factors
Velicer MAP
[1] 0.02 0.02 0.02 0.02 0.01 0.02 0.02 0.02
Very Simple Structure Complexity 1
[1] 0.49 0.54 0.57 0.58 0.53 0.54 0.52 0.52
Very Simple Structure Complexity 2
[1] 0.00 0.63 0.69 0.74 0.73 0.72 0.70 0.69
\end{Soutput}
\end{Schunk}
\end{scriptsize}
\subsubsection{Parallel Analysis}
\label{sect:fa.parallel}
An alternative way to determine the number of factors is to compare the solution to random data with the same properties as the real data set. If the input is a data matrix, the comparison includes random samples from the real data, as well as normally distributed random data with the same number of subjects and variables. For the BFI data, parallel analysis suggests that 6 factors might be most appropriate (Figure~\ref{fig:parallel}). It is interesting to compare \pfun{fa.parallel} with the \fun{paran} from the \Rpkg{paran} package. This latter uses smcs to estimate communalities. Simulations of known structures with a particular number of major factors but with the presence of trivial, minor (but not zero) factors, show that using smcs will tend to lead to too many factors.
\begin{figure}[htbp]
\begin{scriptsize}
\begin{center}
\begin{Schunk}
\begin{Sinput}
> fa.parallel(bfi[1:25],main="Parallel Analysis of a Big 5 inventory")
\end{Sinput}
\begin{Soutput}
Parallel analysis suggests that the number of factors = 6 and the number of components = 6
\end{Soutput}
\end{Schunk}
\includegraphics{factor-046}
\caption{Parallel analysis compares factor and principal components solutions to the real data as well as resampled data. Although vss suggests 4 factors, MAP 5, parallel analysis suggests 6. One more demonstration of Kaiser's dictum.}
\label{fig:parallel}
\end{center}
\end{scriptsize}
\end{figure}
A more tedious problem in terms of computation is to do parallel analysis of \iemph{polychoric} correlation matrices. This is done by \pfun{fa.parallel.poly}. By default the number of replications is 20. This is appropriate when choosing the number of factors from dicthotomous or polytomous data matrices.
\subsection{Factor extension}
Sometimes we are interested in the relationship of the factors in one space with the variables in a different space. One solution is to find factors in both spaces separately and then find the structural relationships between them. This is the technique of structural equation modeling in packages such as \Rpkg{sem} or \Rpkg{lavaan}. An alternative is to use the concept of \iemph{factor extension} developed by \citep{dwyer:37}. Consider the case of 16 variables created to represent one two dimensional space. If factors are found from eight of these variables, they may then be extended to the additional eight variables (See Figure~\ref{fig:fa.extension}).
\begin{figure}[htbp]
\begin{center}
\begin{scriptsize}
\begin{Schunk}
\begin{Sinput}
> v16 <- sim.item(16)
> s <- c(1,3,5,7,9,11,13,15)
> f2 <- fa(v16[,s],2)
> fe <- fa.extension(cor(v16)[s,-s],f2)
> fa.diagram(f2,fe=fe)
\end{Sinput}
\end{Schunk}
\includegraphics{factor-047}
\end{scriptsize}
\caption{Factor extension applies factors from one set (those on the left) to another set of variables (those on the right). \pfun{fa.extension} is particularly useful when one wants to define the factors with one set of variables and then apply those factors to another set. \pfun{fa.diagram} is used to show the structure. }
\label{fig:fa.extension}
\end{center}
\end{figure}
Another way to examine the overlap between two sets is the use of \iemph{set correlation} found by \pfun{set.cor} (discussed later).
\section{Psychometric Theory}
The \Rpkg{psych} package has been developed to help psychologists do basic research. Many of the functions were developed to supplement a book (\url{http://personality-project.org/r/book} An introduction to Psychometric Theory with Applications in \R{} \citep{revelle:intro} More information about the use of some of the functions may be found in the book .
For more extensive discussion of the use of \Rpkg{psych} in particular and \R{} in general, consult \url{http://personality-project.org/r/r.guide.html} A short guide to R.
\section{SessionInfo}
This document was prepared using the following settings.
\begin{tiny}
\begin{Schunk}
\begin{Sinput}
> sessionInfo()
\end{Sinput}
\begin{Soutput}
R Under development (unstable) (2013-01-13 r61644)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] sem_3.1-0 matrixcalc_1.0-3 mvtnorm_0.9-9994 GPArotation_2012.3-1
[5] MASS_7.3-23 psych_1.3.01
loaded via a namespace (and not attached):
[1] tools_3.0.0
\end{Soutput}
\end{Schunk}
\end{tiny}
\newpage
%\bibliography{/Volumes/WR/Documents/Active/book/all}
%\bibliography{all}
\begin{thebibliography}{}
\bibitem[\protect\astroncite{Bechtoldt}{1961}]{bechtoldt:61}
Bechtoldt, H. (1961).
\newblock An empirical study of the factor analysis stability hypothesis.
\newblock {\em Psychometrika}, 26(4):405--432.
\bibitem[\protect\astroncite{Blashfield}{1980}]{blashfield:80}
Blashfield, R.~K. (1980).
\newblock The growth of cluster analysis: {Tryon, Ward, and Johnson}.
\newblock {\em Multivariate Behavioral Research}, 15(4):439 -- 458.
\bibitem[\protect\astroncite{Blashfield and Aldenderfer}{1988}]{blashfield:88}
Blashfield, R.~K. and Aldenderfer, M.~S. (1988).
\newblock The methods and problems of cluster analysis.
\newblock In Nesselroade, J.~R. and Cattell, R.~B., editors, {\em Handbook of
multivariate experimental psychology (2nd ed.)}, pages 447--473. Plenum
Press, New York, NY.
\bibitem[\protect\astroncite{Bliese}{2009}]{bliese:09}
Bliese, P.~D. (2009).
\newblock Multilevel modeling in r (2.3) a brief introduction to r, the
multilevel package and the nlme package.
\bibitem[\protect\astroncite{Cattell}{1966}]{cattell:scree}
Cattell, R.~B. (1966).
\newblock The scree test for the number of factors.
\newblock {\em Multivariate Behavioral Research}, 1(2):245--276.
\bibitem[\protect\astroncite{Cattell}{1978}]{cattell:fa78}
Cattell, R.~B. (1978).
\newblock {\em The scientific use of factor analysis}.
\newblock Plenum Press, New York.
\bibitem[\protect\astroncite{Cohen}{1982}]{cohen:set}
Cohen, J. (1982).
\newblock Set correlation as a general multivariate data-analytic method.
\newblock {\em Multivariate Behavioral Research}, 17(3).
\bibitem[\protect\astroncite{Cohen et~al.}{2003}]{cohen:03}
Cohen, J., Cohen, P., West, S.~G., and Aiken, L.~S. (2003).
\newblock {\em Applied multiple regression/correlation analysis for the
behavioral sciences}.
\newblock L. Erlbaum Associates, Mahwah, N.J., 3rd ed edition.
\bibitem[\protect\astroncite{Cooksey and Soutar}{2006}]{cooksey:06}
Cooksey, R. and Soutar, G. (2006).
\newblock Coefficient beta and hierarchical item clustering - an analytical
procedure for establishing and displaying the dimensionality and homogeneity
of summated scales.
\newblock {\em Organizational Research Methods}, 9:78--98.
\bibitem[\protect\astroncite{Cronbach}{1951}]{cronbach:51}
Cronbach, L.~J. (1951).
\newblock Coefficient alpha and the internal structure of tests.
\newblock {\em Psychometrika}, 16:297--334.
\bibitem[\protect\astroncite{Dwyer}{1937}]{dwyer:37}
Dwyer, P.~S. (1937).
\newblock The determination of the factor loadings of a given test from the
known factor loadings of other tests.
\newblock {\em Psychometrika}, 2(3):173--178.
\bibitem[\protect\astroncite{Everitt}{1974}]{everitt:74}
Everitt, B. (1974).
\newblock {\em Cluster analysis}.
\newblock John Wiley \& Sons, Cluster analysis. 122 pp. Oxford, England.
\bibitem[\protect\astroncite{Fox et~al.}{2012}]{sem}
Fox, J., Nie, Z., and Byrnes, J. (2012).
\newblock {\em {sem: Structural Equation Models}}.
\bibitem[\protect\astroncite{Grice}{2001}]{grice:01}
Grice, J.~W. (2001).
\newblock Computing and evaluating factor scores.
\newblock {\em Psychological Methods}, 6(4):430--450.
\bibitem[\protect\astroncite{Guilford}{1954}]{guilford:54}
Guilford, J.~P. (1954).
\newblock {\em Psychometric Methods}.
\newblock McGraw-Hill, New York, 2nd edition.
\bibitem[\protect\astroncite{Guttman}{1945}]{guttman:45}
Guttman, L. (1945).
\newblock A basis for analyzing test-retest reliability.
\newblock {\em Psychometrika}, 10(4):255--282.
\bibitem[\protect\astroncite{Hartigan}{1975}]{hartigan:75}
Hartigan, J.~A. (1975).
\newblock {\em Clustering Algorithms}.
\newblock John Wiley \& Sons, Inc., New York, NY, USA.
\bibitem[\protect\astroncite{Henry et~al.}{2005}]{henry:05}
Henry, D.~B., Tolan, P.~H., and Gorman-Smith, D. (2005).
\newblock Cluster analysis in family psychology research.
\newblock {\em Journal of Family Psychology}, 19(1):121--132.
\bibitem[\protect\astroncite{Holm}{1979}]{holm:79}
Holm, S. (1979).
\newblock A simple sequentially rejective multiple test procedure.
\newblock {\em Scandinavian Journal of Statistics}, 6(2):pp. 65--70.
\bibitem[\protect\astroncite{Holzinger and Swineford}{1937}]{holzinger:37}
Holzinger, K. and Swineford, F. (1937).
\newblock The bi-factor method.
\newblock {\em Psychometrika}, 2(1):41--54.
\bibitem[\protect\astroncite{Horn}{1965}]{horn:65}
Horn, J. (1965).
\newblock A rationale and test for the number of factors in factor analysis.
\newblock {\em Psychometrika}, 30(2):179--185.
\bibitem[\protect\astroncite{Horn and Engstrom}{1979}]{horn:79}
Horn, J.~L. and Engstrom, R. (1979).
\newblock Cattell's scree test in relation to bartlett's chi-square test and
other observations on the number of factors problem.
\newblock {\em Multivariate Behavioral Research}, 14(3):283--300.
\bibitem[\protect\astroncite{Jennrich and Bentler}{2011}]{jennrich:11}
Jennrich, R. and Bentler, P. (2011).
\newblock Exploratory bi-factor analysis.
\newblock {\em Psychometrika}, pages 1--13.
\newblock 10.1007/s11336-011-9218-4.
\bibitem[\protect\astroncite{Jensen and Weng}{1994}]{jensen:weng}
Jensen, A.~R. and Weng, L.-J. (1994).
\newblock What is a good g?
\newblock {\em Intelligence}, 18(3):231--258.
\bibitem[\protect\astroncite{Loevinger et~al.}{1953}]{loevinger:53}
Loevinger, J., Gleser, G., and DuBois, P. (1953).
\newblock Maximizing the discriminating power of a multiple-score test.
\newblock {\em Psychometrika}, 18(4):309--317.
\bibitem[\protect\astroncite{MacCallum et~al.}{2007}]{maccallum:07}
MacCallum, R.~C., Browne, M.~W., and Cai, L. (2007).
\newblock Factor analysis models as approximations.
\newblock In Cudeck, R. and MacCallum, R.~C., editors, {\em Factor analysis at
100: Historical developments and future directions}, pages 153--175. Lawrence
Erlbaum Associates Publishers, Mahwah, NJ.
\bibitem[\protect\astroncite{Martinent and Ferrand}{2007}]{martinent:07}
Martinent, G. and Ferrand, C. (2007).
\newblock A cluster analysis of precompetitive anxiety: Relationship with
perfectionism and trait anxiety.
\newblock {\em Personality and Individual Differences}, 43(7):1676--1686.
\bibitem[\protect\astroncite{McDonald}{1999}]{mcdonald:tt}
McDonald, R.~P. (1999).
\newblock {\em Test theory: {A} unified treatment}.
\newblock L. Erlbaum Associates, Mahwah, N.J.
\bibitem[\protect\astroncite{Mun et~al.}{2008}]{mun:08}
Mun, E.~Y., von Eye, A., Bates, M.~E., and Vaschillo, E.~G. (2008).
\newblock Finding groups using model-based cluster analysis: Heterogeneous
emotional self-regulatory processes and heavy alcohol use risk.
\newblock {\em Developmental Psychology}, 44(2):481--495.
\bibitem[\protect\astroncite{Nunnally}{1967}]{nunnally:67}
Nunnally, J.~C. (1967).
\newblock {\em Psychometric theory}.
\newblock McGraw-Hill, New York,.
\bibitem[\protect\astroncite{Nunnally and
Bernstein}{1984}]{nunnally:bernstein:84}
Nunnally, J.~C. and Bernstein, I.~H. (1984).
\newblock {\em Psychometric theory}.
\newblock McGraw-Hill, New York,, 3rd edition.
\bibitem[\protect\astroncite{Pedhazur}{1997}]{pedhazur:97}
Pedhazur, E. (1997).
\newblock {\em Multiple regression in behavioral research: explanation and
prediction}.
\newblock Harcourt Brace College Publishers.
\bibitem[\protect\astroncite{Revelle}{1979}]{revelle:iclust}
Revelle, W. (1979).
\newblock Hierarchical cluster-analysis and the internal structure of tests.
\newblock {\em Multivariate Behavioral Research}, 14(1):57--74.
\bibitem[\protect\astroncite{Revelle}{2012}]{psych}
Revelle, W. (2012).
\newblock {\em psych: Procedures for Personality and Psychological Research}.
\newblock Northwestern University, Evanston.
\newblock R package version 1.2.8.
\bibitem[\protect\astroncite{Revelle}{prep}]{revelle:intro}
Revelle, W. ({in prep}).
\newblock {\em An introduction to psychometric theory with applications in
{R}}.
\newblock Springer.
\bibitem[\protect\astroncite{Revelle et~al.}{2011}]{rcw:methods}
Revelle, W., Condon, D., and Wilt, J. (2011).
\newblock Methodological advances in differential psychology.
\newblock In Chamorro-Premuzic, T., Furnham, A., and von Stumm, S., editors,
{\em Handbook of Individual Differences}, chapter~2, pages 39--73.
Wiley-Blackwell.
\bibitem[\protect\astroncite{Revelle and Rocklin}{1979}]{revelle:vss}
Revelle, W. and Rocklin, T. (1979).
\newblock {Very Simple Structure} - alternative procedure for estimating the
optimal number of interpretable factors.
\newblock {\em Multivariate Behavioral Research}, 14(4):403--414.
\bibitem[\protect\astroncite{Revelle et~al.}{2010}]{rwr:sapa}
Revelle, W., Wilt, J., and Rosenthal, A. (2010).
\newblock Personality and cognition: The personality-cognition link.
\newblock In Gruszka, A., Matthews, G., and Szymura, B., editors, {\em Handbook
of Individual Differences in Cognition: Attention, Memory and Executive
Control}, chapter~2, pages 27--49. Springer.
\bibitem[\protect\astroncite{Revelle and Zinbarg}{2009}]{rz:09}
Revelle, W. and Zinbarg, R.~E. (2009).
\newblock Coefficients alpha, beta, omega and the glb: comments on {Sijtsma}.
\newblock {\em Psychometrika}, 74(1):145--154.
\bibitem[\protect\astroncite{Schmid and Leiman}{1957}]{schmid:57}
Schmid, J.~J. and Leiman, J.~M. (1957).
\newblock The development of hierarchical factor solutions.
\newblock {\em Psychometrika}, 22(1):83--90.
\bibitem[\protect\astroncite{Shrout and Fleiss}{1979}]{shrout:79}
Shrout, P.~E. and Fleiss, J.~L. (1979).
\newblock Intraclass correlations: Uses in assessing rater reliability.
\newblock {\em Psychological Bulletin}, 86(2):420--428.
\bibitem[\protect\astroncite{Smillie et~al.}{2012}]{smillie:jpsp}
Smillie, L.~D., Cooper, A., Wilt, J., and Revelle, W. (2012).
\newblock Do extraverts get more bang for the buck? refining the
affective-reactivity hypothesis of extraversion.
\newblock {\em Journal of Personality and Social Psychology}, 103(2):306--326.
\bibitem[\protect\astroncite{Sneath and Sokal}{1973}]{sneath:73}
Sneath, P. H.~A. and Sokal, R.~R. (1973).
\newblock {\em Numerical taxonomy: the principles and practice of numerical
classification}.
\newblock A Series of books in biology. W. H. Freeman, San Francisco.
\bibitem[\protect\astroncite{Sokal and Sneath}{1963}]{sokal:63}
Sokal, R.~R. and Sneath, P. H.~A. (1963).
\newblock {\em Principles of numerical taxonomy}.
\newblock A Series of books in biology. W. H. Freeman, San Francisco.
\bibitem[\protect\astroncite{Spearman}{1904}]{spearman:rho}
Spearman, C. (1904).
\newblock The proof and measurement of association between two things.
\newblock {\em The American Journal of Psychology}, 15(1):72--101.
\bibitem[\protect\astroncite{Steiger}{1980}]{steiger:80b}
Steiger, J.~H. (1980).
\newblock Tests for comparing elements of a correlation matrix.
\newblock {\em Psychological Bulletin}, 87(2):245--251.
\bibitem[\protect\astroncite{Thorburn}{1918}]{thornburn:1918}
Thorburn, W.~M. (1918).
\newblock The myth of occam's razor.
\newblock {\em Mind}, 27:345--353.
\bibitem[\protect\astroncite{Thurstone and Thurstone}{1941}]{thurstone:41}
Thurstone, L.~L. and Thurstone, T.~G. (1941).
\newblock {\em Factorial studies of intelligence}.
\newblock The University of Chicago press, Chicago, Ill.
\bibitem[\protect\astroncite{Tryon}{1935}]{tryon:35}
Tryon, R.~C. (1935).
\newblock A theory of psychological components--an alternative to "mathematical
factors.".
\newblock {\em Psychological Review}, 42(5):425--454.
\bibitem[\protect\astroncite{Tryon}{1939}]{tryon:39}
Tryon, R.~C. (1939).
\newblock {\em Cluster analysis}.
\newblock Edwards Brothers, Ann Arbor, Michigan.
\bibitem[\protect\astroncite{Velicer}{1976}]{velicer:76}
Velicer, W. (1976).
\newblock Determining the number of components from the matrix of partial
correlations.
\newblock {\em Psychometrika}, 41(3):321--327.
\bibitem[\protect\astroncite{Zinbarg et~al.}{2005}]{zinbarg:pm:05}
Zinbarg, R.~E., Revelle, W., Yovel, I., and Li, W. (2005).
\newblock Cronbach's {$\alpha$}, {Revelle's} {$\beta$}, and {McDonald's}
{$\omega_H$}): Their relations with each other and two alternative
conceptualizations of reliability.
\newblock {\em Psychometrika}, 70(1):123--133.
\bibitem[\protect\astroncite{Zinbarg et~al.}{2006}]{zinbarg:apm:06}
Zinbarg, R.~E., Yovel, I., Revelle, W., and McDonald, R.~P. (2006).
\newblock Estimating generalizability to a latent variable common to all of a
scale's indicators: A comparison of estimators for {$\omega_h$}.
\newblock {\em Applied Psychological Measurement}, 30(2):121--144.
\end{thebibliography}
\printindex
\end{document}