---
title: "Factor analysis"
---

##  Vs. principal components 


* Principal components: 


  * Purely mathematical.

  * Find eigenvalues, eigenvectors of correlation matrix.

  * No testing whether observed components reproducible, or even probability model behind it.


* Factor analysis: 


  * some way towards fixing this (get test of appropriateness)

  * In factor analysis, each variable modelled as: "common factor" (eg. verbal ability) and "specific factor" (left over).

  * Choose the common factors to "best" reproduce pattern seen in correlation matrix.

  * Iterative procedure, different answer from principal components.



##  Packages

```{r bFactor-1, warning=F, message=F}
library(ggbiplot)
library(tidyverse)
library(conflicted)
conflict_prefer("mutate", "dplyr")
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")
conflict_prefer("arrange", "dplyr")
```

 


##  Example


* 145 children given 5 tests, called PARA, SENT, WORD, ADD and DOTS. 3 linguistic tasks (paragraph comprehension, sentence completion  and word meaning), 2 mathematical ones (addition and counting dots).

* Correlation matrix of scores on the tests:


```

para 1     0.722 0.714 0.203 0.095
sent 0.722 1     0.685 0.246 0.181
word 0.714 0.685 1     0.170 0.113
add  0.203 0.246 0.170 1     0.585
dots 0.095 0.181 0.113 0.585 1

```



* Is there small number of underlying "constructs" (unobservable) that explains this pattern of correlations?



##  To start: principal components

Using correlation matrix. Read that first:

```{r kids-scree,message=F}
my_url <- "http://ritsokiguess.site/datafiles/rex2.txt"
kids <- read_delim(my_url, " ")
kids
```


## Principal components on correlation matrix


Turn into R `matrix`, using column `test` as column names:

```{r}
kids %>% 
column_to_rownames("test") %>% 
as.matrix() -> m
```

Principal components:

```{r bFactor-2}
kids.0 <- princomp(covmat = m) 
```

I used `kids.0` here since I want `kids.1` and `kids.2` later.


##  Scree plot
```{r bFactor-3, fig.height=3.5}
ggscreeplot(kids.0)
```

 


##  Principal component results


* Need 2 components. Loadings:

\footnotesize
```{r bFactor-4}
kids.0$loadings
```
\normalsize

## Comments

* First component has a bit of everything, though especially the
first three tests.

* Second component rather more clearly `add` and `dots`.

* No scores, plots since no actual data.

- See how factor analysis compares on these data.


##  Factor analysis


* Specify number of factors first, get solution with exactly
that many factors.

* Includes hypothesis test, need to specify how many children
wrote the tests.

* Works from correlation matrix via `covmat` or actual
data, like `princomp`.

* Introduces extra feature, *rotation*, to make
interpretation of loadings (factor-variable relation) easier.



##  Factor analysis for the kids data


* Create "covariance list" to include number of children who
wrote the tests.

* Feed this into `factanal`, specifying how many factors (2).

- Start with the matrix we made before.

```{r bFactor-5 }
m
ml <- list(cov = m, n.obs = 145)
kids.2 <- factanal(factors = 2, covmat = ml)
```

 



##  Uniquenesses

```{r bFactor-6 }
kids.2$uniquenesses
```


* Uniquenesses say how "unique" a variable is (size of
specific factor). Small
uniqueness means that the variable is summarized by a factor (good).

* Very large uniquenesses are bad; `add`'s uniqueness is largest but not large enough to be worried about.

* Also see "communality" for this idea, where *large* is good and *small* is bad.



##  Loadings

\footnotesize
```{r bFactor-7}
kids.2$loadings
```
\normalsize

* Loadings show how each factor depends on variables. Blanks
indicate "small", less than 0.1.

## Comments

* Factor 1 clearly the "linguistic" tasks, factor 2 clearly
the "mathematical" ones.

* Two factors together explain 68\% of variability (like
regression R-squared).
  
- Which variables belong to which factor is *much* clearer than with principal components.

##  Are 2 factors enough? 
```{r bFactor-8 }
kids.2$STATISTIC
kids.2$dof
kids.2$PVAL
```

 

P-value not small, so 2 factors OK.


##  1 factor

```{r bFactor-9 }
kids.1 <- factanal(factors = 1, covmat = ml)
kids.1$STATISTIC
kids.1$dof
kids.1$PVAL
```

 

1 factor rejected (P-value small). Definitely need more than 1.

