reduction and is also primarily used for high dimensional problems and
visualization. In brief, UMAP is also interested in proceeding locally to
learn the structure in the original high dimensional space and then project it
in a lower dimensional subspace. Yet, a key difference between UMAP and t-SNE
is the ability to reproduce a projection solution. Whereas t-SNE is interested
in assuming similarities equate to probabilities based on the overlaying of
t-distributions in the lower dimensional projection process, UMAP avoids the
notion of probabilistic similarities altogether and instead relies on smooth
geometric distances. More firmly rooted in mathematical theory and
specifically manifold learning, UMAP seeks to first learn the shape and
structure of the high dimensional manifold, and then project it locally.
Though the goal was similar to the goal in t-SNE, UMAP updates the manifold
learning process by searching two regions in the higher dimensional space:
between points, and then around regions of points.

UMAP essentially conducts two searches in order to first learn the high
dimensional structure and then be able to project it locally. First, the
algorithm searches for local connection in a pointwise fashion across the full
data space. As observations exist within this search region, they are
connected to each other. Then, the second search is in the ambient space,
previously referenced in the LLE section. This search involves taking the
nerve of the now-learned higher dimensional manifold. Practically, this means
the algorithm places a fuzzy secondary search region around each point, such
that densely populated regions have a smaller sized neighborhood (defined by
the radius of the region), compared to less densely populated regions along
the manifold, which have a larger radius around them to ensure the manifold is
connected, at least locally. This assumption of at least local connection is a
core assumption in UMAP, and one that is potentially problematic, as perhaps
not all observations truly belong to the same manifold. Still, the current
implementation of UMAP defends this mathematically in McInnes, Healy and
Melville (2018), so we continue assuming it is true. The goal with these
search regions, though, is to balance between local and global structure.
Formally, the clearest look at UMAP’s update of t-SNE is in the cost function.
UMAP’s cost function updates Equation 4.2,

$\sum_{x\in\mathcal{X}}x_{h}\text{log}(\frac{x_{h}}{x_{l}})(1-x_{h})\text{log}(\frac{1-x_{h}}{1-x_{l}}).$
(4.3)

Now, the first term, $x_{h}\text{log}(\frac{x_{h}}{x_{l}})$ which was also in
the calculation of cost in t-SNE, measures the local structure, the addition
of the second term, $(1-x_{h})\text{log}(\frac{1-x_{h}}{1-x_{l}})$, in
concerned with the learning the full shape of the manifold.

In addition to updating the cost function, we also update the definitions,
where we are no longer interested in calculating conditional probabilities for
$x_{h}$ and $x_{l}$. But instead, we are interested in calculating $n$
pointwise similarity scores based on smoothed spatial proximity to each other
McInnes, Healy and Melville (2018), e.g.,

$x_{j|i}=exp[(-d(i,j)-\rho_{i})/\sigma_{i}],$ (4.4)

where $\rho$ and $\sigma$ are hyperparameters that control the smoothing
process. Thus, as with t-SNE, there are several hyperparameters to tune, which
apply globally to the model. For example, others include $k$, which is the
size of the neighborhood in the first search region, epochs, which is the
number of times the algorithm sees the data, and also distance, which is a
distance metric (e.g., Euclidean).

To demonstrate the tradeoff in learned high dimensional structure and its
impact on the local projection of that structure, consider the following
demonstration holding the neighborhood size small and fixed at 5, and letting
the times the algorithm sees the data (epochs) vary to be relatively few times
at 20 to many at 500. The intuition is the true structure is more likely to be
learned when the number of epochs is high. The code for this demonstration
relies on the umap package to implement UMAP with the ANES data. Indeed, there
are other packages to fit and diagnose a UMAP model, one of which is covered
at the end of this Section.

Figure 4.3: Comparing UMAP Fits Across Epochs: 500 vs. 20, with $k=5$

From Figure 4.3, we are again able to see clear segmenting on a partisan
basis. But for the purposes of comparing the number of epochs, note the
strange “islands” of observations to the right of the main cloud of data
points in both plots, but especially the right plot in Figure 4.3. Likely, the
reason for this strange pattern is that that algorithm was not given enough
time to learn the true structure of the manifold, such that it had a harder
time of representing it in lower dimensional space. To verify this, and
similar to the grid search exercise with t-SNE, we can update the code to
conduct a grid search across several values of $k$ and the number of epochs. I
essentially update the t-SNE grid search code, and store results after several
fits of UMAP. Then, results are plotted in Figure 4.4, with values of $k$ and
columns showing varying number of epochs.

Figure 4.4: UMAP Grid Search Across $k$ and Epochs

Of note, the UMAP grid search took only about 5.5 minutes to complete, which
is about 3.5 times faster than the t-SNE algorithm on the same data. This
marked increase in computational efficiency is one of the biggest improvements
over t-SNE.

The results in Figure 4.4 confirm our expectations that the strange island of
observations off to the left originally in Figure 4.3 are indeed due to
insufficient time for the algorithm to learn the true structure of the latent
manifold. We see those points slowly migrate over to the full cloud of
observations as the number of epochs increases (i.e., moving from top to
bottom) across Figure 4.4. Given these patterns, though no formal rule exists
for selecting optimal hyperparameter values, I recommend adjusting values
until patterns stabilize. Too much tuning and we are in threat of overfitting
to the training/sampled data. Yet, too little tuning, and we run the risk of
learning or picking up on noise in the data. The trick is to find the sweet
spot across hyperparameter values, such as when the islands in Figure 4.3 stop
shifting as the number of epochs increases for example. Thus, the grid search
is immensely helpful in this regard by providing several visualizations in a
single plot to be able to observe shifts back-to-back across different
hyperparameter values.

In sum, the results across all techniques substantively suggest that members
of all parties, while not extremely separated in the projection space, are
indeed consistently distinct from each other, with minimal overlap between
Democratic respondents compared to Non-Democrats.

###  4.4 A Tidy Tangent on UMAP via uwot::umap() and embed

For the sake of thoroughness, consider a very different, but tidy approach to
UMAP. To do so, I use the tidymodels package, which requires use of two other
packages: uwot and embed, instead of the formerlly used umap package. Here, I
replicate the initial UMAP fit in the left plot of Figure 4.3, where $k=5$ and
the number of epochs $=500$.

The main value, and thus justification of this tidy tangent, is to demonstrate
that different packages and approaches can be used for the same task. This may
come in handy when certain programming grammar is more comfortable to certain
programmers, compared to others. In this case, I use the tidyverse approach,
syntax, and packages (with the exception of uwot). As noted at the outset of
the Element, tidy programming is built on the principle of adaptable code
meant first and foremost to be readable by humans. Related is the ability to
pipe/stack multiple functions, even plotting functions, which results in a
single chunk of code to prepare data, fit a model, and visualize the results.
Such an approach is ultimately simpler and much more efficient.

