Boston Housing

• In an earlier tutorial we considered the Boston Housing data.
• Each observation is a town (suburb) in the Boston metropolitan area.
• There are 14 variables measuring demographic information as well as other factors that may influence house price.

PCA on Boston Housing

#First load required packages
library(tidyverse)
Boston<-readRDS('Boston.rds')
Boston%>%
  column_to_rownames('Town')%>% 
  prcomp(scale.=TRUE)->pcaout
screeplot(pcaout,type = 'l')

Scree Plot

Biplot

biplot(pcaout)

Discussion

• Nearly 60% of the variation of all variables in explained by just 2 PCs.
• Can these PCs be interpreted.
• Sometimes they can but in this example it is difficult.
• This is not surprising, PCA just finds a linear combination that maximises variances.
• To obtain factors with some interpretation we need a more detailed model.

Factor Model

• The factor model is defined as

$$y_{ij}={\lambda }_{j1}{f}_{1i}+{\lambda }_{j2}{f}_{2i}+\dots +{\xi }_{ij}$$

• Or in matrix form

$$y_i=\Lambda f_i+{\xi }_i$$

• $y$ are observed data, $\Lambda$ / $\lambda$ are coefficients, $f$ are latent factors, $\xi$ are error terms.
• The intercept is left out for simplicity.

Notation

• The subscript $i$ denotes the $i^{th}$ cross sectional unit (in the Boston data the town).
• The subscript $j$ denotes the variable (e.g. teacher ratio, distance from downtown etc.)
• The dimensions of $y_i$ and ${\xi }_i$ are $p\times 1$ (or $14\times 1$ in the Boston data).
• If there are $r$ factors then $f_i$ is $r\times 1$ and $\Lambda$ is $p\times r$.
• Verify that all matrix multiplication is conformable.

Regression

• This is similar to a regression model. However
  • In a regression model there are $x$ on the right hand side that are observed.
  • In a factor model these are replaced with $f$ that are unobserved.
• How can we estimate this model?

Assumptions: Errors

• Each idiosyncratic error has its own variance.
  • These variances are called unique variance or uniquenesses
• The idiosyncratic errors are uncorrelated with each other.
  • This is a crucial assumption
• Together these imply that $\text{Var-Cov}\left(\xi \right)=\Psi$ is diagonal.

Assumptions: Factors

• The factor and idiosyncratic errors are uncorrelated.
  • This is similar to regression
• Each factor has a variance of 1.
  • This is harmless since the factor is latent.
• The factors are uncorrelated with each other
  • We relax this assumption later on.
• These imply that $\text{Var-Cov}\left(f\right)=I$.

Estimation

• In general these assumptions imply that

$$E\left(y{y}^{\prime }\right)=\Sigma =\Lambda {\Lambda }^{\prime }+\Psi$$

• The variance is decomposed into two parts
• Part explained by common factors $\Lambda {\Lambda }^{\prime }$ .
  • This is often called the communality or common variance.
• Part unexplained by common factors $\Psi$.
  • This is often called the uniqueness or unique variance.

Estimation

• It is straightforward to estimate $\Sigma$ with its sample equivalent $S$
• We can then choose values $\hat{\Lambda }$ and $\hat{\Psi }$ so that $\hat{\Lambda }{\hat{\Lambda }}^{\prime }+\hat{\Psi }$ is close to $S$.
• There are many ways to do this
• Maximum likelihood estimation is one of the most popular.

Estimation issues

• Using Maximum likelihood estimation does require a distributional assumption about the data.
• The most common assumption is that the data are normally distributed.
• Even when this assumption does not hold the maximum likelihood estimate is still quite robust as long as the data do not differ too much from normality.

Number of factors

• There are a number of strategies for selecting factors
  • Scree plot
  • Kaiser rule
  • Hypothesis tests

Heywood cases

• In some rare cases the maximum likelihood converges to an estimate where the unique variances are negative.
• These are known as Heywood cases
• Since a variance is cannot be negative this is usually caused by
  • Selecting too many factors
  • Too small a sample size.

Using R

• Many packages in R do factor analysis
• We use factanal from the stats package
• First step use the following code

#First load required packages
Boston<-readRDS('Boston.rds')
Boston%>%
  column_to_rownames('Town')%>% 
  factanal(factors = 2,scores = 'none',
           rotation = 'none')->facto

Output

facto$loadings