## Places rated, again

- Read data, transform, rerun principal components, get biplot: 

```{r bFactor-10, message=FALSE, fig.height=6}
my_url <- "http://ritsokiguess.site/datafiles/places.txt"
places0 <- read_table(my_url)
places0 %>% 
mutate(across(-id, \(x) log(x))) -> places
places %>% select(-id) -> places_numeric
places.1 <- princomp(places_numeric, cor = TRUE)
g <- ggbiplot(places.1, labels = places$id,
       labels.size = 0.8)
```

- This is all exactly as for principal components (nothing new here).

## The biplot

```{r bFactor-11, fig.height=8}
g
```

## Comments

- Most of the criteria are part of components 1 *and* 2. 
- If we can rotate the arrows counterclockwise:
  - economy and crime would point straight up
    - part of component 2 only
  - health and education would point to the right
    - part of component 1 only
- would be easier to see which variables belong to which component.
- Factor analysis includes a rotation to help with interpretation.

## Factor analysis

- Have to pick a number of factors *first*. 
- Do this by running principal components and looking at scree plot.
- In this case, 3 factors seemed good (revisit later):

```{r bFactor-12}
places.3 <- factanal(places_numeric, 3, scores = "r")
```

- There are different ways to get factor scores. These called "regression" scores.

## A bad biplot 

```{r bFactor-13, fig.height=4}
biplot(places.3$scores, places.3$loadings,
  xlabs = places$id)
```


## Comments

- I have to find a way to make a better biplot!
- Some of the variables now point straight up and some straight across (if you look carefully for the red arrows among the black points).
- This should make the factors more interpretable than the components were.

## Factor loadings

\footnotesize
```{r bFactor-14}
places.3$loadings
```
\normalsize

## Comments on loadings

- These are at least somewhat clearer than for the principal components:
- Factor 1: health, education, arts: "well-being"
- Factor 2: housing, transportation, arts (again), recreation: "places to be"
- Factor 3: climate (only): "climate"
- In this analysis, economic factors don't seem to be important.

## Factor scores

- Make a dataframe with the city IDs and factor scores:

```{r bFactor-15}
cbind(id = places$id, places.3$scores) %>% 
as_tibble() -> places_scores
```

- Make percentile ranks again (for checking):

```{r bFactor-16}
places %>% 
mutate(across(-id, \(x) percent_rank(x))) -> places_pr
```

## Highest scores on factor 1, "well-being":

- for the top 4 places:

```{r bFactor-17}
places_scores %>% 
slice_max(Factor1, n = 4)
```

## Check percentile ranks for factor 1

```{r bFactor-18}
places_pr %>% 
select(id, health, educate, arts) %>% 
filter(id %in% c(213, 65, 234, 314))
```

- These are definitely high on the well-being variables.
- City #213 is not so high on education, but is highest of all on the others.

## Highest scores on factor 2, "places to be":

```{r bFactor-19}
places_scores %>% 
slice_max(Factor2, n = 4)
```

## Check percentile ranks for factor 2

```{r bFactor-20}
places_pr %>% 
select(id, housing, trans, arts, recreate) %>% 
filter(id %in% c(318, 12, 168, 44))
```

- These are definitely high on housing and recreation.
- Some are (very) high on transportation, but not so much on arts.
- Could look at more cities to see if #168 being low on arts is a fluke.

## Highest scores on factor 3, "climate":

```{r bFactor-21}
places_scores %>% 
slice_max(Factor3, n = 4)
```

## Check percentile ranks for factor 3

```{r bFactor-22}
places_pr %>% 
select(id, climate) %>% 
filter(id %in% c(227, 218, 269, 270))
```

This is very clear.

## Uniquenesses

- We said earlier that the economy was not part of any of our factors:

```{r}
#| echo: false 
wid <- getOption("width")
options(width = 60)
```


```{r bFactor-23}
places.3$uniquenesses
```

```{r}
#| echo: false 
options(width = wid)
```


- The higher the uniqueness, the less the variable concerned is part of any of our factors (and that  maybe another factor is needed to accommodate it).
- This includes economy and maybe crime.


## Test of significance

We can test whether the three factors that we have is enough, or whether we need more to describe our data:

```{r bFactor-24}
places.3$PVAL
```

- 3 factors are not enough. 
- What would 5 factors look like?

## Five factors