A second order benefit of this tangent is to demonstrate that, though some
packages differ in programming philosophy, the results are stable across
different constructions. I encourage readers to inspect package documentation
for umap and uwot (for the tidy implementation) packages to understand how
they differ. With that, consider the following code that uses similar syntax
as in the feature engineering code in Section 2 of the Element, to set up the
embedding process via UMAP. Then, I pipe the projection to a ggplot and get a
virtually identical rendering of the projection subspace in Figure 4.5.

Figure 4.5: A Tidy UMAP Solution

The take away from Figure 4.5 is that, though placement and scales are
different compared to the umap package approach, the algorithm still learns
the underlying manifold along which parties are clearly separated.
Substantively, given UMAP’s nonparametric and unsupervised architecture, this
offers an additional point of validation of results across several techniques
that these ANES preference data are relatively stable and reasonably well-
separated along a partisan dimension. A possible take away could be that the
many dimensions used in these types of surveys might be redundant, as we are
able to pick up on a lot of information just by looking to party affiliation
of the respondent. Of course, a generalization of this sort is not possible
from these data alone, but the possibility exists and deserves a closer look.

###  4.5 Suggestions for Further Reading

Wattenberg, Viégas and Johnson (2016) offer a widely used, practical guide to
t-SNE, which gives a nice complement to the coverage of t-SNE in this Element.
Further, Waggoner (2021) presents an application of UMAP in American
Congressional policymaking related to understanding COVID-19 policy structure.
Ordun, Purushotham and Raff (2020) directly compare UMAP to topic models and
digraphs using COVID-19 Twitter data. These two applications offer a “real-
world” view of these methods and their value in practice.

##  5 Neural Network-Based Approaches

In this final section, we shift gears to a very different approach to
dimension reduction at a conceptual level, though technically we are still
interested in some weighted representation and transformation of input
features to give the low dimensional representation of the high dimensional
data space. The conceptual shift is in how to think about the data first, and
then how to conduct the dimension reduction exercise. That is, we now consider
dimension reduction as a neural map. Within the neural network-based
approaches, which are still unsupervised, there is further nuance, mostly
centering on the number of hidden layers. If we have zero hidden layers, for
example, we are in the world of self-organizing maps. But if we have 10 hidden
layers, we are in the world of deep autoencoders. But before we get too far
down the path, we will start with a high-level look at a neural network in
order to place unsupervised dimension reduction in a neural network-based
context.

###  5.1 A Basic Neural Network Architecture

Neural networks are widely used computational modeling tools, which can give
extremely accurate predictions. Though mostly used for supervised learning
tasks, their architecture (stacking layers for processing tasks) can be
adapted to an unsupervised framework for dimension reduction. But to
understand their construction in such a light, an understanding of the basic
architecture is required. As a caveat, this section is only interested in
presenting the very basic architecture. Indeed, many complications and nuanced
extensions exist. The goal at present, rather, is to introduce unfamiliar
readers to the structure of a neural network and how it processes data to give
predictions in order to understand them in in the dimension reduction context.
For a more exhaustive treatment of neural networks and in a deep learning
framework, see Goodfellow et al. (2016).

The simplest neural network is comprised of three layers: input, hidden,
output. The input layers is the raw data matrix we have seen several times to
this point, $\mathbf{X}$, which is simply an $N\times P$ matrix of raw input
features. The output layer is the response feature/outcome/dependent variable.
This is where predictions of the output are produced. The hidden layer is
sandwiched between the input and output layer and is where the processing of
data occurs. Before getting to processing of data, it’s important to note that
a neural network is connected by weights and biases.555The term bias here is
not the same as in supervised prediction tasks referring to mistakes. A weight
is the value attached to the edges in the network connected the nodes that
comprise each layer, and the biases are another type of weights that each node
is assigned. Upon initialization, all biases are randomly assigned, amounting
to a less-than-intelligent brain, akin to making a random guess about some
outcome before seeing and thinking about any data or patterns therein. The
task of the neural network, then, is to pass raw data through the network one
or a few observations (called “batches”) at a time. The network processes the
data in the hidden layer on the forward pass through the network, and then
renders a “decision.” This prediction is compared to the labeled outcome
feature values and the error is recorded. The data are then passed backward
through the network through a process called backpropagation. During
backpropagation, weights and biases are iteratively updated in an incremental
way, in reaction to the error recorded after the forward pass (or “forward
propagation) of the data. Then, new batches are passed through the network,
but this time using the new weights and biases learned during the previous
iteration. Error is recorded after forward propagation, and then during
backpropagation weights and biases are updated once more. This process is
repeated until the network has seen all of the data once. Seeing the full data
one time is called an epoch.666We encountered this term in the previous
section on UMAP. The definition is the same: the number of times the algorithm
sees the full data set So, there can be multiple iterations in a single epoch.
We then have the choice of continuing to learn until we reach some threshold,
e.g., reaching a prespecified number of epochs (via tuning a hyperparameter)
or once some acceptable level of error is reached (similar to the stopping
criterion in k-means clustering). Importantly, there are many decisions to
make when training a neural network (e.g., how large should the updating steps
be in backpropagation, where smaller steps mean more precise learning but
potential failing to converge and larger steps mean quicker convergence, but
potential a poorly performing model that fails to generalize well. This
decision is controlled by the learning rate hyperparameter, often referred to
as $\eta$). There are many flavors of neural networks, such as a fully
connected artificial feedforward neural network (this is the basis type
previous described) or a convolutional neural network (most often used for
image recognition). And further, a key decision which may place you in an
entirely different subfield is how many hidden layers should you include to
process and learn from the data? As the hidden layers increase, the depth of
the network also increases, eventually placing you in the deep learning world.
In short, the point is neural networks can get extremely complex really
quickly.

###  5.2 Self-Organizing Maps

Self-organizing maps (SOM) were first introduced in Kohonen (1982). Based
originally in topological data analysis, and closely related to clustering
techniques, SOM have not always been thought of as a special case of a neural
network (or equivalently, constrained k-means clustering), as neural networks
have only recently become known by such terminology. In brief, SOM can be
thought of as a special case of a neural network, where we are interested in
mapping an input layer to the output layer. But the special case is defined by
the lack of a hidden layer to process the data in the SOM. Rather, the input
layer is directly connected to the output layer through a series of weighted
edges, which are iteratively updated as the structure is learned from the
input layer. Let’s pull this apart a bit more.

As with a basic neural network, the input layer consists of the raw input
features. The output layer in a SOM is a fixed lattice or grid with nodes set
by the researcher a priori. For example, a SOM might be mapped onto a
$10\times 10$ lattice, where there are 100 total nodes in the output layer.
The goal is to place observations into nodes, such that observations closer to
each other in high dimensional space are represented closer together in the
lower dimension map. For example, there should be spatially similar
observations placed together in a single node or in a very close, neighboring
node. The SOM algorithm is interested, then, in learning the patterns of
similarity across observations in the high dimensional space and mapping them
based on these similarities in a smooth, lower dimensional subspace, which is
the fixed lattice. We will come back to this more in a bit.

