%\VignetteIndexEntry{Recovering a Basic Space in R}
%\VignetteKeywords{multivariate}
%\VignettePackage{basicspace}

\documentclass[12pt]{article}
\usepackage{fullpage}

\usepackage{Sweave}
\usepackage{amsmath}
\usepackage{amscd}

\begin{document}

\title{Recovering a Basic Space in R}
%\author{James Lo}
\maketitle

\section{Introduction}

<<echo=false>>=
options(width=60)
@

This R (R Development Core Team, 2010) package contains software designed to recover a latent dimension (i.e. a basic space)
from a numeric matrix of data. The initial development was conducted by Poole (1998), who wrote the Fortran
executable that was previously used by social scientists to analyze survey data in publications such
as Palfrey and Poole (1987) and Saiegh (2009).  The principal purpose of the package is to facilitate
the analysis of self-placement and/or perceptual survey data by scaling stimuli into a common space. Stated
differently, assume that various political actors (i.e. the stimuli) occupy positions on a continuous
left-right political spectrum, and our objective is to recover estimates of those positions.  This package
facilitates the estimation of these positions under the assumption that survey respondents report their
personal estimates of these positions with different levels of bias and scale attenuation. While intended
for use with survey data, this procedure has potentially broader applications outside of social science
because mathematically it can be thought of as a special case of singular value decomposition on a numeric
matrix with missing data.\footnote{When used to analyze survey data, this property implies that
survey scales can be treated as a continuous rather than ordinal variable. Researchers who believe that
their scales are ordinal should therefore not be using the model presented here.}

This vignette proceeds in four steps. First, we begin with a description of the math that underlies the
basic space estimator. We then provide three examples. First, we show a Monte Carlo analysis that suggests
our estimator produces an accurate decomposition of our simulated data matrix even with 30 per cent of the
data missing. Secondly, we show how the procedure can be applied to self-placement survey data from the 
1980 National Election Survey. We conclude with an application of the model to perceptual data from
the 1980 National Election Study, where various political candidates are ranked along a 7 point
liberal-conservative scale.


\section{Model}