\footnotesize
```{r bFactor-25}
places.5 <- factanal(places_numeric, 5, scores = "r")
places.5$loadings
```
\normalsize

## Comments 1/2

- On (new) 5 factors:
- Factor 1 is health, education, arts: same as factor 1 before.
- Factor 2 is housing, transportation, arts, recreation: as factor 2 before.
- Factor 3 is economy.
- Factor 4 is crime.
- Factor 5 is climate and housing: like factor 3 before.

## Comments 2/2

- The two added factors include the two "missing" variables.
- Is this now enough?

```{r bFactor-26}
places.5$PVAL
```

- No. My guess is that the authors of Places Rated chose their 9 criteria to capture different aspects of what makes a city good or bad to live in, and so it was too much to hope that a small number of factors would come out of these.


##  A bigger example: BEM sex role inventory


* 369 women asked to rate themselves on 60 traits, like "self-reliant" or "shy".

* Rating 1 "never or almost never true of me" to 7 ``always or
almost always true of me''.

* 60 personality traits is a lot. Can we find a smaller number
of factors that capture aspects of personality?

* The whole BEM sex role inventory on next page.



##  The whole inventory


![](bem.png){width=450px}




##  Some of the data


\scriptsize
```{r bem-scree, message=F}
my_url <- "http://ritsokiguess.site/datafiles/factor.txt"
bem <- read_tsv(my_url)
bem
```
\normalsize
 


##  Principal components first
\ldots to decide on number of factors:
```{r bFactor-27 }
bem.pc <- bem %>%
select(-subno) %>%
princomp(cor = T)
```

 


##  The scree plot
```{r genoa,fig.height=3.7}
(g <- ggscreeplot(bem.pc))
```




* No obvious elbow.




##  Zoom in to search for elbow

Possible elbows at 3 (2 factors) and 6 (5):

```{r bem-scree-two,fig.height=3.3,warning=F}
g + scale_x_continuous(limits = c(0, 8))
```





##  but is 2 really good?

```{r bFactor-28, include=FALSE}
options(width = 80)
```

\scriptsize
```{r bFactor-29 }
summary(bem.pc)
```
\normalsize

```{r bFactor-30, include=FALSE}
options(width = 60)
```

##  Comments


* Want overall fraction of variance explained (``cumulative
proportion'') to be reasonably high.

* 2 factors, 28.5\%. Terrible!

* Even 56\% (10 factors) not that good!

* Have to live with that.




##  Biplot

```{r bem-biplot,fig.height=3.5}
ggbiplot(bem.pc, alpha = 0.3)
```




##  Comments


* Ignore individuals for now.

* Most variables point to 1 o'clock or 4 o'clock.

* Suggests factor analysis with rotation will get interpretable
factors (rotate to 12 o'clock and 3 o'clock, for example).

* Try for 2-factor solution (rough interpretation, will be bad):

```{r bFactor-31 }
bem %>%
select(-subno) %>%
factanal(factors = 2) -> bem.2
```



* Show output in pieces (just print `bem.2` to see all of it).



##  Uniquenesses, sorted

\scriptsize
```{r bFactor-32, echo=-1}
options(width = 60)
sort(bem.2$uniquenesses)
```
\normalsize
 
## Comments

* Mostly high or very high (bad).

* Some smaller, eg.: Leadership ability (0.409),
Acts like leader (0.417),
Warm (0.476),
Tender (0.493).

* Smaller uniquenesses captured by one of our two factors.

- Larger uniquenesses are not: need more factors to capture them.



##  Factor loadings, some

\scriptsize
```{r bFactor-33}
bem.2$loadings
```
\normalsize



##  Making a data frame
There are too many to read easily, so make a data frame. A
bit tricky:

\footnotesize
```{r bFactor-34}
bem.2$loadings %>% 
unclass() %>% 
as_tibble() %>% 
mutate(trait = rownames(bem.2$loadings)) -> loadings
loadings %>% slice(1:8)
```
\normalsize
 


##  Pick out the big ones on factor 1

Arbitrarily defining $>0.4$ or $<-0.4$ as "big":

\scriptsize
```{r bFactor-35}
loadings %>% filter(abs(Factor1) > 0.4) 
```
\normalsize



##  Factor 2, the big ones

\footnotesize
```{r bFactor-36}
loadings %>% filter(abs(Factor2) > 0.4)
```
\normalsize
 


##  Plotting the two factors
- A bi-plot, this time with the variables reduced in size. Looking for
unusual individuals.