Central to organization of the map is learning. Learning in SOM is based on
the idea of competitive learning, where these fixed nodes are, in effect,
competing with each other to represent some portion of the data in the two-
dimensional output layer/lattice. This type of learning is in contrast to more
classic approach to learning based on error-minimization as we have seen with
many of the techniques to this point. Importantly, it’s useful to point out
similarity to several other techniques covered and mentioned. First, SOM is
neural network-based in that we are placing data in layers in an attempt to
map one layer onto another layer. But SOM is also constrained k-means as
mentioned above, where similar observations are clustering together. Yet,
further still, by projecting the high dimensional input layer onto a two-
dimensional, fixed lattice for the purpose of mapping, we are also engaging in
dimension reduction. Thus, SOM can be thought of as a combination of neural
networks, clustering, and dimension reduction wrapped into a single technique.
This will be useful later in the section at different points.

Put simply, a SOM learns and gives a solution based on structuring the nodes
in the output layers into clusters of nodes, where clusters are formed based
on spatial proximity and similarities across observations. As with a neural
network, the SOM algorithm is initialized with random weights connecting the
input layer to the output. Then, the weighted raw feature values are projected
onto the lower dimensional lattice, and are in the truest sense, mapped onto
this space. Then, through a series of iterations, the weights are updated
based on spatial similarity. The result then is some nodes are similar to each
other, relative to other nodes that are more different and spatially distant.
As learning is competitive, and as competition amongst nodes to represent some
large chunk of the data space is what drives the algorithm, the “winning” node
that represents the bulk of the observations is called the best matching unit
(BMU). The BMU defines the center of the neighborhood of a cluster of similar
observations. There can be multiple BMUs in a single grid, suggesting multiple
clusters within a common data space. For our example using the ANES data and
searching for latent partisan separation in the projection space, we might
expect roughly two, possibly three clusters surrounding two or three BMUs in
the output layer. We will come back to this in the application.

####  5.2.1 Defining the Steps of the Algorithm

The configuration of the topology of a SOM output layer is defined by
iterating across three main steps. The first step is competition, where nodes
neurons compete to represent the input space as previously mentioned. The
smallest distance between the input values and connection weights gives the
BMU. The second step is cooperation, where the winning node becomes the center
of the neighborhood of nodes in the output layer. This step of cooperation is
critical to SOM learning and representation being treated as a neural network,
based on something called Hebb’s Rule, which most simply is, neurons that fire
together are wired together (Hebb, 1949). Hebb’s Rule, which postulates
learning in the human brain based on synaptic connections between external
stimuli and information processing in the neuron, results in smooth learning
within some neighborhood of neurons. That is, neurons in the human brain do
not learn or are stimulated in isolation, but rather do so in small
neighborhoods. There may be a single, central neuron that is most stimulated
by the raw input (e.g., the BMU in SOM), but there is also a neighborhood of
surrounding neurons that reflect the same stimulation, only to smaller
degrees. This notion of smooth learning within neurons is modeled in the
second step of the SOM algorithm. Where the BMU is the representative, winning
node representing the most similarity of all other nodes for some part of the
data, the surrounding nodes are also representing the same similarity, by
virtue of their proximity to the BMU, only to a lesser degree. This process is
controlled by a (usually Gaussian) decay function, e.g.,

$d_{ij}=e^{(\frac{-d_{ij}^{2}}{2\sigma^{2}})},$ (5.1)

where $d_{ij}$, between the winning neuron, $i$, and the neighbor unit, $j$.
$\sigma$ defines the size of the neighborhood. Of note, $\sigma$ is a
hyperparameter set by the researcher, which is often referred to as $R$ for
“radius,” given that we are defining a neighborhood or region around the BMU.
We start with a large amount of learning at the outset to pick up on clear
patterns, and slowly refine the learning as the iterations increase. That is,
the learning rate and amount of cooperation with neighborhoods of nodes
exponentially decays as the patterns of clusters become clearer, e.g.,
$e^{\frac{-M}{V}}$, where $M$ is the number of iterations, and $V$ is a
constant.

The third and final step is learning, where nodes in the neighborhood of the
BMU participate and update together, such that the higher the weights, the
greater the chance the BMU continues to be the BMU in subsequent iterations.
Put differently, based on the previous two steps of competition and
cooperation, the patterns in the data, and the local groups existing across
the lattice become increasingly clearer as the algorithm sees the data
multiple times. The seeing of the data and updating is indeed learning.
Recall, as we are interested in spatial similarities, which defines a well-
fitting solution, the scoring function is calculated as the distance between
measured values from the input layer to their projected positions in the
output layer.

To recap, the first step of competition means that for each output neuron,
$j$, we compute the distance between each weight vector and input vector,
$d(\boldsymbol{W},\boldsymbol{X})$. We find the BMU, $i$, as the node that
minimizes this distance over all output nodes. The next step of cooperation is
defined by identifying all nodes in the output layer, $j$, within the
neighborhood of $i$, defined by the neighborhood size, $R$ or $\sigma$.
Finally, in the third step, for all nodes in the neighborhood of $i$, we
update the connection weights, $w_{*}\in\boldsymbol{W}$, accordingly,

$w_{new}=w_{old}+\eta(\boldsymbol{X}-w_{old}),$ (5.2)

where $\eta$ is the learning rate, $0<\eta<1$. The algorithms stops when some
criterion is met, which is usually when either the weights stabilize, or
cluster configurations are unchanged. The result is a lower dimensional (that
is, 2D) representation of the high dimensional input space.

####  5.2.2 Applying SOM to the ANES Data

To begin the application of SOM using the ANES data, we need to scale that
data, as we are working with weighted versions of raw input features as with
PCA and other similar techniques in this Element. Upon scaling, we set up the
output layer, which recall is a two-dimensional lattice. For our case we will
set it up to be $10\times 10$, and rectangular (instead of hexagonal). For
most of the application, we will use the kohonen package. For some of the
visual and data organization features, we will use core tidyverse packages,
e.g., ggplot2 and dplyr.

With the grid constructed, and neighborhoods within the grid as Gaussian for
the second cooperation step previous discussed, we are ready to fit the
algorithm. Upon passing the input data and the output grid we previously
defined in search_grid, there are a few additional hyperparameters we need to
set prior to running: alpha (learning rate, which exponentially decays as
previous discussed, from 0.1 to 0.001), radius (neighborhood size around
neurons), rlen (number of iterations), and dist.fcts (distance metric). These
hyperparameters apply globally to the model, and will constrain the training
process. The total training process took about 13 seconds to locally run.

Once fit, there are a variety of approaches to inspecting output from a SOM. A
good place to start is to visualize the training progress. That is, how long
did it take over the iterations, rlen, for the weights to stabilize, and thus
the error to flatten out? We can inspect this by calling the value, $changes
from our som_fit object, shown in Figure 5.1.

Figure 5.1: Training Progress for the SOM Fit

In Figure 5.1 we can see that the algorithm stabilized around 350 iterations
or so, of 500 that we set. Importantly, several choices can impact the
performance of the algorithm. For example, the size of the grid, where smaller
grids offer less space to map the data to, giving the possibility of a poor
solution unable to detect separation. Alternatively, setting rlen too small
will not allow the algorithm time to sufficiently learn the latent patterns in
the data, and thus insufficient time to stabilize and find an optimal
solution. These and all other choices need to be made carefully.