The exposition of the model presented here closely follows Poole (1998). Consider a matrix of survey
data $X_0$ with $n$ respondents and $m$ issue scales, with individuals on the rows and issues on the
columns. Some cells of the matrix $X_0$ are missing, and we let $X$ denote the version of $X_0$ that
has no missing data. In each cell $x_{ij}$, respondent $i$ (\textit{i=1, \ldots, n}) reports their position on issue scale $j$ (\textit{j=1, \ldots, m}), with some responses missing.\footnote{The `0' subscript
indicates that some elements are missing from the matrix.} Now let $\Psi_{ik}$ be the $i$th
individual's position on the $k$th basic dimension (\textit{k=1, \ldots, s}), $W$ be an $m$ by $s$
matrix of weights that map individual positions from the basic space to the issue dimensions, $c$
be a vector of issue dimension intercept terms of length $m$, $J_n$ by an $n$ length vector of ones,
and $E_0$ be error terms in the data matrix. The model that we seek to estimate is:

\begin{equation}
X_0 = [\Psi W' + J_nc']_0 + E_0
\label{one}
\end{equation}

Without loss of generality, we also assume that $E_0$ is drawn from a symmetric distribution with
mean 0 and that the centroid of the basic space coordinates is at the origin (i.e. $J_n'\Psi=0$).
Substituting into the model equation, this implies that $J_n'[X-J_nc'] = 0_m'$, where $0_m$ is an $m$
length vector of zeroes. Then in the situation where $X_0$ has no missing data, the parameters
of interest can all be recovered using singular value decomposition.  To see why this is true, recall
that for an $n$ by $m$ matrix of real elements with $n \ge m$, there exists an $n$ by $m$ orthogonal
matrix $U$, an $m$ by $m$ orthogonal matrix $V$, and an $m$ by $m$ matrix $\Lambda$ such that:

\begin{equation}
X = U \Lambda V'
\label{two}
\end{equation}

\noindent where $\Lambda$ is a diagonal matrix of singular values.\footnote{A more general
form of this equation can be written in which $\Lambda$ is instead an $n$ by $m$ matrix and $U$ is
$n$ by $n$.} To solve [1], set $c$ equal to the column means of $X$, or $c_j = \frac{\sum\limits_{i=1}^n x_{ij}}{n} = \bar{x}_j$. Then using [2], the 
singular value decomposition of $X - J_nc'$ can be expressed as:

$$X - J_nc'\ =\ U \Lambda V'\ =\ \Psi W'$$

This implies that in the absence of missing data, one solution for $\Psi$ and $W$ is:

$$\Psi = U\Lambda^{0.5}$$

$$W = V\Lambda^{0.5}$$

\noindent with $\Lambda^{0.5}$ being a diagonal matrix where diagonal elements are
the square roots of $\Lambda$. While other solutions to this problem exist, Eckart
and Young (1936) have shown that the least squares approximation in $s$ dimensions
of a matrix $A$ can be found by using only the first $s$ singular values of $A$ along the
diagonal of $\Lambda$ and remultiplying $U \Lambda V'$.

In the presence of missing data in data matrix $X_0$, the use of singular value decomposition
to solve for $W$ and $\Psi$ is no longer possible, and we instead estimate $\hat{W}$ and
$\hat{\Psi}$ using an alternating least squares (ALS) technique that is similar to the procedures
used in Carroll and Chang (1970) and Takane, Young, and de Leeuw (1977). The objective
function to be minimized is the sum of the squared deviations across all cells in $A$ after
the columns have been adjusted for column means, or:

$$\xi = \displaystyle\sum\limits_{i=1}^n \sum\limits_{j=1}^{m_i} \{[ \sum\limits_{k=1}^s \Psi_{ik}W_{jk}] + c_j -x_{ij} \}^2$$

In minimizing this objective function, two constraints from the earlier analysis with no missing
data are applied. First, we exploit the the fact that $\Psi$ and $W$ are orthogonal matrices, which implies
that $\Psi'\Psi = W'W$.\footnote{More specifically, $\Psi'\Psi = \Lambda^{0.5}U'U\Lambda^{0.5} = \Lambda^{0.5}I_m\Lambda^{0.5} = \Lambda = W'W.$} Secondly, following our earlier restriction that $J_n'[X-J_nc'] = 0_m'$, $J_n'U = J_n'\Psi = 0_m'$ as well.
These restrictions produce the Lagrangian multiplier problem:


$$\mu = \xi = 2\gamma'['\Psi'J_n] + tr[\Phi(\Psi'\Psi - W'W)]$$

\noindent where $\Phi$ is a symmetric $s$ by $s$ matrix of Lagrangian multipliers and $\gamma$ is an $s$ length
vector of Lagrangian multipliers. Since all Lagrangian multipliers are zero,\footnote{See Appendix A in Poole (1998)
for a full proof that all Lagrangian multipliers are zero.} the partial derivatives of $\xi$ are:

\begin{equation}
\frac{\partial \mu}{\partial \Psi_{ik}} = 2\displaystyle\sum\limits_{j=1}^{m_i} [(\sum\limits_{l=1}^s w_{jl}\psi_{jl})+ c_j - x_{ij}]w_{jk}
\label{three}
\end{equation}

\begin{equation}
\frac{\partial \mu}{\partial w_{jk}} = 2\displaystyle\sum\limits_{i=1}^{n_j} [(\sum\limits_{l=1}^s w_{jl}\psi_{jl})+ c_j - x_{ij}]\psi_{jk}
\label{four}
\end{equation}

\begin{equation}
\frac{\partial \mu}{\partial c_j} = 2\displaystyle\sum\limits_{i=1}^{n_j} [(\sum\limits_{l=1}^s w_{jl}\psi_{jl})+ c_j - x_{ij}]
\label{five}
\end{equation}

Let $W^*$ be an $m_i$ by $s$ matrix with appropriate rows corresponding to missing entries in $X_0$ removed, $x_{0i}$ be the
length $m_i$ row of $X_0$, and $c_0$ be the length $m_i$ vector of constants corresponding to the elements of $x_{0i}$. Then
if $W^*$'$W^*$ exists, the $ith$ row of $\Psi$ can be estimated by setting (\ref{three}) to zero, collecting the $s$ partial derivatives
of the $i$th row of $\Psi$ into a vector and solving for $\psi_j$ as:

\begin{equation}
\hat{\psi}_i = (W^*\mathrm{'}W^*)^{-1} W^*\mathrm{'}[x_{0i} - c_{0}]
\label{six}
\end{equation}

\noindent which can of course by estimated using ordinary least squares. Similarly, let $\Psi_j^* = [\Psi_0|J_0]$ be an
$n_j$ by $s+1$ matrix with the appropriate rows corresponding to missing data removed and bordered by ones, $w_j$ be
the $s$ length vector of row $j$ in $W$, $c_j$ be the $j$th element of $c$, and $x_{0j}$ be the $j$th column of $X_0$.
Then if $\Psi_{j}^*$$'\Psi_j^*$ exists, $w_j$ and $c_j$ can be jointly estimated by combining (4) and (5) as:

\begin{equation}
\frac{\hat{w}_j}{\hat{c}_j} = (\Psi_j^*\mathrm{'}\Psi_j^*)^{-1} \Psi_j^*\mathrm{'}x_{0j}
\label{seven}
\end{equation}

Equations (\ref{six}) and (\ref{seven}) represent the core set of equations that are used to solve for $\hat{W}$,
$\hat{c}$, and $\hat{\Psi}$. Once a set of starting values has been generated, (\ref{six}) and (\ref{seven}) are
iterated until convergence is complete. Generation of appropriate start values is conducted one dimension at a time,
and a more detailed justification of the procedure can be found in Poole (1998). On the first dimension, start
values are generated by using the following three equations:

\begin{equation}
\hat{c}_{j} = \frac{\displaystyle\sum\limits_{j=1}^{n_j} x_{ij}}{n_j} = \bar{x}_j
\label{eight}
\end{equation}

\begin{equation}
w_{j1} = diag(\Gamma)
\label{nine}
\end{equation}

\noindent where $\Gamma$ is an $m$ by $m$ diagonal matrix with diagonal elements either set to 1 or -1
that maximizes the number of positive elements in the $m$ by $m$ covariance matrix
$\Gamma[X_0 - J_nc']\mathrm{'}[X_0 - J_nc']\Gamma$. $\Gamma$ is found by a simple iterative process
similar to that used to speed eigenvector/eigenvalue decomposition (Poole, 1998). Given equations (\ref{eight})
and (\ref{nine}), starting values for $\psi$ are:

\begin{equation}
\hat{\psi}_{i1} = \frac{\displaystyle\sum\limits_{j=1}^{m_i} \hat{w}_{j1}(x_{ij}-\hat{c}_j)}{m_i}
\label{ten}
\end{equation}

If more than one dimensions are to be estimated ($s > 1$), start values for other dimensions can be generated
simply by replacing the data matrix $X_0$ with the matrix of residuals $E_{0s}$ in equations (\ref{nine}) and (\ref{ten}). However, no further estimation of start values for $\hat{c}$ is required. The matrix of residuals to be used for generating start values on dimension $s$ is:

$$E_{0s} = X_0 - \displaystyle\sum\limits_{j=1}^{s} \hat{\Psi}_s\hat{w}_s\mathrm{'} - J_n\hat{c}'$$

This residual matrix allows the generation of higher-dimension start values by iterating $\Gamma$ to
maximize the positive elements in $E_{0s}$. The starting values are now:

%Equation 12 in paper
\begin{equation}
\hat{\psi}_{is} = \frac{\displaystyle\sum\limits_{j=1}^{m_i} \hat{w}_{js}e_{1(s-1)ij}}{\sum\limits_{j=1}^{m_j} \hat{w}_{js}^2}
\label{eleven}
\end{equation}

\noindent where the initial $\hat{w_{js}}$ values of +1s and -1s are used to obtain $\hat{\psi}_{is}$ starting values. The starting values of $\hat{w_{js}}$ are now:

%Equation 11 in paper
\begin{equation}
\hat{w}_{js} = \frac{\displaystyle\sum\limits_{i=1}^{n_j} \hat{\psi}_{is}e_{(s-1)ij}}{\sum\limits_{i=1}^{n_j} \hat{\psi}_{is}^2}
\label{twelve}
\end{equation}


Summarizing the preceding discussion in full, the basic space technique decomposes an $n$ by $m$ matrix $X_0$ with
$n > m$ following equation (\ref{one}). In the absence of missing data, this decomposition can be solved using
singular value decomposition as in equation (\ref{two}). The basic space technique can therefore be thought of as
a generalization of singular value decomposition to matrices with missing data. Estimation of (\ref{one}) proceeds
in three steps. In the first stage, starting values on the first dimension are generated for $\hat{c}_{j}$, $w_{j1}$,
and $\hat{\psi}_{i1}$ by iterating equations (\ref{eight})-(\ref{ten}) until convergence. In the second stage, if the
number of dimensions to be estimated $s > 1$, higher dimensional starting values for $\hat{\psi}_{is}$ and $\hat{w}_{is}$
are generated dimension by dimension using equations (\ref{eleven})-(\ref{twelve}). Finally, the starting values generated 
in the preceding two stages are improved by iterating equations (\ref{six})-(\ref{seven}) until convergence.

\section{Monte Carlo Test}

In this section, we present the first of three motivating examples. We begin with a Monte Carlo example that
tests the basic space technique against simulated data. Four key variables should be set in each simulation:
the number of respondents $N$ (set here to $N=1000$), the number of issue scales (also referred to as stimuli,
and set here to $M=20$), the number of explantory dimensions (set here to $s=2$), the
fraction of observations that are missing (set here as 0.3), and the distribution of error terms (set here as
random uniform draws from -0.5 to 0.5). These variables can be changed for other simulations, but the
restriction that $N>M$ must be hold true. In cases where $M>N$ please refer to the third example that uses
\emph{blackbox\_transpose()} and \emph{aldmck()}.

<<one>>=
set.seed(1231)
library(basicspace)
N <- 1000
M <- 20
s <- 2
fraction.missing <- 0.3
E <- matrix(runif(N*M, min=-0.5, max=0.5),nrow=N,ncol=M)
@

To generate the $X$ matrix (i.e. the matrix in (\ref{one}) before missing values are introduced), separately
generate the matrices that produce the singular value decomposition of $X$ following (\ref{two}). Also
generate the $J_n$ and $c$ vectors from (1). While $X$ can be generated directly in one step, creating the
components separately enjoys two significant advantages. First, recovery of the true values of $\Psi$ and $W$
is simplified. Secondly, the creation of $\Lambda$ separately allows us to more easily tune the dimensionality
of the matrix as desired.

<<two>>=
U <- matrix(runif(N*s),nrow=N,ncol=s)
D <- diag(seq(from=2.1, by=-0.2,length.out=s))
V.prime <- matrix(runif(s*M),nrow=s,ncol=M)
c <- rnorm(M)
Jn <- rep(1,N)
@

With the intermediate matrices just generated, we can produce our $X$ matrix by using equation (\ref{one})
and the true $\Psi$ and $W$ matrices using: $\Psi = U\Lambda^{0.5}$ and $W = V\Lambda^{0.5}$. 

<<three>>=
X.true <- U %*% D %*% V.prime + Jn %o% c
X.0 <- X.true + E
Psi.true <- U %*% sqrt(D)
W.true <- t(V.prime) %*% sqrt(D)
@

$X_0$ is simply the $X$ matrix with missing data values included completely at random, so we insert
our missing data code (999 in the example) into the appropriate fraction of values as follows:

<<four>>=
missing <- sample(1:(N*M), round(fraction.missing*N*M))
X.0[missing] <- 999
@

The final step before estimation is to assign row and column names to the data set prior to input. In
most applications these names are generally pulled from a survey, but they can also be generated
manually:

<<five>>=
rownames(X.0) <- paste("Legis", 1:N, sep="")
colnames(X.0) <- paste("V", 1:M, sep="")
@

Estimation of the Monte Carlo data after formatting is trivial. The function that applies the basic space
decomposition described in this paper is \emph{blackbox()}.  It takes four arguments: the matrix to
be decomposed, a vector of missing data values, a boolean flag indicating whether verbose output is
desired, the number of dimensions to estimate, and the minimum number of issue scales that an individual
needs to provide responses to if they are to be included in the estimation.

<<six>>=
result <- blackbox(X.0,missing=c(999),verbose=TRUE,dims=2,minscale=8)
names(result)
@

The output object contains multiple data frames summarizing the results of the estimation. The key data
frames are \textbf{stimuli}, which contain estimates of $\hat{W}$ and $\hat{c}$, as well as $\textbf{individuals}$,
which contain estimates of $\hat{\Psi}$.  The other quantities are fit statistics described in greater
detail in the standard documentation for the function.

With the estimates complete, we are now able to test the recovery of our parameters of interest. In general,
scaling problems are not fully identified. Stated differently, given $X = \Psi W'$, $\Psi$ and $W'$ are not
unique solutions because $X = \Psi K K^{-1} W'$ for any conformable and invertible matrix $K$, so $X$ can
always be decomposed instead as $X = \Psi^* W^{*'}$ where $\Psi^* = \Psi K$ and $W^{*'} = K^{-1} W'$. 
When evaluating parameter fit, we are therefore largely concerned with finding monotonic relationships between
the true and estimated parameters of interest.  Figure \ref{fig:seven} compares the true vs. estimated values of
$\Psi$ across two dimensions, and the results suggest a reasonable model fit.

\begin{figure}
\begin{center}
<<label=seven,fig=TRUE,echo=TRUE,width=14,height=7>>=
par(mfrow=c(1,2))

plot(Psi.true[,1], result$individuals[[2]]$c1, xlim=c(0, 1.5), ylim=c(-0.65, 0.6),
pch=20,cex=0.4, cex.lab=1.6, bty="n",
xlab="True Psi, first dimension", ylab="Recovered Psi, first dimension")

plot(Psi.true[,2], result$individuals[[2]]$c2, xlim=c(0, 1.5), ylim=c(-0.65, 0.6),
pch=20,cex=0.4, cex.lab=1.6, bty="n",
xlab="True Psi, second dimension", ylab="Recovered Psi, second dimension")
@
\end{center}
\caption{Plots of True vs. Estimated $\Psi$ scores, first and second dimension.}
\label{fig:seven}
\end{figure}

Figure \ref{fig:six} shows the results for the same procedure applied to $W$. In Figure
\ref{fig:eight} we repeat this analysis for $c$, which is a column mean that is only estimated in
one dimension.  In both cases the estimates for $\hat{W}$ and $\hat{c}$ are a monotonic transformation of the
true parameters as expected.\footnote{In other estimates, the relationship may only be affine because $X = \Psi W'$
implies $X = -(\Psi) -(W')$ as well.}

\begin{figure}
\begin{center}
<<label=six,fig=TRUE,echo=TRUE,width=14,height=7>>=
par(mfrow=c(1,2))

plot(W.true[,1], result$stimuli[[2]]$w1, ylim=c(1, 2.6),
pch=20, cex=1.5, cex.lab=1.6, bty="n",
xlab="True W, first dimension", ylab="Recovered W, first dimension")

plot(W.true[,2], result$stimuli[[2]]$w2, ylim=c(-2, 2),
pch=20,cex=1.5, cex.lab=1.6, bty="n",
xlab="True W, second dimension", ylab="Recovered W, second dimension")

@
\end{center}
\caption{Plots of True vs. Estimated $W$ scores, first and second dimension.}
\label{fig:six}
\end{figure}


\begin{figure}
\begin{center}
<<label=eight,fig=TRUE,echo=TRUE,width=4,height=4>>=

par(mfrow=c(1,1))

plot(c, result$stimuli[[2]]$c,
pch=20,cex=1.2, cex.lab=1.1, bty="n",
xlab="True C", ylab="Recovered C")
@
\end{center}
\caption{Plot of True vs. Estimated $c$ scores.}
\label{fig:eight}
\end{figure}

Finally, we pool our estimates of $\hat{W}$, $\hat{Psi}$, and $\hat{c}$ together to
estimate the full matrix $\hat{X}$ following Equation ($\ref{one}$).  While social
scientists are principally concerned with estimation of $\hat{W}$ and $\hat{Psi}$,
others seeking to conduct singular value decomposition of matrices with missing
data may find $\hat{X}$ to be of value. One obvious application of $\hat{X}$ is
its potential use as an imputation tool for missing data.\footnote{The simulation
presented here simulates missing data under the Missing Completely at Random (MCAR)
assumption --- nevertheless, there is no reason to think that this would not work
under conditions where data are instead Missing at Random (MAR).} To test the viability
of this idea, we separately plot the true values of $X$ against the estimated values
of $\hat{X}$ separately for the cells retained in the estimation, and compared
those results to estimates of $\hat{X}$ in cells that were discarded prior to estimation
to simulate the missing data mechanism. Figure \ref{fig:nine} presents our results for retained vs.
imputed X. What is particularly notable about this result is the close similarity
between these plots --- the imputed values not only appear reasonable (i.e. line up
with the true values along a $45^{\circ}$ line), but imputed values do not appear to have
significantly higher mean squared error than the values that were retained (i.e. variance
along the $45^{\circ}$ line is similar in both plots). While further tests are necessary,
these results suggest that the use of the techniques demonstrated here may have greater
applicability beyond survey research.  Further discussion of imputation can be found in the
unpublished appendix to Poole (1998).


\begin{figure}
\begin{center}
<<label=nine,fig=TRUE,echo=TRUE,width=14,height=7>>=

W.hat <- cbind(result$stimuli[[2]]$w1, result$stimuli[[2]]$w2)
Psi.hat <- cbind(result$individuals[[2]]$c1, result$individuals[[2]]$c2)
X.hat <- Psi.hat %*% t(W.hat) + Jn %o% result$stimuli[[2]]$c

par(mfrow=c(1,2))

plot(X.true[missing], X.hat[missing],
pch=20,cex=0.4, cex.lab=1.2, bty="n",
xlab="True X, missing values", ylab="Recovered X, missing values")

plot(X.true[!(1:(N*M) %in% missing)], X.hat[!(1:(N*M) %in% missing)],
pch=20,cex=0.4, cex.lab=1.2, bty="n",
xlab="True X, nonmissing values", ylab="Recovered X, nonmissing values")
@
\end{center}
\caption{Plots of True vs. Estimated $X$ scores for missing vs. nonmissing values.}
\label{fig:nine}
\end{figure}


\section{Recovering a Basic Space with the 1980 National Election Study Issue Scales}

In this section we present an application of the basic space model to a common problem in social science.
One issue of interest to political scientists in particular are empirical models of voting (also
known as ideal point models), that allow legislator locations in an abstract policy or ideological space
to be inferred from their roll call votes. The recovered scores have wide applicability to the study
of Congress (Poole and Rosenthal 1997; McCarty et al. 2006), elections (Herron and Lewis 2007),
courts (Martin and Quinn, 2002), and in non-legislative voting bodies such as the United
Nations (Voeten 2000).\footnote{For a more extensive review of applications of spatial modelling
in the social sciences, see Poole (2005).}  The most prominent ideal point model in the political
science literature is NOMINATE (Poole and Rosenthal 1997), a model that estimates the policy
preferences of legislators using observed roll call votes as the primary source of data.\footnote{The
package \emph{wnominate} on CRAN contains software used to estimate NOMINATE scores.}  However, the
use of models such as NOMINATE may not always be possible because roll call data is often not
available or recorded.  In such instances, the basic space model shown here presents an attractive
alternative estimator.\footnote{For a more comprehensive review of the advantages and disadvantages
of different data sources for spatial models, see Saiegh (2009).}

We present a simple example that applies the basic space model to a set of issue scales from the
1980 National Election Study. This survey contains N=1614 respondents who were asked to place themselves
on scales about desired levels of defense spending, inflation, tax cuts, abortion, liberal-conservative
scales, the role of women, the role of government in providing jobs, busing, and other similar issues.
We assume that each respondent has a location in a common ideological space and attempt to recover
estimates of those locations, which is represented as $\Psi$ in (\ref{one}). In providing responses
to the issue scales, each respondent reports their true ideological position $\Psi$, modified by a
stretch parameter $W$ and an additive intercept parameter $c$ along with noise $E_0$. The data is simply
stored in a standard matrix or data frame with respondents on the rows and survey questions (i.e. stimuli)
on the columns as follows:

<<ten>>=
data(Issues1980)
Issues1980[1:10,1:4]
@

Virtually all surveys contain missing data, and for the two survey questions about abortion, `7' is
used as a missing data code. However, many of the other scales in this data set us 7 point scales, so
we need to recode the missing data for those questions.

<<eleven>>=
Issues1980[Issues1980[,"abortion1"]==7,"abortion1"] <- 8
Issues1980[Issues1980[,"abortion2"]==7,"abortion2"] <- 8
@

Estimation of the scores is now trivial using the \emph{blackbox()} function, which takes the same
arguments already described in the Monte Carlo example:

<<twelve>>=
result <- blackbox(Issues1980,missing=c(0,8,9),verbose=FALSE,dims=3,minscale=8)
@

Objects of class \emph{blackbox()} can also be summarized using the \emph{summary()} function, although
the summaries largely provide only summaries of the stimuli. For each dimension estimated, the summary
provides the intercept ($c$) and stretch ($w_1 \ldots w_3$) parameters for each question, as well
as the number of respondents and various fit statistics.

<<thirteen>>=
summary(result)
@

When using \emph{blackbox()} for applied research, the researcher's principal goal is the recovery of
the individual parameters stored as the \emph{individuals} data frame. These typically represent our
estimate of the individual's ideological measure in the basic space. Due to the issues with model identification
discussed earlier, these measures are measured only up to an affine transformation of the true space.
In particular, the rotation of the estimate is not specified, so if the ideological measure is to be
substantively measured as a liberalism/conservatism score, its rotation should be validated so that it
can be transformed if necessary.  Here we conduct such a check by correlating our recovered scores with
self-reported liberal-conservative scores, where higher scores indicate higher levels of conservatism.
The correlation is negative, suggesting that as the recovered scores increase, the respondents become more liberal.
Since the norm in political science research is to orient liberal-conservative scores to increase as
conservatism increases, the researcher may wish to rotate the scores (i.e. by multiplying them by -1) before
using them for auxiliary analyses. 

<<fourteen>>=
cor(result$individuals[[1]]$c1, Issues1980[,"libcon1"])
@

\section{Recovering a Basic Space with the 1980 National Election Study Liberal-Conservative Scales}

In our previous example applying the basic space model to analyze respondent self-placement
on issues scales, we considered an example where the bias and stretch parameters $c$ and $w$ were
estimated for the column parameters.  However, we may instead wish to estimate a version of the model
where $c$ and $w$ are estimated for the row parameters (i.e. the survey respondents) instead.  This is
simply a transposed version of the basic space model, where $m > n$ instead of $n > m$.  A transposed
model may be reasonable in cases where the survey data to be used is perceptual data of the stimuli
of interest.  In this example we analyze perceptual data from the 1980 National Election Study. A total
of N=888 respondents were asked to place six stimuli (Carter, Reagan, Kennedy, Anderson, the Republicans,
and the Democrats) on a 7 point liberal-conservative scale. Our objective is to estimate the locations
of the six stimuli in the basic space, which each respondent perceives with some bias and stretch parameter.
The data is input in a manner identical to before, with survey respondents on the rows and stimuli on
the columns. One very important difference between \emph{blackbox()} and \emph{blackbox\_transpose()} is that
in most survey data sets, the number of respondents is very large relative to the number of stimuli.
This typically means that \emph{blackbox\_transpose()} takes much longer to estimate because each respondent
estimates both a bias $c$ and stretch $W$ parameter. To estimate the 1980 liberal-conservative placements
using \emph{blackbox\_transpose()}, we simply load the data and call the estimator as follows:

<<blackbt>>=
data(LC1980)
LCdat=LC1980[,-1]
LCdat[1:10,]
result <- blackbox_transpose(LCdat,missing=c(0,8,9),dims=3,minscale=5,verbose=TRUE)
@

In an effort to simplify interpretation of results from \emph{blackbox\_tranpose()}, we include two plot functions.
These functions plot the location of the stimuli against a probability and cumulative distribution plot of
locations of the population weights. These weights represent reasonable estimates of the two-dimensional
ideal points, as shown in Abrajano and Poole (2010).

\begin{figure}
\begin{center}
<<label=bbt,fig=TRUE,echo=TRUE,width=14,height=7>>=

par(mfrow=c(1,2))

plot(result)
plotcdf.blackbt(result)
@
\end{center}
\caption{Blackbox Transpose PDF and CDF plots.}
\label{fig:bbt}
\end{figure}

We can also produce summary reports of the stimuli as follows:

<<blackfinal>>=
summary(result)
@


\section{Scaling Perceptual Data using Aldrich and McKelvey's Esimator}

The transposed basic space model is a generalization of a model developed by Aldrich and McKelvey (1977). The transposed
basic space model can be applied to matrices with missing values and can estimate more than one dimension. For historical purposes, we include
the original Aldrich-McKelvey estimator with this package. In the basic model, the estimator assumes that:

$$Y_{ij} = Z_j + \epsilon_{ij}$$

\noindent where $Z_j$ is the true location of $j$ amd $\epsilon_{ij}$ is a random variable with mean 0, positive variance that is independent of $i$ and $j$ (homoskedastic), and zero covariance across the $i$'s and $j$'s. Aldrich and McKelvey then introduce two distortion parameters, $c_i$ and $w_i$, that transform the perceived candidate position into a reported candidate position $Z_{ij}$, according to:

$$Z_{ij}= \frac{1}{w_i}(Y_{ij} - c_i)$$

A least-squares minimization procedure is then used to obtain estimates of $\{Z_j\}_{j=1}^{J}$ and $\{w_i,c_i\}_{i=1}^{I}$.

We begin by reestimating the earlier results using the 1980 Liberal-Conservative scales using the Aldrich-McKelvey estimator. While the \emph{aldmck()} function accepts nearly identical arguments the \emph{blackbox\_transpose()}, one notable difference appears by default.  \emph{aldmck()} also accepts a column in the data matrix, specified by the \emph{respondent} argument, that specifies the respondent's self placement on the issue scale. The reported respondent rating is then transformed into an ideology score by applying the respondent's personal stretch and bias parameters to that score. Note that the results largely correspond to those shown earlier with \emph{blackbox\_transpose()}.

<<aldmck>>=
data(LC1980)
result <- aldmck(data=LC1980, polarity=2, respondent=1, missing=c(0,8,9),verbose=TRUE)
summary(result)
@



\begin{figure}
\begin{center}
<<label=aldmck_plot,fig=TRUE,echo=TRUE,width=14,height=14>>=

plot.aldmck(result)
@
\end{center}
\caption{Aldrich-McKelvey plots.}
\label{fig:aldmck_plot}
\end{figure}



Estimation of uncertainty for estimates using Aldrich-McKelvey can be obtained via the non-parametric bootstrap (Efron and Tibshirani 1994). To simulate 20 samples from the 1980 Liberal-Conservative scales, we do the following:

<<aldmck_boot>>=
Ntrials <- 20
results <- vector("list", Ntrials)
for(i in 1:Ntrials) results[[i]] <- aldmck(data=LC1980[sample(1:nrow(LC1980),nrow(LC1980),replace=TRUE),], polarity=2, respondent=1, missing=c(0,8,9),verbose=FALSE)
@

The Aldrich-McKelvey function is primarily intended for scaling perceptual data from surveys, though it can also be used to replicate previously published Monte Carlo results from Palfrey and Poole (1987). Palfrey and Poole find that the Aldrich-McKelvey algorithm is robust to the assumption of homoskedastic error, and test this by replacing the homoskedastic error term $\epsilon_{ij}$ with a respondent-specific $\epsilon_{i}$. In this example we replicate their result in a single trial, and show that the estimated rank ordering of the stimuli $\hat{Z_j}$ is identical to the true $Z_j$.

<<aldmck_MC>>=
Nstimuli <- 6
Nresp <- 500

Z_j <- rnorm(6)
Z_j <- (Z_j - mean(Z_j))/sd(Z_j)

respondent.sd <- runif(Nresp,min=0.3,max=0.9)
error_heteroskedastic <- matrix(NA,Nresp,Nstimuli)
for(i in 1:Nresp) error_heteroskedastic <- rnorm(Nstimuli,sd=respondent.sd)

w_i <- runif(Nresp,min=0,max=1)
c_i <- rnorm(Nresp)
Y_ij <- rep(1,500) %o%  Z_j
Y_ij <- Y_ij + error_heteroskedastic
Z_ij <-  1/w_i %o% rep(1,Nstimuli) * (Y_ij - c_i %o% rep(1,Nstimuli))

result <- aldmck(Z_ij,polarity=6,missing=c(999))

rank(Z_j)
rank(result$stimuli)
@

Although we only show one trial in this paper, the result shown here is robust to tests over multiple simulations.\footnote{Replication of the Monte Carlo test with separate groups of informed and uniformed individuals, from Palfrey and Poole (1987, pg. 515), can be conducted by simply changing the respondent's error deviations in the code above to have 250 respondents with $\sigma_i=0.3$ and 250 respondents with $\sigma_i=0.9$.}





\section{Conclusion}

Social scientists often wish to infer the locations of voters and legislators in an abstract policy
or ideological space.  In situations where this cannot be accomplished with a preference data-oriented
estimator such as NOMINATE (Poole and Rosenthal, 1997), the use of perceptual data-oriented
estimators such as the basic space technique described here are a useful alternative. In fact, it is
somewhat surprising that given the abundance of perceptual data questions found in most political science
surveys, articles that use perceptual data are still uncommon in the social science literature. By providing an
R (R Development Core Team, 2010) package that facilitates the analysis of perceptual data in a popular
statistics environment, we hope to encourage a renewed interest in this important literature.

\newpage
\begin{thebibliography}{3}

\bibitem{aldmck} Aldrich, John H. and Richard D. McKelvey (1977). ``A Method of Scaling with Applications to the
1968 and 1972 Presidential Elections.'' \emph{American Political Science Review} 71(1), 111-130.

\bibitem{Carroll} Carroll, Douglas and Jih-Jie Chang (1970). ``Analysis of Individual Differences in Multidimensional
Scaling Via an N-Way Generalization of ``Eckart-Young'' Decomposition.'' \emph{Psychometrika} 35:283-319.

\bibitem{EckYoung} Eckart, Carl and Gale Young (1936). ``The Approximation of One Matrix By Another of Lower Rank.''
\emph{Psycometrika} 1, 211-218.

\bibitem{EfronTibshirani} Efron, Bradley and Robert Tibshirani (1993). \emph{An Introduction to the Bootstrap.} New York: Chapman and Hall.

\bibitem{HerronLewis} Herron, M. C. and J. B. Lewis (2007). ``Did Ralph Nader Spoil Al Gore's Presidential Bid?
A Ballot-Level Study of Green and Reform Party Voters in the 2000 Presidential Election.''
\emph{Quarterly Journal of Political Science} 2, 205-226.

\bibitem{MartinQuinn} Martin, A. D. and K. Quinn (2002). ``Dynamic Ideal Point Estimation via Markov Chain Monte
Carlo for the U.S. Supreme Court, 1953-1999.'' \emph{Political Analysis} 10, 134-153.

\bibitem{MPR} McCarty, Nolan., Keith T. Poole, and Howard Rosenthal (2006). \emph{Polarized America: The Dance of
Ideology and Unequal Riches}. Cambridge: MIT Press.

\bibitem{PalfreyPoole} Palfrey, Thomas R. and Keith T. Poole (1987). ``The Relationship between Information,
Ideology, and Voting Behavior.'' \emph{American Journal of Political Science} 31(3), 511-530. 

\bibitem{Poole98} Poole, Keith T. (1998). ``Recovering a Basic Space From a Set of Issue Scales.''
\emph{American Journal of Political Science} 42(3), 954-993. 

\bibitem{Poole05} Poole, Keith T. (2005). \emph{Spatial Models of Parliamentary Voting.}
Cambridge: Cambridge University Press.

\bibitem{PR1997} Poole, Keith T. and Howard Rosenthal (1997).
\emph{Congress: A Political-Economic History of Roll Call Voting.}
New York: Oxford University Press.

\bibitem{R} R Development Core Team (2010). \emph{R: A Language and Environment for Statistical Computing.}
Vienna: R Foundation for Statistical Computing \emph{http://www.R-project.org}.

\bibitem{Saiegh} Saiegh, Sebastian (2009). ``Recovering a Basic Space from Elite Surveys: Evidence from Latin
America.'' \emph{Legislative Studies Quarterly}  34 (1), 117-145.

\bibitem{Takane} Takane, Yoshio, Forrest Young, and Jan de Leeuw (1977). ``Nonmetric Individual Differences
Multidimensional Scaling: An Alternative Least Squares Method with Optimal Scaling Features.'' \emph{Psychometrika}
42, 7-67.

\bibitem{Voeten} Voeten, Erik (2000). ``Clashes in the assembly.'' \emph{International Organization} 54 (2), 185-215.


\end{thebibliography}

\end{document}