- Have to run `factanal` again to get factor scores for plotting.

```{r biplot-two-again, eval=F}
bem %>% select(-subno) %>% 
factanal(factors = 2, scores = "r") -> bem.2a
biplot(bem.2a$scores, bem.2a$loadings, cex = c(0.5, 0.5))
```



- Numbers on plot are row numbers of `bem`
data frame.


##  The (awful) biplot

```{r biplot-two-ag,fig.height=4,echo=F}
bem.2a <- factanal(bem[, -1], factors = 2, scores = "r")
biplot(bem.2a$scores, bem.2a$loadings, cex = c(0.5, 0.5))
```




##  Comments


* Variables mostly up ("feminine") and right ("masculine"),
accomplished by rotation.

* Some unusual individuals: 311, 214 (low on factor 2), 366
(high on factor 2),
359, 258
(low on factor 1), 230 (high on factor 1).



##  Individual 366

\tiny
```{r bFactor-37}
bem %>% slice(366) %>% glimpse()
```
\normalsize

## Comments



* Individual 366 high on factor 2, but hard to see which traits should have high scores
(unless we remember).

- Idea 1: use percentile ranks as before.

* Idea 2: Rating scale is easy to interpret. So
*tidy* original data frame to make easier to look
things up.



##  Tidying original data

\scriptsize
```{r bFactor-38}
bem %>%
ungroup() %>% 
mutate(row = row_number()) %>%
pivot_longer(c(-subno, -row), names_to="trait", 
             values_to="score") -> bem_tidy
bem_tidy
```
\normalsize
 


##  Recall data frame of loadings

\footnotesize
```{r bFactor-39}
loadings %>% slice(1:10)
```
\normalsize
 

Want to add the factor scores for each trait to our tidy data frame
`bem_tidy`. This is a left-join (over), matching on the column
`trait` that is in both data frames (thus, the default):


##  Looking up loadings

\scriptsize
```{r bFactor-40}
bem_tidy %>% left_join(loadings) -> bem_tidy
bem_tidy %>% sample_n(12)
```
\normalsize
 


##  Individual 366, high on Factor 2
So now pick out the rows of the tidy data frame that belong to
individual 366 (`row=366`) and for which the `Factor2`
score exceeds 0.4 in absolute value (our "big" from before):

\scriptsize
```{r bFactor-41}
bem_tidy %>% filter(row == 366, abs(Factor2) > 0.4)
```
\normalsize

As expected, high scorer on these.


##  Several individuals
Rows 311 and 214 were *low* on Factor 2, so their scores should
be low. Can we do them all at once?


\scriptsize
```{r bFactor-42}
bem_tidy %>% filter(
row %in% c(366, 311, 214),
abs(Factor2) > 0.4
)
```
\normalsize
 

Can we display each individual in own column?


##  Individual by column
Un-`tidy`, that is, `pivot_wider`:

\tiny
```{r bFactor-43}
bem_tidy %>%
filter(
  row %in% c(366, 311, 214),
  abs(Factor2) > 0.4
) %>%
select(-subno, -Factor1, -Factor2) %>%
pivot_wider(names_from=row, values_from=score)
```
\normalsize

366 high, 311 middling, 214 (sometimes) low.


##  Individuals 230, 258, 359
These were high, low, low on factor 1. Adapt code:

\tiny
```{r bFactor-44}
bem_tidy %>%
filter(row %in% c(359, 258, 230), abs(Factor1) > 0.4) %>%
select(-subno, -Factor1, -Factor2) %>%
pivot_wider(names_from=row, values_from=score)
```
\normalsize



##  Is 2 factors enough?
Suspect not:
```{r bFactor-45 }
bem.2$PVAL
```

 

2 factors resoundingly rejected. Need more. Have to go all the way to
15 factors to not reject:

```{r bFactor-46 }
bem %>%
select(-subno) %>%
factanal(factors = 15) -> bem.15
bem.15$PVAL
```



Even then, only just over 50\% of variability explained.

## What's important in 15 factors?

- Let's take a look at the important things in those 15 factors.

- Get 15-factor loadings into a data frame, as before:  

\small
```{r bFactor-47}
bem.15$loadings %>% 
unclass() %>% 
as_tibble() %>% 
mutate(trait = rownames(bem.15$loadings)) -> loadings
```
\normalsize
 

- then show the highest few loadings on each factor.