We can plot SOM results in many other ways such as summed distances between
neighbors in the grid space (plot(..., type $=$ "dist.neighbours") or the mean
distances between units and observations (plot(..., type $=$ "quality").
Though these diagnostic plots are useful, for the sake of space I point
readers to the documentation (e.g., ?plot.kohonen), and focus instead on tying
results back to our substantive example of understanding latent partisan
structure in feeling thermometers. To this end, I will focus on the codes
produced by the trained SOM.

Codes, which are sometimes called weight vectors (called via $codes), show the
representation of each feature in each node. These, which are akin to feature
loadings in PCA, give us a clear understanding of the organization of the SOM
solution in the output layer. Some features will contribute to different
locations on the grid, compared to other features. By inspecting the code,
then, we can better understand which features are contributing to the
segmentation of the space. The standard output from a SOM fit in R is
something called a fan plot, where the size of each fan blade, one for each
feature in each node, indicates magnitude (bigger = greater magnitude). Thus,
there are individual fan plots in each node in the output from a basic call of
the plot() function. Yet, this plot is less intuitive and does not fully
describe the value of the SOM solution for helping toward our substantive goal
of understanding whether separation along a partisan dimension characterizes
these feeling thermometer data. Therefore, a more informative approach is to
fit a clustering algorithm to the codes output from a SOM. The intuition here
is to explore whether codes are naturally grouped along a substantive
dimension. For our case, of course, this is along a partisan dimension.

To do this, we require a few steps. First, for ease of plotting, store the
point colors used throughout in two objects: point_colors and neuron_colors,
with the former darker shades corresponding with the individual observations
and the latter lighter shades corresponding with the nodes/neurons in the
output layer (the background in the grid shown in the Figures below). Then, I
fit a k-means algorithm to the codes data, searching for two clusters given
the dichotomous party affiliation feature used throughout the Element. Then,
derive cluster labels (1 or 2) for the clusters, and plot accordingly varying
color by party affiliation. The expectation is that the majority of
observations (points) should correspond to the color of the nodes found from
the k-means clustering solution. If points are scattered and not clearly
coordinating, then this would suggest there is not clear separation in the
output layer.

Figure 5.2: Searching for 2 Clusters via K-Means

Indeed, Figure 5.2 shows clear separation in the output layer along a partisan
dimension. Red points (Non-Democratic) are mostly grouped on the red side of
the output layer, compared to blue points (Democratic) being mostly groups on
the blue side of the output layer. There is not perfect separation here but
based on the results throughout this Element showing separation between
parties, but still some blending across parties, these results corroborate
this pattern.

But note, k-means clustering is more of a brute force approach to clustering,
where an observation is assigned to one and only one cluster and must be
assigned to a cluster (rather than left unclustered). But given the blending
of several observations across parties we have seen to this point, such a hard
partitional approach to clustering may not be the best approach to clustering.
Rather, a soft clustering approach like fuzzy c-means clustering, which is
based on a majority decision from a fractional assignment of observations
theoretically belonging to both clusters, just to varying degrees, may be more
appropriate. Note, for readers unclear on these approaches to clustering in
the context of social science problems, I point them to the recent Element in
this series, Waggoner (2020), which details these two and several other
algorithms for partitioning some data space in an unsupervised way. We can now
update our solution, but this time using fuzzy c-means via the fcm() function
from the ppclust package. See the results in Figure 5.3.

Figure 5.3: Searching for 2 Clusters via Fuzzy C-Means

The results in Figure 5.3 look slightly more like the patterns seen elsewhere,
though they are not extremely different from the patterns in Figure 5.2
showing the k-means version. Separation for most of the observations is clear,
and there are some overlapping cases in both the k-means and fuzzy c-means
(FCM) versions. Even still, FCM is more theoretically motivated given the
potential of observations in the middle of the two dimensions plausibly being
“assigned” to either party. Thus, the FCM approach would be more appropriate
in this case.

Substantively, codes also provide a sense of the specific features that are
attracted to similar nodes, such that we are able to pick up on grouping
across features in the output layer. For example, similar features (e.g.,
feelings toward Barack Obama and Bernie Sanders) might heavily characterize
(or “weight”) a node based on similarity, compared to other features (e.g.,
feelings toward Donald Trump). The result is a picture of the relationships
between both features and nodes across the full input space.

Given the high dimensional space with which we are working, I demonstrate this
view of feature-level code explorations by a sampling of a feature-by-feature
basis and tie the results back into a substantive understanding and
exploration of the feeling thermometer space from the ANES. Such a use of
codes from our fit SOM results in a slightly nuanced view of these weight
vectors compared to more common interpretation previously discussed.
Substantively, we can think of weight vectors/codes like we might correlations
between features and nodes. Features that are more similar will trend toward
nodes in a positive direction, suggesting similar grouping and thus similarity
in latent structure across those features. This is compared to features that
are more different from each other, which will trend negatively. No
relationship would result in a uniformly distributed cloud of observations
with little to no slope.

Consider the direction of the codes between feelings toward Trump and Obama.
We might expect a negative relationship across these features. Inversely, we
can see that features that we might expect to be more similar to each other or
picking up common structure to be trending in a positive direction. For
example, we can see this positive trend in feelings toward Bernie Sanders and
Barack Obama. See both of these examples inf Figures 5.5 and 5.5 with an
overlaid loess smoother to aid visualization of the trend.

Figure 5.4: SOM Codes: Trump and Obama

Figure 5.5: SOM Codes: Sanders and Obama

These simpler expectations are clearly negative and positive as we might
expect across political actors in Figures 5.5 and 5.5, respectively. Though
more overt, the same logic holds for base expectations across institutions, as
we might expect negative relationships between feelings toward the UN and the
NRA, but positive between ICE and the NRA. As a final “sanity” check, this is
exactly what we see in Figures 5.7 and 5.7, suggesting the SOM is picking up
this latent structure in line with substantive expectations.

Figure 5.6: SOM Codes: UN and NRA

Figure 5.7: SOM Codes: ICE and NRA

In sum, the task has remained constant, which is to learn patterns in the data
and then project these patterns in a lower dimensional subspace. Yet, the
views of SOM as a blend of clustering, neural network, and dimension
reduction, mean that simpler interpretation is not as obvious as in PCA, for
example, where we clearly get a lower dimensional solution with PCA scores for
the newly calculated PCs/features. Rather in SOM, we are certainly reducing
the dimensionality of the input space, but we are doing so in a neural
network-based way, and of equal importance to dimension reduction in a typical
SOM application is the clustering of observations in the output layer. This
potential downside in interpretation is outweighed by the extremely
flexibility and wider value from the solution, compared to more nuanced
dimension reduction techniques covered to this point.

It is worth concluding with a word of caution when using SOM. Like t-SNE, the
solutions will never be identical, as the process of weight-updating and
learning is stochastic. Thus, different implementation of SOM can potentially
give different solutions, which limits generalizability and applicability of
this method. Further, SOM are entirely dependent on the shape and size of the
output layer, which recall is a grid/lattice pre-set by the researcher. Set
too small, and no structure can be revealed. Set too large, and parts of the
grid will be left unpopulated by observations, potentially implying much
greater separation than may truly exist in the data. Thus, selecting the
“right” grid size is not always an apparent choice, and requires careful
validation and transparent reporting. I recommend a similar approach when
reporting SOM results as the Risk Averse Workflow discussed at the conclusion
of the LLE section. Still, though, if situated in a broader context of
dimension reduction and compared to several other techniques, SOM can be a
very powerful tool for nonparametrically reducing the complexity of a high
dimensional data space.

###  5.3 Autoencoders

We come now to the final technique: autoencoders. Autoencoders are built on a
similar intuition as PCA, and several other techniques covered to this point,
but all from a neural network-based framework. In brief, autoencoders taken an
input data matrix, create an encoded version of that data matrix by forcing
information loss (much like in PCA), and then in the final stage they attempt
to reconstruct the original input space, but only on the restricted, encoded
version. The goal, then, is to create a reduced version of the input space and
then essentially assess the quality of that encoded version by testing its
ability to inform a good reconstruction of the original high dimensional input
space.

A typical autoencoder topology, then, consists of three layers, like the basic
feedforward neural network we discussed at the outset of this section: input
layer, hidden layer, and an output layer. As before, the input layer is the
matrix of raw data inputs including all features, which in our case is the set
of feeling thermometers from the ANES. The hidden layer is also where data is
processed. But in an autoencoder architecture, “processing” information
equates to compressing the high dimensional data space to be on a lower
dimensional subspace. That is, the hidden layer is of a size less than the
size of the input layer, such that a bottleneck essentially forces the input
layer to be at least of size $p-1$, though in practice the hidden layer is
much smaller than the input layer. This process of data compression in the
hidden layer results in a new, lower dimensional version of the data, as we
saw from PCA, for example. But the output layer is the same size as the input
layer, where the task of the autoencoder in the latter half of the network is
to reconstruct the original input space on the basis of the compressed version
emerging from the hidden layer.

Thus, there are two steps to a basic implementation of an autoencoder: first,
encode; second, decode. The encoding step is the data compression/bottleneck
step. The second decoding step is reconstruction step. As with other neural
networks, once the batch of data has made its way through the full network,
error is calculated prior to backpropagation. Yet, at this point readers might
be wondering: I thought this was unsupervised, so how do we conceptualize and
calculate error? To address this, think about the decoding step; it is
concerned with reconstructing the input layer. Thus, the target is a perfect
replication of the input layer. So the input layer acts as our “labeled data,”
thereby giving a ground truth to which we are able to compare the decoded
version of the data space. This is not a ground truth as we might consider in
a typical supervised task like classification, where the ground truth is a
labeled class, $\\{0,1\\}$ or $\\{\textrm{no},\textrm{yes}\\}$, or
$\\{\textrm{Republican},\textrm{Democrat}\\}$. Rather, the “ground truth” for
an autoencoder is simply the original input space. In this way, we are able to
calculate error (comparing the decoded version to the original version), which
is called reconstruction error, while still remaining firmly in an
unsupervised framework. Put simply, then, the task of the autoencoder is quite
simply to learn effective patterns or structure underlying some input set of
data. This is exactly the goal of unsupervised dimension reduction as we have
addressed throughout the Element. The main difference with autoencoders
compared to other dimension reduction techniques, then, is that it is rooted
in a neural network-based approach to computational modeling (different from
LLE or PCA, e.g.), but still very different from a feedforward neural network
for a supervised problem (e.g., regression or classification).

To reiterate, we are interested in training the network to minimize
reconstruction error. By backpropagating through the network and iteratively
updating weights based on learned structure through the reconstruction error,
the autoencoder is, in a very true sense, learning from the data in an
unsupervised way.

Vastly important to understand with autoencoders is the role of the hidden
layer. As previously noted, the hidden layer forces information loss at the
encoding stage, and thus constrains the picture of the high dimensional data
the decoder gets to see. By creating a lower dimensional representation and
then forcing the reconstruction of the original high dimensional input space,
the resultant bottleneck is intentionally making the decoder work to learn
underlying non-random structure in the data. Indeed, if there were loss of
information, then the task of the decoder would be to simply multiply the
input space by 1, which is a useless task, as there would be nothing to
decode.

From an autoencoder, as with PCA, we get a set of newly constructed features.
Though called principal components in PCA, these are now called either codings
or deep features in autoencoders, depending on the depth of the network. And
by depth, I mean increasing the number of hidden layers, or opportunities to
process and learn from the data. The more hidden layers, usually the better
the solution. However, too many hidden layers, and we are in danger of
overfitting, or learning too much nuance from the training set, such that we
are left with poorly generalizable solution. Still, deep autoencoders of this
sort can be extremely powerful.

A final point of clarification is the number of nodes or neurons in the hidden
layer. If there are more of these than in the input layer, then this would be
considered an overcomplete autoencoder. This is compared to an undercomplete
autoencoder, where the number of nodes in the hidden layer is less than the
size of the input layer, which in our case is 35. So at a minimum, as noted
above, the size of the hidden layer should be $p-1$. There are much fewer uses
of an overcomplete autoencoder compared to the undercomplete version, because
of the very nature of the use of these in practice. That is, if our goal is to
force information loss, and thereby compress the data in the hidden layer to
make the task of the decoder harder to learn better, then it follows that the
size of the hidden layer should indeed be smaller to create this information
loss. Alternatively, in the overcomplete version of an autoencoder, some nodes
may be virtually empty or non-representative, which amounts to cluttering the
encoding step. Though this could be a justifiable decision to make the task of
the decoder more difficult in a different way, the loss of information is
consistently considered the core of the logic of the autoencoder. Thus, I
recommend undercomplete autoencoders in most applications, unless
justification otherwise is extremely well-motivated.

Extending common notation (e.g., Goodfellow et al. (2016)), consider a brief
formalization of the generic form for the encoding step,

$M=f(\boldsymbol{X}),$ (5.3)

where $M$ records the codings from the original input, $\boldsymbol{X}$. Then,
the decoder’s goal is to reconstruct the input, $\boldsymbol{X}$, based on the
codings, or the “representation layer” from 5.3. We then define the decoder
as,

$\boldsymbol{X^{\prime}}=f(M),$ (5.4)

where $\boldsymbol{X^{\prime}}$ is the reconstruction of the inputs, based on
the $M$ from 5.3. To compare $\boldsymbol{X}$ and $\boldsymbol{X^{\prime}}$,
recall we must record the error from the reconstruction attempt. Thus, the
goal is to minimize the reconstruction error, in an attempt to get as close as
possible to the original input space, but on the basis of the codings from the
representation layer,

$\mathrm{min}\mathcal{L}=f(\boldsymbol{X},\boldsymbol{X^{\prime}}),$ (5.5)

or more compactly,

$\mathcal{L}(\hat{\boldsymbol{X}},\boldsymbol{X}).$ (5.6)

####  5.3.1 Applying Autoencoders to the ANES Data

For this final application using the 2019 ANES data, I will first walk through
constructing an autoencoder using the machine learning engine, H20. But,
extending beyond past applications in this Element, I will then use the
extracted features (“deep features”) in a simple supervised task to
demonstrate the value and ease of such an extension. This application, then,
is intended to give a realistic, albeit limited, picture of a common use of an
autoencoder in a social science task.

We start with setting a few things up. In the code for Section 6, we start by
initializing the H20 session, set our party affiliation feature aside for
coloring as usual, but also for our supervised task at the end of the section.
Also, to speed up processing time, as these types of neural networks can take
a long time, especially as they and the data grow increasingly complex.
Finally, I divide the data into training, testing, and validation sets. This
is a common approach in machine learning research to train, test and tune
learners, but not threaten the generalizability of the learner. For more, see
Friedman, Hastie and Tibshirani (2001). Now, with the data and H20 session set
up, we are ready to build a shallow (non-deep) undercomplete autoencoder.

The syntax to build an autoencoder using the H20 engine includes several
tunable hyperparameters and several default values. Of note is the autoencoder
argument, which when set to TRUE allows us to build an autoencoder. Also, we
have specified that we want a single hidden layer with 16 nodes, per
discussion above on ensuring this is an undercomplete autoencoder. Once
trained on 60% of out data, we are also telling the h2o.deeplearning()
function that we want to “test,” generate predictions using the held out 20%
of the data. This is the argument that allows for calculation of
reconstruction error. Once run, we need to extract the deep features/codings
learned from the model. That is, we want to focus less on reconstruction at
this stage, and instead extract the codings from the representation/hidden
layer. This will help generate interpretable and comparable plots in line with
the other techniques and approach throughout the Element. To do so, we rely on
the h2o.deepfeatures() function to extract the deep features.

The returned output is a bit odd. The deep feature labels, e.g., DF.L1.C1 are
read: “data frame, layer number, column number.” For our purposes, which is a
simple application, we only a single data frame and a single hidden layer
(i.e., layer = 1 in the previous function call). Thus, the only value in the
titles that will change is the $C_{*}$. As we specified a compressed hidden
layer with 16 nodes/neurons, then the output returned 16 column vectors,
giving a deep feature matrix of size $1910\times 16$, which is the size of the
training set. For substantive purposes, we can plot these deep features
against each other and color by party affiliation as before to get a first
look as to whether they are picking up on partisan separation in the
projection space. Due to limited space, the code to produce these plots is
omitted, but can easily be built by repurposing some of the earlier plot code
in the Element. The output is in Figure 5.8.

Figure 5.8: Deep Feature Plots

Surprisingly, the patterns suggest less separation than prior plots from other
techniques we have seen, though separation across all deep features is clear.
Democrats are indeed distinct in their feelings toward a battery of issues and
actors compared to Non-Democrats. Importantly, these patterns are only from
part (60%) of the data, and could be the reason for the different patterns. I
encourage readers to plot across the validation set to compare, or even change
the portions of data splits in the sets.

But to demonstrate a fuller application, we will turn now to a supervised
task. That is, rather than color points by party affiliation, we will predict
party affiliation as a function of the deep features extracted from our
autoencoder. To do so, we first extract the features, and then use these as
input in a deep neural network. Importantly, the model at present now is not
unsupervised as we have covered to this point in this Element. But it is
supervised, where we have labeled output, party affiliation. We are attempting
to classify respondents’ party affiliations based on the deep features we
extracted from our autoencoder fit above. So there are two neural networks at
play: the first (autoencoder) giving us our set of input features, and the
second (deep feedforward neural network based on the code below, hidden = c(8,
8) denoting two hidden layers, each with 8 neurons) for predicting party
affiliation as a function the deep features. A common and simple way to assess
the quality of a binary choice model of this sort is to generate a confusion
matrix showing the true (false) positive predictions (Democrat) relative to
the true (false) negative predictions (Non-Democrat). These results are shown
in Table 5.1.

Table 5.1: Confusion Matrix For Deep Neural Network Predicting Party Affiliation | Observed
---|---
Predicted |  | D | ND |
D′ |  210 83.7% |  61 16.7% |
ND′ |  41 16.3% |  304 83.3% |

The confusion matrix in Table 5.1 suggests our network did a very good job at
predicting party affiliation as a function the deep features from the
autoencoder based on the “true” rate reported in the main diagonal of the
matrix. This also suggests that the autoencoder found deep features that are
picking up on substantive, latent structure given their ability to predict
party affiliation at such high rates (i.e., high true positive and true
negative rates, relative to false versions of both).

A final approach to explore results from an autoencoder is exploring feature
importance. Feature importance is often measure by a feature’s contribution to
overall fit of the model. For our purposes, we explore feature importance from
the deep network that was used for the prediction task. Substantively, this
will help by telling use which deep feature from the autoencoder was driving
the prediction solution. There are two main ways to think about feature
importance: relative importance (summing to 1.0 across all features) and raw
percentage (contributed by a given feature). I explore this and present the
results for each in Figures 5.10 and 5.10.

Figure 5.9: Relative Feature Importance

Figure 5.10: Raw Feature Importance

Across both of these plots of feature importance in Figures 5.10 and 5.10,
deep features 12 and 13 are the most important by a wide margin. Note, this
pair was not a combination of plotted features in Figure 5.8. With this, as a
final inspection of the results, we will now plot these against each other as
they were most influential in the classification task. Thus, we might expect
to see the clearest separation between party affiliation across these two deep
features. I plot these against each other using both the training and
validation sets below in Figure 5.11, as patterns should be similar across
each.

Figure 5.11: Deep Features 12 and 13 by Party Affiliation

Now, the patterns in both training and validation sets in Figure 5.11 show a
pattern that looks quite similar to patterns from other techniques throughout
the Element, where there is very clear separation between the parties in the
projection space, but still some blending toward the middle of both
dimensions. This suggests, then, that deep features 12 and 13 in this case can
be thought of as the first two “principal components” from a PCA fit. Thus,
the output of an autoencoder does not automatically sort the deep features
based on any metric, whether something akin to PVE in PCA, or feature
importance as we have explored here. As such, it is up to the researcher to
dive into the results and effectively pull them apart and explain patterns in
a substantive light.

Importantly, as with all other techniques covered in this Element,
autoencoders are not without limitations. As the task of an autoencoder is to
reproduce the original input space based on a simplified or “encoded” version
of that space, the solution may very easily perform poorly on new data. That
is, by training and tuning an autoencoder to learn defining aspects of the
latent structure to aid in reconstruction of the high dimensional space, it is
in effect formalizing a process of overfitting to the training data. The
training data in this scenario is the full input space, and the decoded output
layer gives a learned version of that space based on the encoded version of
the same space. For our present application working with nationally
representative survey data, we might bypass this weakness by trusting that the
ANES sample gives a realistic look at the American population. Thus, an
“overfit” autoencoder would not be as much of a worry as the solutions should
work relatively well on some other sample drawn from the same population using
the same ANES sampling techniques. Yet, from a methodological perspective, we
may be interested in a technique that gives a more generalizable solution,
especially in cases leveraging worse or different data than a nationally
representative American survey. For example, such an approach to guard against
this threat of overfitting that is innate to all vanilla autoencoders, is a
recent innovation called a denoising autoencoder. The task of the denoising
autoencoder is to first add statistical noise to the original input layer.
This initial step has the effect of making the job of the decoder much harder,
as it must attempt reconstruction based on a noisy version of the already-
reduced input space. And the increased difficulty of decoding in this scenario
should, in theory, allow for a more generalizable solution (Goodfellow et al.,
2016, 504-505). The point here is that though sophisticated and often
accurate, even complex methods like autoencoders have weaknesses and
limitations that must be considered when used in a research project. At a
minimum these limitations should be mentioned in the reporting of the results.

###  5.4 Suggestions for Further Reading

As with the previous section, the value of these methods is clearest when
looking at implementation in applied settings. For example, Kourtit, Nijkamp
and Arribas (2012) use self-organizing maps to explore smart cities, and Li et
al. (2018) derive a clever application of an autoencoder to generate poetry.

##  6 Final Thoughts on Dimension Reduction

It has been my goal in this Element to present a modern picture of dimension
reduction, with application in the social sciences. The modern part of this
goal focusing on recent methodological development, though, is only possible
in light of the simple, but powerful original approach to dimension reduction:
principal components analysis (PCA).

Building on PCA and the shared goal of learning from data to produce a simpler
projection, LLE is based on a similar linear construction, but is much more
flexible to model nonlinearities in data. Yet, t-SNE and UMAP offer very
different approaches to capturing global structure and projecting it locally.
The extensions of t-SNE and UMAP allow for even more flexibility in model
complex structures beyond what either LLE or PCA could handle, and they do so
in a fundamentally different (graph-based) way. Further still, another
approach to dimension reduction is through a neural network architecture
focusing on both grouping patterns in data and also minimizing reconstruction
error. Both SOM and autoencoders, though differing in construction, take
advantage of an iterative process for filtering data to learn and project
structure. Through coverage of these techniques, this Element has given a
framework for thinking about and applying dimension reduction to address a
host of problems all stemming from data complexity and size. This framework
guides method selection, which might be dependent on the goals of the project
(e.g., feature extraction versus representation) or conceptualizations of how
best to think about similarity, structure, and complexity (e.g., maximizing
variance versus manifold reconstruction).

All techniques covered in this Element have been extended over the years in
numerous ways, with deepened complexity in their own rite. That is, we can
always go deeper and get more complex, dependent on the nature and scope of
the problem at hand. For example, there are many versions of autoencoders
beyond even the deep version we briefly touched on in the previous section.
There are variational, sparse, and denoising autoencoders, to name a few. Each
of these and the many other extensions were derived to address unique problems
in extremely clever ways. As noted in the previous section, the denoising
autoencoder is interested in making the decoding task harder to (hopefully)
give a more reliable solution. It does so by adding some extra “noise” to the
original input layer. The task of the decoder, after the encoding stage, is
not only to now reconstruct the encoded input space, but also to separate the
noise from the signal, such that the added noise should be absent from the
decoded solution. Such extensions are worth mentioning to not only demonstrate
the rapidly developing nature of the field, but to also encourage readers to
think critically about the best technique for the problem at hand. And
importantly, if a technique does not exist, then readers should feel empowered
to follow similar past approaches and build on the firm foundations, such as
the autoencoder, to develop the needed extension to solve some problem. To
underscore this point, it is important to note that the autoencoder is itself
an extension of an existing method: the restricted Boltzmann machine (RBM).
The RBM, which was an extension of the Boltzmann machine, was essentially
attempting to leverage a neural network framework to develop a version of PCA,
where the goal was the output from the encoding stage. The RBM algorithm stops
short of decoding the encoded layer, and instead outputs the reduced set of
features, similar to the result we get from a PCA solution. Thus, by adding
the inverse of the encoding stage to the backend of the RBM, the autoencoder
was developed.

Perhaps an even broader, but still related conclusion is that we have only
scratched the surface regarding coverage of techniques. Though much more
ground could be covered in building a dimension reduction toolbox, instead
this Element has introduced a framework for thinking about dimension reduction
in a holistic light. To this end, we have covered dimension reduction on the
basis of calculating and extracting new features on the basis of learned
trends in the data (maximizing variance via PCA). We have also covered
dimension reduction for the purpose of learning and representing the full
latent structure of the input space instead of just a summary in a linear way
(via LLE) and in a nonlinear way, for both feature extraction and
visualization (via t-SNE and UMAP). Further, we have covered dimension
reduction in a way that emulates the process of neurons firing and learning
within neighborhoods on the basis of raw stimuli (via SOM). And finally, we
covered dimension reduction for learning the structure of a space by creating
and solving a problem based on an imperfect version of the original space (via
autoencoders). Though the general goal remained constant, the selection and
implementation of the technique resulted in unique ways to conceptualize and
solve a common problem. And through the power of visualization, we
demonstrated throughout that the patterns uncovered, regardless of variance in
method-specific mechanics, were largely consistent suggesting the latent
structure was indeed real, rather than a function of noise or randomness in
the data. If the latter were the case, then different approaches to solving
the dimension reduction problem would have given different versions of the
same data space.

Related, and perhaps most importantly, coverage and presentation of the
techniques in this Element was based in the recognition that method selection
always flows from a theoretical assumption of the world, if even that
assumption is not formalized. That is, our conceptualization of both the input
data space and also the ideal way to make that input space more
understandable, is manifested in the selected method. As with supervised
modeling tasks, in unsupervised dimension reduction the mere selection of a
technique is itself a theoretical assertion at some level. For example,
deciding whether to create a lower dimensional version of the full space is
based in the assumption that the contours of the space (e.g., UMAP) reveal a
more useful version of the data than a global summary of part of the space
(e.g., PCA). Or alternatively, autoencoders do not rely on the local linearity
assumption of the manifold as LLE does, so the modeling process is not bound
by protecting against a violation of linearity. The result is may be better
performance by the algorithm, but the sacrifice may be in interpretability.
LLE is much simpler to interpret than a deep autoencoder. As a result, there
is an implicit tradeoff in model complexity versus interpretability in
unsupervised learning as there is in supervised learning. There should, at a
minimum, be an appreciation of this reality along with the fact that our
unique biases influence every model we build and every choice we make in the
research process. This is not a normatively “bad” or “good” thing; it is
merely a reality. So, by selecting a model, we are in a sense codifying our
assumptions in tandem with making a judgement on the
complexity/interpretability tradeoff. This codification through model
selection should be appreciated, well-motivated, and justified throughout the
modeling process, whether supervised or unsupervised. As such, I urge readers
to take great care in justifying selection of a technique, and follow the Risk
Averse Workflow addressed at the conclusion of the LLE section.

Along with presentation of this material in a “modern” context means that much
work is currently being done to push the field of dimension reduction forward
through development of new techniques and offering new ways to think about
addressing new problems. Of note is the issue of big data. How best to process
and treat big data in the age-old service of dimension reduction is a
particularly vexing issue (Zhang and Yang, 2018). Some work has been done on
this front. For example, Krishnan, Samaranayake and Jagannathan (2018) develop
a multistep dimension reduction approach to effectively deal with data in
batches. Related, dimension reduction is not beholden to numeric data. A host
of techniques, some entirely different, are used to reduce the dimensionality
and complexity of text data. For example, vector space models like word2vec
can handle massive data of sizes in the billions (Mikolov et al., 2013). The
idea behind word2vec and other vector space embedding models is to reduce
dimensionality by learning and reducing complexity based on semantic
similarity across a text corpus.

In sum, regardless of the type of data, the size of the data, or the
techniques used, dimension reduction offers researchers a powerful set of
tools for making complex data spaces more manageable, interpretable, and
simpler. This Element has shown a few of some of the most widely used, modern
approaches for doing so, along with tips on implementation of the methods in
R. Taken with my earlier Element on clustering (Waggoner, 2020), it is my hope
that the value of unsupervised machine learning for uncovering non-random
structure and so learning from data is clear for a variety of applications
that may have been previously unclear. As a result, readers are encouraged to
dive deeper into this exciting, rapidly developing world and perhaps even
develop algorithms of their own.

## References

  * (1)
  * ANES (2019) ANES, American National Election Studies. 2019.  Pilot Study.
  * Armstrong et al. (2014) Armstrong, David A, Ryan Bakker, Royce Carroll, Christopher Hare, Keith T Poole, Howard Rosenthal et al. 2014.  Analyzing spatial models of choice and judgment with R.  CRC Press.
  * Chang and Yeung (2006) Chang, Hong and Dit-Yan Yeung. 2006.  “Robust locally linear embedding.” Pattern recognition 39(6):1053–1065.
  * Friedman, Hastie and Tibshirani (2001) Friedman, Jerome, Trevor Hastie and Robert Tibshirani. 2001.  The elements of statistical learning.  Springer.
  * Goodfellow et al. (2016) Goodfellow, Ian, Yoshua Bengio, Aaron Courville and Yoshua Bengio. 2016\.  Deep learning.  Massachusetts Institute of Technology Press.
  * Hebb (1949) Hebb, Donald Olding. 1949.  The organization of behavior: a neuropsychological theory.  J. Wiley; Chapman & Hall.
  * James et al. (2013) James, Gareth, Daniela Witten, Trevor Hastie and Robert Tibshirani. 2013\.  An introduction to statistical learning.  Springer.
  * Kennedy and Waggoner (2021) Kennedy, Ryan and Philip D. Waggoner. 2021.  Introduction to R for Social Scientists: A Tidy Programming Approach.  CRC Press.
  * Kohonen (1982) Kohonen, Teuvo. 1982.  “Self-organized formation of topologically correct feature maps.” Biological cybernetics 43(1):59–69.
  * Kourtit, Nijkamp and Arribas (2012) Kourtit, Karima, Peter Nijkamp and Daniel Arribas. 2012.  “Smart cities in perspective–a comparative European study by means of self-organizing maps.” Innovation: The European journal of social science research 25(2):229–246.
  * Krishnan, Samaranayake and Jagannathan (2018) Krishnan, Raghavan, VA Samaranayake and Sarangapani Jagannathan. 2018.  “A multi-step nonlinear dimension-reduction approach with applications to big data.” IEEE Transactions on Knowledge and Data Engineering 31(12):2249–2261.
  * Li et al. (2018) Li, Juntao, Yan Song, Haisong Zhang, Dongmin Chen, Shuming Shi, Dongyan Zhao and Rui Yan. 2018.  Generating classical chinese poems via conditional variational autoencoder and adversarial training. In Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing.  pp. 3890–3900.
  * Maaten and Hinton (2008) Maaten, Laurens van der and Geoffrey Hinton. 2008.  “Visualizing data using t-SNE.” Journal of machine learning research 9(Nov):2579–2605.
  * McInnes, Healy and Melville (2018) McInnes, Leland, John Healy and James Melville. 2018.  “Umap: Uniform manifold approximation and projection for dimension reduction.” arXiv preprint:1802.03426 .
  * Mikolov et al. (2013) Mikolov, Tomas, Kai Chen, Greg Corrado and Jeffrey Dean. 2013.  “Efficient estimation of word representations in vector space.” arXiv preprint arXiv:1301.3781 .
  * Ordun, Purushotham and Raff (2020) Ordun, Catherine, Sanjay Purushotham and Edward Raff. 2020.  “Exploratory analysis of covid-19 tweets using topic modeling, umap, and digraphs.” arXiv preprint arXiv:2005.03082 .
  * Riemann (1873) Riemann, B. 1873.  “On the Hypotheses that Lie at the Bases of Geometry (1854).” English Translation by WK Clifford, Nature 8\.
  * Roweis and Saul (2000) Roweis, Sam T and Lawrence K Saul. 2000.  “Nonlinear dimensionality reduction by locally linear embedding.” science 290(5500):2323–2326.
  * Rubin (1976) Rubin, Donald B. 1976.  “Inference and missing data.” Biometrika 63(3):581–592.
  * Saul and Roweis (2003) Saul, Lawrence K and Sam T Roweis. 2003.  “Think globally, fit locally: unsupervised learning of low dimensional manifolds.” Journal of machine learning research 4(Jun):119–155.
  * Schölkopf, Smola and Müller (1997) Schölkopf, Bernhard, Alexander Smola and Klaus-Robert Müller. 1997\.  Kernel principal component analysis. In International conference on artificial neural networks.  Springer pp. 583–588.
  * Van Der Maaten (2014) Van Der Maaten, Laurens. 2014.  “Accelerating t-SNE using tree-based algorithms.” The Journal of Machine Learning Research 15(1):3221–3245.
  * Waggoner (2019) Waggoner, Philip D. 2019.  “amerika: American Politics-Inspired Color Palette Generator. R package version 0.1.0.”.
  * Waggoner (2020) Waggoner, Philip D. 2020.  Unsupervised Machine Learning for Clustering in Political and Social Research.  Cambridge University Press.
  * Waggoner (2021) Waggoner, Philip D. 2021.  “Pandemic Policymaking.” Journal of Social Computing 2(1):14–26.
  * Wattenberg, Viégas and Johnson (2016) Wattenberg, Martin, Fernanda Viégas and Ian Johnson. 2016.  “How to use t-SNE effectively.” Distill 1(10):e2.
  * Wickham et al. (2019) Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, Alex Hayes, Lionel Henry, Jim Hester et al. 2019.  “Welcome to the Tidyverse.” Journal of Open Source Software 4(43):1686.
  * Zhang and Yang (2018) Zhang, Tonglin and Baijian Yang. 2018.  “Dimension reduction for big data.” Statistics and Its Interface 11(2):295–306.
  * Zou, Hastie and Tibshirani (2006) Zou, Hui, Trevor Hastie and Robert Tibshirani. 2006.  “Sparse principal component analysis.” Journal of computational and graphical statistics 15(2):265–286.