## 
## Loadings:
##         Factor1 Factor2
## CRIM     0.548         
## ZN      -0.401   0.455 
## INDUS    0.752  -0.513 
## CHAS            -0.120 
## NOX      0.630  -0.628 
## RM      -0.280   0.371 
## AGE      0.531  -0.576 
## DIS     -0.611   0.573 
## RAD      0.923   0.282 
## TX       0.979   0.103 
## PTRATIO  0.472   0.251 
## B       -0.459         
## LSTAT    0.532  -0.384 
## MEDV    -0.622   0.297 
## 
##                Factor1 Factor2
## SS loadings      5.071   2.080
## Proportion Var   0.362   0.149
## Cumulative Var   0.362   0.511

Output

• An advantage of printing the loadings like this is that values close to zero are surpressed.
• This will help with the interpretation of factors.
• For 2 factors, it can be useful to also plot the factors.
• To prepare the data use the tidy function in the broom package.

Loadings

library(broom)
fa_df<-tidy(facto) #Get into data frame

Plotting

The plot is clearer if arrows are used

ggplot(fa_df,aes(x=fl1,y=fl2,
                 label=variable))+
  geom_segment(aes(xend=fl1,
                   yend=fl2,x=0,y=0),
               arrow = arrow())+
  geom_text(color='red',nudge_y = -0.05)

Plotting

Interpretation

It is difficult to interpret these Factors

• Factor 1 seems to take everything into account except the Charles river dummy.
• Factor 2 takes everything into account except crime and race.
• It would be easier to interpret the factors if Factor 1 loaded onto a small set of variables and Factor 2 loaded onto a different small set of variables.
• Can we do this?

Rotations

Recall the model is

$$y_i=\Lambda f_i+{\xi }_i$$

Assume there is an r × r rotation matrix $R$. Since $R^{\prime }R=I$ the model above is equivalent to

$$y_i=\Lambda R^{\prime }R{f}_i+{\xi }_i$$

The rotation trick

Grouping parts together we have

$$y_i=\left(\Lambda R^{\prime }\right)\left(R{f}_i\right)+{\xi }_i$$

Now we have new loadings $\tilde{\Lambda }=\Lambda R^{\prime }$ and new factors $\widetilde{f_i}=R{f}_i$

• All rotated versions of the loadings and factors explain the data equally well and satisfy all assumptions of the model.

Varimax

• Some rotated versions of the factors may be easier to interpret.
• Generally if there are many zero loadings, then the factors are easy to interpret.
• An algorithm known as varimax tries to find a rotation with as many loadings close to zero as possible.
• It can be implemented using the rotation='varimax' option in factanal.

Varimax in R

Boston%>%
  column_to_rownames('Town')%>% 
  factanal(factors = 2,scores = 'none',
           rotation = 'varimax')->facto_vari

Loadings

Varimax

Oblique Rotation

• An orthogonal rotation did not work so instead of considering a matrix where $RR=I$, consider a matrix $G{G}^{-1}=I$ $$y_i=\Lambda G{G}^{-1}{f}_i+{\xi }_i$$
• Now we have new loadings $\tilde{\Lambda }=\Lambda G$ and new factors $\widetilde{f_i}=G^{-1}{f}_i$
• By setting rotation='promax' in factanal, an oblique 'rotation' can be carried out.

Varimax in R

Boston%>%
  column_to_rownames('Town')%>% 
  factanal(factors = 2,scores = 'none',
           rotation = 'promax')->facto_promax

Loadings

Promax

Possible Interpretation

• Factor 1 is
  • Positively correlated with age.
  • Negatively correlated with distance.
• Factor 1 is a geographic factor.
• Factor 2 is
  • Positively correlated with crime, pupil-teacher ratio.
  • Negatively correlated with the race variable.
• Factor 2 is a socioeconomic factor.

More than 2 factors

• If there are more than 2 factors look at the loadings matrix.
• The pattern of zeros should give some clue to the interpretation of the factors.
• Also look for large loadings (in absolute value)

Oblique rotation

• Oblique rotations will lead to factors that are correlated with one another.
• This is not the case for orthogonal factors.
• Other rotation options are available by downloading the package GPArotation
• Orthogonal Rotations: Varimax, Quartimax, Equimax
• Oblique Rotations: Promax, Oblimin, Quartimin, Simplimax

Factor scores

• The factor scores themselves can be estimated using a variety of methods. Two are available as options in the factanal function.
  • Regression Scores
  • Bartlett’s Scores
• Bartlett’s scores are unbiased estimates
• These can be implemented setting scores='regression' or scores='Bartlett' in factanal.

Estimation alternatives

Other estimation methods can also be used for the factor model.

• One example is Principal Axis Factoring, which is available for R using the psych package.
• Principal Axis Factoring does not require the normality assumption and can be adapted for item response data such as Likert scales.

Extended topics

• What we have discussed today is often called exploratory factor analysis.
• In many social sciences the latent variables may themselves influence other observed variables.
• Such models are called structural equation models.
• They can also be estimated by maximum likelihood.