##  Factor 1 (of 15)

\footnotesize
```{r bFactor-48}
loadings %>%
arrange(desc(abs(Factor1))) %>%
select(Factor1, trait) %>%
slice(1:10)
```
\normalsize
 
Compassionate, understanding, sympathetic, soothing: thoughtful of
others. 



##  Factor 2

\footnotesize
```{r bFactor-49}
loadings %>%
arrange(desc(abs(Factor2))) %>%
select(Factor2, trait) %>%
slice(1:10)
```
\normalsize
 

Strong personality, forceful, assertive, dominant: getting ahead. 



##  Factor 3

\footnotesize
```{r bFactor-50}
loadings %>%
arrange(desc(abs(Factor3))) %>%
select(Factor3, trait) %>%
slice(1:10)
```
\normalsize
 

Self-reliant, self-sufficient, independent: going it alone.



##  Factor 4

\footnotesize
```{r bFactor-51}
loadings %>%
arrange(desc(abs(Factor4))) %>%
select(Factor4, trait) %>%
slice(1:10)
```
\normalsize
 

Gentle, tender, warm (affectionate): caring for others.



##  Factor 5

\scriptsize
```{r bFactor-52}
loadings %>%
arrange(desc(abs(Factor5))) %>%
select(Factor5, trait) %>%
slice(1:10)
```
\normalsize
 

Ambitious, competitive (with a bit of risk-taking and individualism):
Being the best.



##  Factor 6

\scriptsize
```{r bFactor-53}
loadings %>%
arrange(desc(abs(Factor6))) %>%
select(Factor6, trait) %>%
slice(1:10)
```
\normalsize
 

Acts like a leader, leadership ability (with a bit of Dominant):
Taking charge.



##  Factor 7

\footnotesize
```{r bFactor-54}
loadings %>%
arrange(desc(abs(Factor7))) %>%
select(Factor7, trait) %>%
slice(1:10)
```
\normalsize
 
Happy and cheerful.


##  Factor 8

\footnotesize
```{r bFactor-55}
loadings %>%
arrange(desc(abs(Factor8))) %>%
select(Factor8, trait) %>%
slice(1:10)
```
\normalsize
 
Affectionate, flattering: Making others feel good.



##  Factor 9

\footnotesize
```{r bFactor-56}
loadings %>%
arrange(desc(abs(Factor9))) %>%
select(Factor9, trait) %>%
slice(1:10)
```
\normalsize
 

Taking a stand.



##  Factor 10

\footnotesize
```{r bFactor-57}
loadings %>%
arrange(desc(abs(Factor10))) %>%
select(Factor10, trait) %>%
slice(1:10)
```
\normalsize

 

Feminine. (A little bit of not-masculine!)



##  Factor 11

\footnotesize
```{r bFactor-58}
loadings %>%
arrange(desc(abs(Factor11))) %>%
select(Factor11, trait) %>%
slice(1:10)
```
\normalsize
 

Loyal.



##  Factor 12

\footnotesize
```{r bFactor-59}
loadings %>%
arrange(desc(abs(Factor12))) %>%
select(Factor12, trait) %>%
slice(1:10)
```
\normalsize
 

Childlike. (With a bit of moody, shy, not-self-sufficient, not-conscientious.)



##  Factor 13

\footnotesize
```{r bFactor-60}
loadings %>%
arrange(desc(abs(Factor13))) %>%
select(Factor13, trait) %>%
slice(1:10)
```
\normalsize
 

Truthful. (With a bit of happy and not-gullible.)



##  Factor 14
\footnotesize

```{r bFactor-61}
loadings %>%
arrange(desc(abs(Factor14))) %>%
select(Factor14, trait) %>%
slice(1:10)
```
\normalsize
 

Decisive. (With a bit of self-sufficient and not-soft-spoken.)



##  Factor 15
\footnotesize
```{r bFactor-62}
loadings %>%
arrange(desc(abs(Factor15))) %>%
select(Factor15, trait) %>%
slice(1:10)
```
\normalsize
 

Not-compassionate, athletic, sensitive: A mixed bag. ("Cares about self"?)


##  Anything left out? Uniquenesses

\scriptsize
```{r bFactor-63}
enframe(bem.15$uniquenesses, name="quality", value="uniq") %>%
  slice_max(uniq, n = 10)
```
\normalsize
 

Uses foul language especially, also loves children and analytical. So
could use even more factors.



