# Script to run principal component analysis (PCA) with FactoMineR.

##### PREPARATIONS ####
####Install packages (only once if you haven't done it before)
#install.packages("FactoMineR")
#install.packages("factoextra")
#install.packages("ggrepel")
#install.packages("data.table")
#install.packages("foreign")

library(FactoMineR)
library(factoextra)
library(tidyverse)
library(foreign)
library(ggrepel)

## For this exercise, we use a dataset from 
## https://stats.oarc.ucla.edu/r/dae/multinomial-logistic-regression/
## The data set contains variables on 200 students. 
## We will build a geometric model of their scores on reading (read),
## writing (write), mathematics (math), science (science) and 
## social studies (socst). Once we have this model, we can use it to predict
## an outcome, such as program selection (prog), or we can use a background 
## variable, such as social class or "SES", to predict the scores.

hsbdemo <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta") 

head(hsbdemo)

pca_hsbdemo <- PCA(
  X = hsbdemo[,c(3, 5:10)], # select scores and...
  quali.sup = c(1:2)) # ses:prog


fviz_eig(
  X = pca_hsbdemo,
  addlabels = TRUE
)

fviz_contrib(
  X = pca_hsbdemo,
  choice = "var",
  axes = c(1)
)

fviz_pca_var(
  X = pca_hsbdemo
)


fviz_ind(
  X = pca_hsbdemo,
  axes = c(1, 2),
  label = "none",
  habillage = "prog",
  addEllipses = TRUE
)

# So, what can we say about social class and selected program?
# Let's look at the eta2.

pca_hsbdemo$quali.sup$eta2