---
title: "Multivariate Data Analysis via Stan Modelling Software"
author: "Sean van der Merwe"
date: "29 November 2021"
output:
  ioslides_presentation:
    self_contained: false
    logo: UFSlogoRunner.png
css: svdmprestemplate.css
---

```{r, include=FALSE}
library(knitr)
library(formatR)
opts_chunk$set(tidy.opts=list(width.cutoff=55),tidy=TRUE)
# l10n_info()
# rmarkdown::render('MDAGstan2021.Rmd')
```

## Enjoying this presentation {.smaller}

[seanvdm.co.za/presentation](https://seanvdm.co.za/presentation)

Please open the presentation on your own computer (type link above) to enjoy the following benefits:

-   Exploring the data yourself

-   Interacting with the graphs yourself

-   Moving at your own pace (press *o* for overview)

-   Maximum quality

    -   Try *f* or *F11* for full screen
    -   Press *w* at any time to change to wide screen mode and back
    -   *Ctrl + Mouse Wheel* to zoom

The presentation and source code will be available for download on my website at [seanvdm.co.za/files/workshops/mdag2021](https://seanvdm.co.za/files/workshops/mdag2021)

# Introduction

## Abstract

Latent variable models are useful when you have excess observed variables whose values are well explained by a smaller number of underlying, but unobserved, variables. The goal is to simplify the multivariate structure, often to the point where it can be more easily visualised. These models have some similarities to Factor Analysis and Principal Components Analysis, but are parametric models. This presentation discusses an implementation of these models, along with various expansions to handle issues typical of observed data. Expansions discussed include Likert style responses, robust fits, missing value imputation, and possible future directions. 

A lot of introductory and supplementary topics are also discussed.

## Illustration

```{r, echo=FALSE, message=FALSE}
library(visNetwork)
cols <- viridisLite::viridis(7)[3:6]
nodes <- data.frame(id = 1:9, 
                  label = c('Data Generating Process', 'Latent Variable 1', 'Latent Variable 2', 'Random Noise', paste('Observed Var', 1:5)),
                    group = rep(c('Hidden', 'Observed'), c(4,5)),
                    shape = rep(c("box", "ellipse"), c(4,5)),
                    color = cols[c(1,2,2,3,4,4,4,4,4)],
                  title = c('🏏🏃⚾🤾', '🏏', '⚾', '☔🌈' ,'Score 1', 'Score 2', 'Score 3', 'Score 4', 'Score 5')
                    )
edges <- data.frame(from = c(1,1,rep(2:4, e=5)),
                    to   = c(2,3,rep(5:9, t=3)),
                    arrows = "to", 
                    length = c(20, 20, rep(100,10), rep(200, 5))
                    )
net1 <- visNetwork(nodes, edges)
net1 <- visGroups(net1, groupname = 'Hidden', color = cols[3], shape = 'box')
net1 <- visGroups(net1, groupname = 'Observed', color = cols[4], shape = 'ellipse')
net1 <- visLegend(net1)
library(widgetframe)
frameWidget(net1)
```


## Despair

-   Have you faced non-Gaussian data structures (zero inflation / censoring / extremes *etc.*) with weird dependence relationships (non-linear mixed effects say)?
-   Maybe data on a scale, such as a Likert scale?
-   Do you maybe have multiple response variables? Perhaps too many and a need for dimension reduction?
-   Have you ever looked at the spread of missing values and wanted to cry in despair?
-   **What would you do if you encountered ALL of these AT THE SAME TIME?**

No need for despair, I'll show you that you can actually analyse such data in a single unified joint model and answer all the desired statistical questions in a sensible way.

## Why Stan?

If all the nuts you ever have to tighten are the same one or two sizes then you buy the best spanners in those sizes. You don't use a monkey wrench that needs to be set each time. But if you're encountering different sizes of nuts, fittings, and other objects that need tightening then you get a monkey wrench. In statistical modelling, Stan is the monkey wrench.

According to <https://mc-Stan.org/>: "Stan is a state-of-the-art platform for statistical modelling and high-performance statistical computation."

You type out your model in Stan language and then Stan does all the work, including: simulating the whole posterior distribution, simulating the approximate posterior, and/or Maximum Likelihood.

## Advantages of Stan

-   It compiles your model into C++ code, so it is very fast
-   Latest simulation methods (HMC+NUTS)
-   Can do usual Bayes, Variational Bayes or Maximum Likelihood as you like - You do not change your model, just one line of code
-   Most importantly, **if you want to change your model then you change only the part of the model that needs changing**. No changing software, no installing new packages, no changing notation for a different 'PROC', etc.

## Where Stan might help {.smaller}

Suppose you are doing a regression, but

-   Perhaps a Normal distribution is not appropriate; maybe a t, gamma, GPD, *etc.* model makes more sense
-   Maybe your observations belong to groups with differing unknown variances
-   Maybe you need to reduce the dimension or work with latent variables
-   Perhaps you have some prior information on one or more parameters
-   Perhaps you have change points, or censoring, or random effects

While tools exist for dealing with each of these issues, each issue has its own tools with their own notations and limitations. Few of these tools can handle any of the other issues at the same time.

With Stan you can handle all these issues **at the same time**.

## Getting help {.larger}

The Stan language is intuitive and sensible, but you will still need to look up most things the first time you want to use them.

The key bookmark is <https://mc-stan.org/users/documentation/> which has everything you need.

I find myself visiting the language functions reference often to see the exact notation for the distribution I want to use.

For the other issues, a standard internet search gets you very far. Stan is popular free software so it has a large and active user base ready to help.

There are many helper packages that can do a lot of (or all) the work for you, including creating the model code for you to play with.

# Running Stan

## Interfaces

Stan is currently available for R, Python, and C++ as far as I can tell, with other platforms being able to call it via one of those three. We will focus on *rstan* as it's the easiest for most statisticians.

There are many ways to use Stan in R:

-   You can type your Stan model in a separate file, then tell rstan to compile that file
-   You can type your Stan model in a Stan block in your R Markdown document and click the play button to compile it
-   You can save the compiled model in a .Rds or Rdata file and then load the compiled model each time you need it
-   You can mix and match the above as it suits you

## R Markdown interface

I will use R Markdown to demonstrate.

I use R Markdown for almost everything, including teaching, assessment, consultation, some research, websites, and even making presentations like this one.

To use R Markdown you open RStudio, click File -\> New File -\> R Markdown...

To put in a Stan model you click on the Insert button at the top right of the code window, then click *Stan*.

On the next slide I'll show a very basic Stan model in RStudio.

```{=html}
<!--
, eval=FALSE, include=FALSE
-->
```
```{stan, output.var='studentt', eval=FALSE, include=FALSE}
//
// This Stan program defines a t fit to a vector of observations.
// The code below is by Sean van der Merwe, UFS
//
// The input data is a vector 'x' and its length 'n'.
data {
  int<lower=0> n;   			    // number of observations
  vector[n] x;      				  // observations
}
// The parameters accepted by the model.
parameters {
  real<lower=0> v;            // degrees of freedom
  real<lower=0> s;            // scale parameter
  real mu;                    // location parameter
}
// The model to be estimated.
model {
  x ~ student_t(v,mu,s);      // likelihood
  v ~ exponential(0.01);      // prior on degrees of freedom
}
```

## First model: t density

Notice the play button in the top right corner (circled in red).

```{r drawing1, out.width = "95%", echo=FALSE}
include_graphics("stanblock.png")
```

## Model breakdown

The critical components are:

1.  The data
2.  The parameters
3.  The likelihood

The optional components are:

-   Priors for some or all parameters - default is uniform on the domain of the parameter
-   Transformed Parameters - for when your model gets complicated
-   Transformed Data - better to do in R in advance
-   Generated Quantities - for creating things you want later

## Installation and Loading

To be able to use *rstan* the way I do:

1.  Install R
2.  Install Rtools (only needed on Windows, get it where you get R)
3.  Install RStudio
4.  In the RStudio Console window, type *install.packages('rstan')*
5.  In your code, type *library(rstan)*

## Number of cores

-   Optionally, set the number of cores for Stan to use when simulating in parallel
-   Easiest is to pick a number: *options(mc.cores = 3)*
-   On most computers you can get a big speed increase easily if you match the number of cores and chains
-   So calculate and store a good number: library(parallel); mycores \<- ceiling(0.75\*detectCores(logical = FALSE))
-   Then use that number everywhere: *options(mc.cores = mycores)* and later on when simulating set *chains = mycores*
-   I like 75% because then I can still use my computer while simulating (*i.e.* watch YouTube while I'm technically working).

## Log returns

Let's apply the t model to some data.

```{r tapplycode, eval=FALSE}
library(openxlsx)
### <b>
sourcedata <- read.xlsx('IBM19891998.xlsx')
lret <- sourcedata$LogReturns
### </b>
dates <- sourcedata$year + sourcedata$month/12 + sourcedata$day/365

cols <- viridisLite::viridis(3)
par(mfrow=c(2,1), mar=c(4,4,2,0.2))
plot(dates, lret, type='l', col=cols[1], lwd=1, xlab='Date', ylab='Log return', main='IBM')
hist(lret, 50, col=cols[2], density = 20, ylab='Frequency', xlab='Log return', main='IBM')
```

## IBM data

```{r tapplyrun, echo=FALSE}
library(openxlsx)
sourcedata <- read.xlsx('IBM19891998.xlsx')
lret <- sourcedata$LogReturns
dates <- sourcedata$year + sourcedata$month/12 + sourcedata$day/365

cols <- viridisLite::viridis(3)
par(mfrow=c(2,1), mar=c(4,4,2,0.2))
plot(dates, lret, type='l', col=cols[1], lwd=1, xlab='Date', ylab='Log return', main='IBM')
hist(lret, 50, col=cols[2], density = 20, ylab='Frequency', xlab='Log return', main='IBM')
```

## Fitting the model

```{r loadstanforreals, include=FALSE}
library(parallel); mycores <- max(1,floor(0.75*detectCores(logical = FALSE))) # Choose the number of computer cores to use, 75% of capacity is a good rule of thumb
suppressPackageStartupMessages(library(rstan,warn.conflicts=FALSE,quietly=TRUE))
options(mc.cores = mycores) 
rstan_options(auto_write = TRUE)
```

```{r runstan1, eval=FALSE}
stan_data <- list(x=lret, n=length(lret)) # Create the data list

# Maximum Likelihood / Posterior Mode estimation:
fitML1 <- optimizing(studentt, stan_data)

# This is the main simulation function:
fitBayes1 <- sampling(studentt, stan_data, chains = mycores) 
# Type ?sampling for a lot of additional options

saveRDS(fitBayes1,'studenttIBM.Rds') # Save model fit

summary(fitBayes1)$summary # See results of simulation

traceplot(fitBayes1) # See simulation paths to check for stability

list_of_draws1 <- extract(fitBayes1) # Extract simulations
```

```{r loadstanstudentt, include=FALSE}
fitBayes1 <- readRDS('studenttIBM.Rds')
list_of_draws1 <- extract(fitBayes1)
```

## t model summary

```{r showfittable, echo=FALSE}
s1 <- summary(fitBayes1)$summary
kable(round(s1[1:3,c(1,2,3,4,8,10)],3))
```

The *Rhat* measure should always be close to 1.

*Rhat* compares the variation between chains to the variation within chains, as well as the variation in different parts of each chain. If these are all the same then that is a good indication of convergence **in distribution**.

## t model trace plot {.smaller}

```{r showfittraceplot, echo=FALSE, fig.width=7, fig.height=2.5}
traceplot(fitBayes1, 'mu')
```

On the trace plot you should see the same patterns on the left as on the right to indicate that you have done enough burn-in.

You should also see the same sort of patterns in each chain.

When the chains are viewed together they should cover the whole posterior distribution.

## Non-standard distributions

We don't need to only use standard distributions in Stan.

For example, we could use a better prior on the degrees of freedom (provided by Prof Abrie van der Merwe):

$$p(\nu)\propto \left[\psi'(0.5\nu)-\psi'(0.5(\nu+1))-2(\nu+3)\nu^{-1}(\nu+1)^{-2}\right]^{0.5}$$

where $\psi'$ is the trigamma function.

```{stan, output.var='nonsense1', include=FALSE, eval=FALSE}
model {
  x ~ student_t(v,mu,s);      // likelihood
  target += 0.5*log(trigamma(0.5*nu) - trigamma((nu+1)*0.5) - 2*(nu+3)/nu/(nu+1)/(nu+1));
  // log prior is added to the log likelihood to create log posterior (target)
}
```

```{r weirdpriorinstan, out.width = "100%", echo=FALSE}
include_graphics("stanblockprior.png")
```

## Working with the results {.smaller}

```{r plotresults1, fig.width=7, fig.height=3.5}
# Posterior Predictive Distribution
postpred1 <- rt(length(list_of_draws1$mu),list_of_draws1$v)*list_of_draws1$s + list_of_draws1$mu
par(mar=c(4.4,4.4,0.4,0.4))
plot(density(postpred1),xlab='Future values',ylab='Posterior density',main='',col='green',lwd=3)
```

## Intervals and p-values

After fitting the model, we can obtain intervals for any quantity of interest as follows:

-   Calculate the quantity for each simulated parameter vector, to get a vector of quantities
-   This vector of quantities can be seen as an approximation to the posterior distribution of the quantity of interest, let us call it *sims*.
-   Calculate a symmetric interval, say 95%, using the built-in quantile function: *quantile(sims,c(0.025,0.975))*
-   Calculate an equivalent to the traditional two-sided p-value using this function:

```{r pvalfunc}
pvalfunc <- function(sims, target=0) { 2*min( mean(sims<target), mean(sims>target) ) }
```

## Shortest interval {.smaller}

For a given coverage, a shorter interval provides more useful information.

If the distribution of the quantity of interest is skew then the symmetric interval is not ideal.

We can calculate the shortest interval, also known as the highest posterior density (HPD) interval, quite easily though if we have simulations (also works for bootstraps).

The process is simple: to get the shortest 95% interval, just work out the length of every 95% interval and see which one is the shortest!

```{r intervalfunc}
shortestinterval <- function(sims, alpha=0.05) { # Coded by Sean van der Merwe, UFS
sorted.sims <- sort(sims)
nsims <- length(sims)
gap <- round(nsims*(1-alpha))
widths <- diff(sorted.sims,gap)
interval <- sorted.sims[c(which.min(widths),(which.min(widths) + gap))]
return(interval) }
```

## Simulate once, calculate many

Once you have a set of simulations, you can calculate estimates, intervals, and probabilities to your heart's content:

1.  Assume the first set of simulated parameters are the true values.
2.  Calculate the quantity of interest using those parameters.
3.  Repeat for simulations 2 to ... and keep all the answers.
4.  Summarise the vector of answers using *mean*, *summary*, *quantile*, *hist*, *etc.*

Example: to calculate $P[\mu^2>5]$ we type

```{r}
mean( list_of_draws1$mu^2 > 5 )
```

```{r, echo=FALSE}
rm(fitBayes1, list_of_draws1, s1, sourcedata, cols, dates, lret, postpred1, mycores)
```

# Binomial responses

## Likert Scale as Binomial

Assume we are measuring something on a scale and we have reason to believe that the underlying process is continuous but forced into the scale at some stage in the data generating process.

Suppose we have Likert scale questions with 5 levels, then we can mark them as 0, 1, 2, 3, and 4.

In the Binomial GLM formulation it is assumed that the responses $y_i\sim Binomial(4, p_i)$, while $logit(p_i)=f(\mathbf{x}_i)$ where $f(\mathbf{x}_i)$ is some function of the explanatory variables corresponding to observation $i$ (linear function for a true GLM). In the simplest case we have a constant/intercept only: $f(\mathbf{x}_i)=\beta_0$, but we could just as easily include predictors such as age or gender in this framework.

## Generate raw data {.smaller}

Suppose we ask some students whether they have a good lecturer.

```{r}
set.seed(12345)
n <- 100
Age <- rpois(n,3) + 18
Gender <- round(runif(n))
Likert <- rbinom(n, 4, plogis( -2 + 0.1*Age + 0.3*Gender))
```

```{r, include=FALSE}
UFScolours <- c(Black="#000000FF", Blue="#0F204BFF", Red="#A71930FF", LightGray="#A7A8AAFF", DarkGray="#8D8D8EFF", NAS="#0039A7FF", EDU="#00675AFF", THEO="#C69317FF", HUM="#EA8400FF", SA="#A40084FF", SC="#9E83B7FF", QWA="#00B140FF", LAW="#BB133EFF", HLTH="#490E6FFF", VeryLight="#D7D8DDFF")


plotly_table_UFS <- function(datafrm, colors=UFScolours, useRowNames = FALSE, rowNamesHeading = " ") {
  # For the colour palette: ensure that the first colour is dark and the last colour is light
  datafrm <- as.data.frame(datafrm)
  nms <- names(datafrm)
  if (useRowNames) {
    datafrm <- data.frame(row.names(datafrm), datafrm, row.names = NULL)
    nms <- c(rowNamesHeading, nms)
    names(datafrm) <- nms
  }
  k <- ncol(datafrm)
  algn <- rep('center', k)
  cell_values <- rbind(t(as.matrix(unname(datafrm))))
  fig <- plot_ly(
  type = 'table',
  header = list(
    values = nms,
    align = algn,
    line = list(width = 1, color = colors[3]),
    fill = list(color = colors[2]),
    font = list(family = "Arial", size = 14, color = "#FFFFFFFF")
  ),
  cells = list(
    values = cell_values,
    align = algn,
    line = list(color = colors[3], width = 1),
    fill = list(color = list(rep(c("#FFFFFFFF", colors[length(colors)]), ceiling(length(cell_values)/2)))),
    font = list(family = "Arial", size = 12, color = colors[1])
  ))
fig
}
```

::: {.centered}
::: {.columns-2}
```{r, echo=FALSE, message=FALSE, fig.width=4, fig.height=3}
library(plotly)
plotdata <- as.data.frame(table(Likert,Gender))
fig <- plot_ly(plotdata, x = ~Gender, y = ~Freq, type = 'bar', color = ~Likert, colors=unname(UFScolours[1:5]) )
fig <- layout(fig, yaxis = list(title = 'Frequency'), barmode = 'stack')
fig
```

```{r, echo=FALSE, message=FALSE, fig.width=4, fig.height=3}
fulldata <- data.frame(Likert=Likert, Age=Age, Gender=Gender, row.names = NULL)
tabulateddata <- aggregate(fulldata, by=as.list(fulldata), FUN='length')
tabulateddata <- tabulateddata[1:4]
names(tabulateddata)[4] <- 'Count'
plotly_table_UFS(tabulateddata)
```
:::
:::

## Response by age

```{r, echo=FALSE}
plotdata <- data.frame(Likert=Likert, Age=factor(Age), Gender=Gender, row.names = NULL)
cols <- viridisLite::viridis(9)
fig <- plot_ly(plotdata, x = ~Age, y = ~Likert, type = 'box', color = ~Age, colors = cols)
fig
```

## All at once

```{r, echo=FALSE, fig.width=8, fig.height=6}
fig <- plot_ly(tabulateddata, type = 'scatter3d', mode = 'markers', x = ~Age, y = ~Likert, z = ~Count, color = ~Gender, colors = UFScolours[2:3],
             text = ~paste('Age:', Age, '<br>Response:', Likert, '<br>Gender:', Gender, '<br>Count:', Count))
layout(fig, scene=list(camera=list(center=list(x=0, y=0, z=-0.3))))
```

## Standard GLM

Here we fit a Binomial GLM to the response in the standard way:

```{r}
y <- Likert/4
s1 <- coef(summary(glm(y~Age+Gender, family = binomial(), weights = rep(4,n))))
```

```{r, echo=FALSE}
plotly_table_UFS(round(s1,3), useRowNames = TRUE)
```

## Stan GLM

There are several easy ways to do a GLM using Stan. The *rstanarm* and *brms* packages are highly recommended as easy interfaces for using Stan and for doing Bayes in general.

However, we are going to do it the hardcore way, for reasons that will make sense later. We will code the model in Stan ourselves.

We will use the [the STAN interface](https://mc-stan.org/) via [rstan](https://mc-stan.org/users/interfaces/rstan) to implement our models.

```{r loadSTAN, warning=FALSE, message=FALSE, class.source = "fold-hide"}
library(parallel)
library(rstan)
mycores <- max(1,floor(detectCores(logical = FALSE)*0.75))
options(mc.cores = mycores)
rstan_options(auto_write = TRUE)
```

## Model definition

```{stan output.var='Binomial1', eval=FALSE}
data {
  int<lower=2> n;            // number of observations
  int<lower=1> m;            // maximum value of observations
  int<lower=0, upper=m> y[n]; // observations
  int<lower=1> k;            // number of explanatory variables, including intercept
  row_vector[k] x[n];            // explanatory matrix
}
parameters {
  vector[k] beta;            // coefficients
}
model {
  for (i in 1:n) {
    y[i] ~ binomial_logit(m, x[i]*beta);
  }
}
```

## Simulation results

```{r, eval=FALSE}
BinResults <- sampling(Binomial1, list(n=n, m=4, y=Likert, k=3, x=model.matrix(~Age+Gender)), chains=mycores)
# plot(BinResults)
```

```{r, eval=FALSE, include=FALSE}
saveRDS(BinResults, 'BinResults.Rds')
```

```{r, include=FALSE}
BinResults <- readRDS('BinResults.Rds')
```

```{r}
BinResSims <- extract(BinResults)
```

```{r, echo=FALSE}
BinResTable <- round(summary(BinResults)$summary[1:3, c(1,2,3,4,6,8,9,10)],3)
BinResDF <- data.frame(c('Intercept', 'Age', 'Gender'), BinResTable) 
names(BinResDF) <- c('Coefficient', colnames(BinResTable))
plotly_table_UFS(BinResDF)
```

## Simulation results visual

```{r, echo=FALSE}
plotdata <- data.frame(Age=BinResSims$beta[,2], Gender=BinResSims$beta[,3], row.names = NULL)
fig <- plot_ly(plotdata, y = ~Age, type = 'box', colors = unname(UFScolours[2:3]), name = "Age Coefficient")
fig <- layout(fig, yaxis = list(title = 'Coefficient Simulation'))
add_trace(fig, y = ~Gender, type = 'box', name = "Gender Coefficient", color=UFScolours[3])
```

```{r, echo=FALSE}
rm(BinResDF, BinResSims, BinResTable, BinResults, fig, fulldata, plotdata, s1, tabulateddata, Age, cols, Gender, Likert, n, y)
```

## Pause, ponder, and question {.larger}

::: {.green}

-   How could your work be improved with access to a modelling tool with such flexibility?

:::

::: {.red}

-   What would be the obstacles to using this?

:::

# PCA

## Generate raw data {.smaller}

In this example we are first going to take some multivariate normal data and apply PCA to it.

```{r, message=FALSE}
library(MASS)
set.seed(12345)
under_Dim <- 2
full_Dim <- 5
true_weights <- matrix(runif(full_Dim*under_Dim,-3,3), under_Dim, full_Dim)
n <- 100
rawZ <- mvrnorm(n, rep(0,under_Dim), diag(rep(1,under_Dim)))
rawX <- matrix(rnorm(n*full_Dim), n, full_Dim) + rawZ%*%true_weights
GroupNum <- (rowSums(rawZ)>0)+1
Group <- c('🐬', '🐪')[GroupNum]
```

We are using `r under_Dim` underlying dimensions and random weights to relate those dimensions to the `r full_Dim` observed dimensions.

```{r, echo=FALSE, message=FALSE, warning=FALSE}
plotly_table_UFS(round(true_weights,3))
```

For fun, we create a grouping variable based on whether the underlying variables are positive on average and see whether we can see the separation later.

## Visualising the raw data {.smaller}

Our raw data has 6 dimensions, which makes it a challenge to visualise. We can try though. In this case I used a logistic transform of the 5th variable for size, and you have to mouse over to see the group. It is not easy to see the multivariate structure.

```{r, echo=FALSE, fig.width=8, fig.height=6, warning=FALSE}
rawplotdata <- rawX
colnames(rawplotdata) <- paste0('V', 1:full_Dim)
rawplotdata <- data.frame(rawplotdata, Group=Group, Size=plogis(rawplotdata[,5])+2)
fig <- plot_ly(rawplotdata, type = 'scatter3d', mode = 'markers', x = ~V1, y = ~V2, z = ~V3, color = ~V4, size=~Size, colors = UFScolours[c(12,10)], text = ~Group)
layout(fig, scene=list(camera=list(center=list(x=0, y=0, z=-0.3))))
```

## Scale data {.smaller}

It is often a good idea to scale the data, although not always. We do so here, but not in all cases going forward.

```{r}
X <- scale(rawX)
round(colMeans(rawX), 4)
round(colMeans(X), 4)
```

The means are now zero and the standard deviations one.

::: {.yellow}

-   Random side note: if you are running my code on your computer then replace the emojis with the usual 'Group A' and 'Group B' before knitting.

:::


## Standard PCA

```{r SvdMbiplotPlotly, echo=FALSE}
biplot.plotly <- function(dataset, grouping_points=rep(1,nrow(dataset)), pointtext=rownames(dataset), colpallette=c('virShuffle', 'rainbow', 'vir')[1], expansion=1.16, export_filename='plotlyexport', export_width=600, export_height=400, export_format='svg') {
  # Custom biplot function by Sean van der Merwe, University of the Free State
  # Takes already scaled matrix as input ideally
  PCA <- prcomp(dataset) # Do PCA
  lam <- PCA$sdev[1:2]*sqrt(nrow(dataset)) # Get scale parameter
  PoV <- (PCA$sdev^2)/sum(PCA$sdev^2)
  scores <- t(t(PCA$x[, 1:2])/lam)
  vars <- t(t(PCA$rotation[, 1:2])*lam)
  groupfac <- factor(grouping_points, levels = unique(grouping_points))
  colsvec <- as.numeric(groupfac)
  rangePoints <- apply(scores,2, function(x){max(abs(c(min(x),max(x))))})
  rangeArrows <- apply(vars,2, function(x){max(abs(c(min(x),max(x))))})
  vars <- t(t(vars)*rangePoints/rangeArrows)
  nArrows <- nrow(vars)
  nGroups <- nlevels(groupfac)
  mix_nums <- function(n) {  mix <- rev(order(abs((1:n)-(n+1)/2))); c(mix[n],mix[-n]) }
  colsarrows <- viridisLite::plasma(nArrows); colsarrows <- colsarrows[mix_nums(nArrows)]
  if (length(colpallette) == 1) {
    Finalcolpallette <- switch(EXPR = colpallette, 
                          rainbow = {rainbow(nGroups)},
                          virShuffle = {pal <- viridisLite::viridis(nGroups+2)[c(-1,-(nGroups+2))]; pal[mix_nums(nGroups)]},
                          vir = {viridisLite::viridis(nGroups+2)[c(-1,-(nGroups+2))]}
                        )
  } else {
    Finalcolpallette <- unname(colpallette)
  }
  if (is.null(pointtext)) {pointtext=rep('', nrow(dataset))}
  plotdata <- data.frame(scores, Group=groupfac, Hover=pointtext, row.names = NULL)
  output <- plot_ly(plotdata, x=~PC1, y=~PC2, color = ~Group, text=pointtext, type='scatter', mode='markers', colors = Finalcolpallette)
  output <- add_annotations(output, text=rownames(vars), ax=vars[,1]*expansion, ay=vars[,2]*expansion, x=vars[,1], y=vars[,2], xref='x', yref='y', axref='x', ayref='y', arrowcolor=colsarrows)
  output <- add_annotations(output, text='', x=vars[,1], y=vars[,2], ax=rep(0,nArrows), ay=rep(0,nArrows), xref='x', yref='y', axref='x', ayref='y', arrowcolor=colsarrows)
  output <- layout(output, xaxis=list(title=paste0('PC1 (',round(PoV[1]*100,1),'%)')), yaxis = list(title=paste0('PC2 (',round(PoV[2]*100,1),'%)')))
  output <- config(output, toImageButtonOptions = list(format=export_format, filename=export_filename, width=export_width, height=export_height))
  output
}
```

```{r, echo=FALSE, results='asis'}
colnames(X) <- varNames <- paste('Var', 1:full_Dim)
biplot.plotly(X, grouping_points = Group, pointtext = Group)
```

## Bayesian Latent Variable Model {.smaller}

```{stan, output.var='PCAbasic1', eval=FALSE}
// This Stan block defines a basic Bayesian PCA, Joseph H. Sakaya & Suleiman A. Khan
// https://www.cs.helsinki.fi/u/sakaya/tutorial/code/pca.R 
data {
	int<lower=0> N;    // Number of samples
	int<lower=0> D;    // The original dimension
	int<lower=0> K;    // The maximum latent dimension, must be < D
	matrix[N, D] X;   }// The data matrix
parameters {
	matrix[N, K] Z;    // The latent matrix
	matrix[D, K] W;    // The weight matrix
	real<lower=0> tau; // Noise term 
	vector<lower=0>[K] alpha; } // ARD prior
transformed parameters{
	vector<lower=0>[K] t_alpha;
	real<lower=0> t_tau;
  t_alpha = inv(sqrt(alpha));
  t_tau = inv(sqrt(tau)); }
model {
	tau ~ gamma(1,1);			
	to_vector(Z) ~ normal(0,1);
	alpha ~ gamma(1e-3,1e-3);				
	for(k in 1:K) W[,k] ~ normal(0, t_alpha[k]);
	to_vector(X) ~ normal(to_vector(Z*W'), t_tau); }
```

. . . . . . . . . . . . . . . . . See <https://www.cs.helsinki.fi/u/sakaya/tutorial/> for more info.

## Simulation and processing {.smaller}

```{r, warning=FALSE}
reduced_Dim <- full_Dim - 1
stan_data <- list(N=nrow(X), D=ncol(X), K=reduced_Dim, X=X)
```

```{r, warning=FALSE, eval=FALSE}
modelfit <- sampling(PCAbasic1, stan_data, pars=c('alpha', 'tau', 'W', 'Z'), iter = 4000, chains = mycores)
# modelfit <- vb(PCAbasic1, stan_data, pars=c('alpha', 'tau', 'W', 'Z')) # , algorithm = 'fullrank')
```

```{r, eval=FALSE, include=FALSE}
saveRDS(modelfit, 'PCAexperiment1.Rds')
```

```{r, include=FALSE}
modelfit <- readRDS('PCAexperiment1.Rds')
```

```{r, fig.width=4, fig.height=2.5}
traceplot(modelfit, c('alpha'))
sims <- extract(modelfit)
```

## The non-identifyability problem {.smaller}

-   In a regular PCA the latent dimensions are not naturally in the ideal order, they are sorted after calculation according to the eigenvalues. Similarly, we must sort the latent dimensions here.

-   However, as can be seen on the trace plot, the latent dimension order is constantly changing across simulations.

-   The full posterior distribution includes all possible orderings of the latent dimensions. 😲

    -   People use variational inference to fit the model instead of simulations, often arriving at only a single ordering; but it doesn't fully solve the problem.
    -   Variational inference can be much faster in large samples though, where it also the most accurate. Stan does variational inference in a very clever way.

-   The latent dimensions can also be flipped arbitrarily, and are in the simulations.

## Dimension order {.smaller}

Here we sort each simulation vector by the order of the weight of the dimensions.

```{r}
nsims <- length(sims$tau)
orders <- apply(sims$alpha, 1, order)
for (sim in seq_len(nsims)) {
  sims$W[sim,,] <- sims$W[sim,,orders[,sim]]
  sims$Z[sim,,] <- sims$Z[sim,,orders[,sim]]
  sims$alpha[sim,] <- sims$alpha[sim,orders[,sim]]
}
```

At this stage the weights per dimension are almost perfectly bimodal. Consider the weight of the first variable in the first dimension:

```{r, fig.height=3, fig.width=5, echo=FALSE}
plot_ly(x=sims$W[,1,1], type='histogram', marker = list(color = UFScolours[2], line = list(color = UFScolours[3], width = 1)))
```

## Dimension orientation

So we pick a direction for each dimension and flip dimensions as needed per simulation.

In this case we use the direction of the first variable as the anchor, but this can be changed, or done multiple times for different anchors to be sure.

```{r}
anchor <- 1
for (dim in seq_len(reduced_Dim)) {
  flip <- sims$W[,anchor,dim] < 0
  sims$W[flip,,dim] <- -sims$W[flip,,dim]
  sims$Z[flip,,dim] <- -sims$Z[flip,,dim]
}
```

## Proportions of variation explained {.smaller}

Once the dimensions are sorted we can calculate the proportion of variation explained by each dimension.

```{r}
MedianVariation <- apply(1/(sims$alpha), 2, median)
PoV_Bayes <- MedianVariation/(sum(MedianVariation) + median(1/sims$tau))
```

The cumulative proportion of variation explained by the model is theoretically `r round(PoV_Bayes,3)*100`.

```{r, echo=FALSE}
PoVfunc <- function(obs, pred1, pred12) {
  varhat <- sd(c(obs))^2
  varhat1 <- sd(c(pred1) - c(obs))^2
  varhat12 <- sd(c(pred12) - c(obs))^2
  PoV <- c(varhat-varhat1, varhat1-varhat12)/varhat
  if (!is.null(ncol(obs))) {
    varhat <- apply(obs, 2, function(clm) {sd(clm)^2})
    varhat1 <- apply(pred1 - obs, 2, function(clm) {sd(clm)^2})
    varhat12 <- apply(pred12 - obs, 2, function(clm) {sd(clm)^2})
    PoV <- c(mean((varhat-varhat1)/varhat), mean((varhat1-varhat12)/varhat))
  }
  PoV
}

Zhat <- apply(sims$Z, 2:3, median)
What <- apply(sims$W, 2:3, median)
PoV_Raw <- PoVfunc(X, Zhat[,1] %*% t(What[,1]), Zhat[,1:2] %*% t(What[,1:2]))
```

The proportion of variance explained by the model, calculated for each variable and averaged, is practically `r round(PoV_Raw,3)*100`.

The reason why we resort to this alternative calculation will become clearer soon.

```{r BayesBiplotPlotly, echo=FALSE}
BayesBiplot.plotly <- function(data_sims, var_sims, grouping_points=rep(1,dim(data_sims)[2]), vartext=NULL, pointtext=NULL, PoV = rep("?",2), colpallette=c('rainbow', 'virShuffle', 'vir')[2], expansion=1.16, export_filename='plotlyexport', export_width=600, export_height=400, export_format='svg', rotate=0.0, rectangles=c('none','standard','full')[1]) {
  # Custom biplot function by Sean van der Merwe, University of the Free State
  # Takes simulated scores and simulated factors
  data_sims <- data_sims[,,1:2]
  var_sims <- var_sims[,,1:2]
  data_ests <- apply(data_sims,2:3,mean)
  var_ests <- apply(var_sims,2:3,mean)
  vmax <- matrix(c(1,0,0,1),2)
  if (is.numeric(rotate)) { 
    vmax <- matrix(c(cos(rotate), -sin(rotate), sin(rotate), cos(rotate)),2) 
  } else {
    if (rotate == 'varimax') { vmax <- varimax(var_ests)$rotmat }
    if (rotate == 'promax') { vmax <- promax(var_ests)$rotmat }
  }
  data_ests <- data_ests %*% vmax
  nsims <- length(data_sims[,1,1])
  for (i in nsims) {
    data_sims[i,,] <- data_sims[i,,] %*% vmax
  }
  var_ests <- var_ests %*% vmax
  data_ints <- apply(data_sims, 2:3, shortestinterval, alpha=0.5)
  if (length(PoV) > 20) {
    PoV <- round(PoVfunc(PoV, data_ests[,1] %*% t(var_ests[,1]), data_ests %*% t(var_ests)),3)*100
  }
  maxes <- apply(data_ests, 2, function(d) { max(abs(d)) })
  for (d in 1:2) { var_ests[,d] <- var_ests[,d]*maxes[d]/max(abs(var_ests[,d]))/1.1 }
  groupfac <- factor(grouping_points, levels = unique(grouping_points))
  groupnum <- as.numeric(groupfac)
  nArrows <- nrow(var_ests)
  nGroups <- nlevels(groupfac)
  mix_nums <- function(n) {  mix <- rev(order(abs((1:n)-(n+1)/2))); c(mix[n],mix[-n]) }
  colsarrows <- viridisLite::plasma(nArrows); colsarrows <- colsarrows[mix_nums(nArrows)]
  colpallette <- switch(EXPR = colpallette, 
                        rainbow = {rainbow(nGroups)},
                          virShuffle = {pal <- viridisLite::viridis(nGroups+2)[c(-1,-(nGroups+2))]; pal[mix_nums(nGroups)]},
                          vir = {viridisLite::viridis(nGroups+2)[c(-1,-(nGroups+2))]}
                        )
  if (is.null(pointtext)) {pointtext=rep('', nrow(data_ests))}
  if (is.null(vartext)) {vartext=paste0('Var', seq_len(nArrows))}
  rectlist <- lapply(seq_len(nrow(data_ests)), function(i) {
    col <- paste0(substr(colpallette[groupnum[i]], 1, 7), "22")
    list(type = "rect", fillcolor = col, line = list(color = col), x0 = data_ints[1,i,1], x1 = data_ints[2,i,1], xref = "x", y0 = data_ints[1,i,2], y1 = data_ints[2,i,2], yref = "y")
  })
  colnames(data_ests) <- c('D1', 'D2')
  plotdata <- data.frame(data_ests, Group=groupfac, Hover=pointtext, row.names = NULL)
  output <- plot_ly(plotdata, x=~D1, y=~D2, color = ~Group, text=~Hover, type='scatter', mode='markers', colors = colpallette)
  if (rectangles=='full') {
    output <- layout(output, xaxis=list(title=paste0('D1 (',PoV[1],'%)')), yaxis = list(title=paste0('D2 (',PoV[2],'%)')) , shapes=rectlist)
  } else {
    output <- layout(output, xaxis=list(title=paste0('D1 (',PoV[1],'%)')), yaxis = list(title=paste0('D2 (',PoV[2],'%)')) )
  }
  output <- add_annotations(output, text=vartext, ax=var_ests[,1]*expansion, ay=var_ests[,2]*expansion, x=var_ests[,1], y=var_ests[,2], xref='x', yref='y', axref='x', ayref='y', arrowcolor=colsarrows)
  output <- add_annotations(output, text='', x=var_ests[,1], y=var_ests[,2], ax=rep(0,nArrows), ay=rep(0,nArrows), xref='x', yref='y', axref='x', ayref='y', arrowcolor=colsarrows)
  output <- config(output, toImageButtonOptions = list(format=export_format, filename=export_filename, width=export_width, height=export_height))
  output
}
```

## Latent variable model results

::: {.columns-2}
```{r, echo=FALSE, fig.width=3.5, fig.height=2.8}
BayesBiplot.plotly(sims$Z, sims$W, PoV = X, rotate = 0, grouping_points = Group)
BayesBiplot.plotly(sims$Z, sims$W, PoV = X, rotate = -0.4, grouping_points = Group)
```

```{r, echo=FALSE, fig.width=3.5, fig.height=2.8}
BayesBiplot.plotly(sims$Z, sims$W, PoV = stan_data$X, rotate = -0.2, grouping_points = Group)
BayesBiplot.plotly(sims$Z, sims$W, PoV = stan_data$X, rotate = -0.6, grouping_points = Group)
```
:::

## More rotation

```{r, echo=FALSE}
FullSpinBiplot <- function(data_sims, var_sims, grouping_points=rep("Data", dim(data_sims)[2]), vartext=NULL, pointtext=NULL, colpallette=c('rainbow', 'virShuffle', 'vir')[2], export_filename='plotlyexport', export_width=600, export_height=400, export_format='svg') {
  # Custom biplot function by Sean van der Merwe, University of the Free State
  # Takes simulated scores and simulated factors
  data_sims <- data_sims[,,1:2]
  var_sims <- var_sims[,,1:2]
  data_ests <- apply(data_sims,2:3,mean)
  var_ests <- apply(var_sims,2:3,mean)
  nGroups <- length(unique(grouping_points))
  mix_nums <- function(n) {  mix <- rev(order(abs((1:n)-(n+1)/2))); c(mix[n],mix[-n]) }
  nArrows <- nrow(var_ests)
  n <- nrow(data_ests)
  colsarrows <- viridisLite::plasma(nArrows); colsarrows <- colsarrows[mix_nums(nArrows)]
  if (is.null(pointtext)) {pointtext=rep('', n)}
  if (is.null(vartext)) {vartext=paste0('Var', seq_len(nArrows))}
  colnames(data_ests) <- c('D1', 'D2')
  colpallette <- switch(EXPR = colpallette, 
                        rainbow = {rainbow(nGroups)},
                          virShuffle = {pal <- viridisLite::viridis(nGroups+2)[c(-1,-(nGroups+2))]; pal[mix_nums(nGroups)]},
                          vir = {viridisLite::viridis(nGroups+2)[c(-1,-(nGroups+2))]}
                        )
  allplotdata <- data.frame(D1 = double(), D2 = double(), Group=character(), Hover=character(), Rotation=double())
  for (Degrees in seq(0, 345, 15)) {
    rotate <- Degrees*pi/180
    vmax <- matrix(c(cos(rotate), -sin(rotate), sin(rotate), cos(rotate)),2)
    rot_data_ests <- data_ests %*% vmax
    rot_var_ests <- var_ests %*% vmax
    maxes <- apply(rot_data_ests, 2, function(d) { max(abs(d)) })
    for (d in 1:2) { rot_var_ests[,d] <- rot_var_ests[,d]*maxes[d]/max(abs(rot_var_ests[,d]))/1.1 }
    colnames(rot_data_ests) <- c('D1', 'D2')
    plotdata <- data.frame(rot_data_ests, Group=grouping_points, Hover=pointtext, Rotation=rep(Degrees, n) , row.names = NULL)
    arrowacc <- 50
    for (v in 1:nArrows) {
      arrowdata <- data.frame(
        D1 = seq(0, rot_var_ests[v,1], l=arrowacc),
        D2 = seq(0, rot_var_ests[v,2], l=arrowacc),
        Group = rep(vartext[v], arrowacc),
        Hover = rep(vartext[v], arrowacc),
        Rotation = rep(Degrees, arrowacc),
        row.names = NULL )
      plotdata <- rbind(plotdata, arrowdata)
    }
    allplotdata <- rbind(allplotdata, plotdata)
  }

  output <- plot_ly(allplotdata, x=~D1, y=~D2, color = ~Group, text=~Hover, frame=~Rotation, type='scatter', mode='markers', colors = c(colpallette, colsarrows))
  output <- config(output, toImageButtonOptions = list(format=export_format, filename=export_filename, width=export_width, height=export_height))
  output
}
```

```{r, echo=FALSE, fig.width=8, fig.height=6}
FullSpinBiplot(sims$Z, sims$W, grouping_points = Group)
```

## Pause, ponder, and question {.larger}

::: {.green3}

-   Why go through all this trouble to get essentially the same results as we got from the simple PCA?

:::

# Combining

## Now for the magic

Now for the magic: what if I had multiple Likert scale responses and wanted to do a PCA style analysis?

```{r}
LikertScale <- c('Strongly Disagree', 'Disagree', 'Neutral', 'Agree', 'Strongly Agree')
nLikert <- length(LikertScale)
```

```{r}
set.seed(12345)
Y <- plogis(rawZ%*%true_weights)
Y <- matrix(rbinom(length(Y), nLikert-1, Y), nrow(X), ncol(X))
```

```{r, echo=FALSE}
plotly_table_UFS(Y)
```

## Multiple responses visual

```{r, echo=FALSE}
plotdata <- data.frame(Likert=factor(rep(LikertScale, times=full_Dim), levels = LikertScale), Question=factor(rep(seq_len(full_Dim), each=nLikert)), Responses = c(apply(Y, 2, table)))
fig <- plot_ly(plotdata, x = ~Question, y = ~Responses, type = 'bar', color = ~Likert, colors=unname(UFScolours[c(6,2,3,4,5)]) )
fig <- layout(fig, yaxis = list(title = 'Frequency'), barmode = 'stack')
fig
```

## Multiple Likert Scale Responses Model {.smaller}

```{stan, output.var='PCAbinomial', eval=FALSE}
// This Stan block defines a Bayesian PCA, Joseph H. Sakaya & Suleiman A. Khan but
// with the dependent variable changed to Binomial (by Sean van der Merwe, UFS)
data {
	int<lower=0> N; // Number of samples
	int<lower=0> D; // The original dimension
	int<lower=0> K; // The maximum latent dimension, must be < D
	int<lower=0> M; // The binomial maximum
	int<lower=0, upper=M> X[N*D]; // The data matrix as vector
}
parameters {
	matrix[N, K] Z; // The latent matrix
	matrix[D, K] W; // The weight matrix
	vector<lower=0>[K] alpha; // ARD prior
}
transformed parameters{
	vector<lower=0>[K] t_alpha;
  t_alpha = inv(sqrt(alpha));
}
model {
	to_vector(Z) ~ normal(0,1);
	alpha ~ gamma(1e-3,1e-3);				
	for(k in 1:K) W[,k] ~ normal(0, t_alpha[k]);
	X ~ binomial_logit(M, to_vector(Z*W'));
}
```

```{r, warning=FALSE, echo=FALSE}
stan_data<-list(N=nrow(Y),D=ncol(Y),K=4,X=c(Y),M=4)
```

```{r, warning=FALSE, eval=FALSE, echo=FALSE}
binfit <- sampling(PCAbinomial, stan_data, pars=c('alpha', 'W', 'Z'), iter = 4000, chains = mycores)
```

```{r, eval=FALSE, include=FALSE}
saveRDS(binfit, 'PCAexperimentBinomial.Rds')
```

```{r, include=FALSE}
binfit <- readRDS('PCAexperimentBinomial.Rds')
```

```{r, echo=FALSE}
sims <- extract(binfit)

nsims <- nrow(sims$alpha)
orders <- apply(sims$alpha, 1, order)
for (sim in seq_len(nsims)) {
  sims$W[sim,,] <- sims$W[sim,,orders[,sim]]
  sims$Z[sim,,] <- sims$Z[sim,,orders[,sim]]
  sims$alpha[sim,] <- sims$alpha[sim,orders[,sim]]
}
for (dim in seq_len(reduced_Dim)) {
  flip <- sims$W[,anchor,dim] < 0
  sims$W[flip,,dim] <- -sims$W[flip,,dim]
  sims$Z[flip,,dim] <- -sims$Z[flip,,dim]
}
```

## Proportion of variation {.smaller}

We calculate the proportion of variation in a crude way, by looking at the square differences of predicted and observed values.

We do this for 0, 1, and 2 dimensions and combine the numbers to get the marginal variation explained.

```{r}
ScaleMax <- 4
varhat <- sd(stan_data$X)^2
Zhat <- apply(sims$Z, 2:3, median)
What <- apply(sims$W, 2:3, median)
varhat1 <- sd((plogis(c(Zhat[,1] %*% t(What[,1])))*ScaleMax) - stan_data$X)^2
varhat2 <- sd((plogis(c(Zhat[,1:2] %*% t(What[,1:2])))*ScaleMax) - stan_data$X)^2
(PoV_Bayes <- c(varhat-varhat1, varhat1-varhat2)/varhat)
```

## Multiple Likert Scale Responses Results

```{r, echo=FALSE}
BayesBiplot.plotly(sims$Z, sims$W, PoV = round(PoV_Bayes*100,1), grouping_points = Group)
```

## Missing values

What if some of the values were missing (MCAR)?

```{r}
nMissing <- 50
nTotal <- stan_data$N*stan_data$D
whichMissing <- sample.int(nTotal, nMissing)
Y2 <- Y
Y2[whichMissing] <- NA_integer_
whichObs <- seq_len(nTotal)[-whichMissing]
```

```{r, echo=FALSE}
plotly_table_UFS(Y2)
```

## Missing data model {.smaller}

```{stan, output.var='PCAbinomialMissing', eval=FALSE}
// This Stan block defines a Bayesian PCA, Joseph H. Sakaya & Suleiman A. Khan but with the dependent variable changed to Binomial and missing values added by Sean  van der Merwe, UFS
data {
	int<lower=0> N; // Number of samples
	int<lower=0> D; // The original dimension
	int<lower=0> K; // The maximum latent dimension, must be < D
	int<lower=0> M; // The binomial maximum
	int<lower = 0> N_obs; // Total number of observed values in matrix
	int<lower=0, upper=M> X[N_obs]; // The data matrix as vector, only observed
	int<lower = 1, upper = N*D> ii_obs[N_obs]; }
  // Indicator vector of which values are observed in the matrix (column wise)
parameters {
	matrix[N, K] Z; // The latent matrix
	matrix[D, K] W; // The weight matrix
	vector<lower=0>[K] alpha; }// ARD prior
transformed parameters {
	vector<lower=0>[K] t_alpha;
	vector[N*D] zw;
  t_alpha = inv(sqrt(alpha));
  zw = to_vector(Z*W'); }
model {
	to_vector(Z) ~ normal(0,1);
	alpha ~ gamma(1e-3,1e-3);				
	for (k in 1:K) W[,k] ~ normal(0, t_alpha[k]);
	X ~ binomial_logit(M, zw[ii_obs]); }
```

## Missing data model fit

```{r}
stan_data <- list(N=nrow(Y2), D=ncol(Y2), K=4, X=c(Y2)[whichObs], M=4, N_obs=(nTotal-nMissing), ii_obs=whichObs)
```

```{r, warning=FALSE, eval=FALSE}
binfit <- sampling(PCAbinomialMissing, stan_data, pars=c('alpha', 'W', 'Z'), iter = 4000, chains = mycores)
```

```{r, eval=FALSE, include=FALSE}
saveRDS(binfit, 'PCAexperimentBinomialmissing.Rds')
```

```{r, include=FALSE}
binfit <- readRDS('PCAexperimentBinomialmissing.Rds')
```

```{r, echo=FALSE}
sims <- extract(binfit)
nsims <- nrow(sims$alpha)
orders <- apply(sims$alpha, 1, order)
for (sim in seq_len(nsims)) {
  sims$W[sim,,] <- sims$W[sim,,orders[,sim]]
  sims$Z[sim,,] <- sims$Z[sim,,orders[,sim]]
  sims$alpha[sim,] <- sims$alpha[sim,orders[,sim]]
}
for (dim in seq_len(reduced_Dim)) {
  flip <- sims$W[,anchor,dim] < 0
  sims$W[flip,,dim] <- -sims$W[flip,,dim]
  sims$Z[flip,,dim] <- -sims$Z[flip,,dim]
}
```

```{r}
varhat <- sd(stan_data$X)^2
Zhat <- apply(sims$Z, 2:3, median)
What <- apply(sims$W, 2:3, median)
varhat1 <- sd((plogis(c(Zhat[,1] %*% t(What[,1])))*ScaleMax)[whichObs] - stan_data$X)^2
varhat2 <- sd((plogis(c(Zhat[,1:2] %*% t(What[,1:2])))*ScaleMax)[whichObs] - stan_data$X)^2
PoV_Bayes <- c(varhat-varhat1, varhat1-varhat2)/varhat
```

## Missing data model results

```{r, echo=FALSE}
BayesBiplot.plotly(sims$Z, sims$W, PoV = round(PoV_Bayes[1:2]*100,1), grouping_points = Group)
```

## Missing values {.smaller}

The model is fitted using all of the observed values and none of the missing values. This is an ideal approach for data Missing Completely at Random (MCAR).

In the case of Missing at Random (MAR) and even Missing Not at Random (MNAR) where the missingness mechanism is based on the latent values, this approach is still probably a great approach, although I haven't tested it. For really complicated missingness, it may be possible to expand the models discussed to include additional external information, but this is beyond the scope of today's presentation.

**These models all naturally produce uncertain predictions of the missing values in the data set as a by-product, so multiple imputation is a breeze:**

```{r}
nImpute <- 10 # Number of new data sets to impute
SimsSelection <- sample(seq_len(nsims), nImpute) # Pick random set of simulations
NewDataSets <- lapply(SimsSelection, function(sim) {
  preddata <- rbinom(n*full_Dim, ScaleMax, plogis(c(sims$Z[sim,,] %*% t(sims$W[sim,,])))) # Predict all values
  preddata[whichObs] <- stan_data$X # Put back values we know for sure
  preddata <- matrix(preddata, n) # Put it in columns
  preddata # Return the filled in data set
})
# Write all data sets to an Excel file with each data set on a sheet
openxlsx::write.xlsx(NewDataSets, 'ImputedData.xlsx')
```

## Pause, ponder, question

Let us pause here a moment and think about what we have discussed. 

::: {.purple}
-   We have an algorithm that can easily produce a joint predictive posterior distribution over all observations, including missing values, regardless of the apparent marginal distributions.

    -   Multiple imputation of any number of plausible datasets is thus a breeze.

-   The only major assumption is the existence of an underlying latent structure of lower dimension that explains the observed variation via linear transformation.
:::

# Over the top fun time

## Let us go nuts with distributions

What if the first variable was normal, but then the second one had some extreme values, so maybe a $t$ density would fit better? What if the third variable was something counted, so maybe a Poisson distribution? The fourth variable can be a True/False question and the last one a Likert scale.

Let's do a latent variable analysis and see if we can still get the underlying latent dimensions out of the chaos. 😀

```{r, echo=FALSE}
set.seed(12345)
M <- rawZ%*%true_weights
M[,2] <- rt(n, 4) + M[,2]
M[,3] <- rpois(n, exp(M[,3]))
M[,4] <- rbinom(n, 1, plogis(M[,4]))
M[,5] <- rbinom(n, nLikert-1, plogis(M[,5]))
colnames(M) <- c('Normal', 't', 'Poisson', 'Bernoulli', 'Binomial')
plotly_table_UFS(round(M,3))
```

```{stan, output.var='PCAcrazy', eval=FALSE, echo=FALSE}
// This Stan block defines a Bayesian PCA, Joseph H. Sakaya & Suleiman A. Khan but
// with the dependent variable changed to Binomial (by Sean van der Merwe, UFS)
data {
	int<lower=0> N; // Number of samples
	int<lower=0> K; // The maximum latent dimension, must be < D
	real x1[N];
	real x2[N];
	int<lower=0> x3[N];
	int<lower=0, upper=1> x4[N];
	int<lower=0, upper=5> x5[N];
}
parameters {
	matrix[N, K] Z; // The latent matrix
	matrix[5, K] W; // The weight matrix
	vector<lower=0>[K] alpha; // ARD prior
	real<lower=0.5> nu;
	real<lower=0> s1;
	real<lower=0> s2;
}
transformed parameters{
	vector<lower=0>[K] t_alpha;
	vector[N*5] zw;
  t_alpha = inv(sqrt(alpha));
  zw = to_vector(Z*W');
}
model {
	to_vector(Z) ~ normal(0,1);
	alpha ~ gamma(1e-3,1e-3);				
	for (k in 1:K) W[,k] ~ normal(0, t_alpha[k]);
	x1 ~ normal(zw[1:N], s1);
	x2 ~ student_t(nu, zw[(N+1):(N*2)], s2);
	x3 ~ poisson_log(zw[(N*2+1):(N*3)]);
	x4 ~ binomial_logit(1, zw[(N*3+1):(N*4)]);
	x5 ~ binomial_logit(5, zw[(N*4+1):(N*5)]);
	nu ~ exponential(1);
	s1 ~ exponential(0.1);
	s2 ~ exponential(0.1);
}
```

```{r, warning=FALSE, echo=FALSE}
crazy_data<-list(N=nrow(M),K=4, x1=M[,1], x2=M[,2], x3=M[,3], x4=M[,4], x5=M[,5])
```

```{r, warning=FALSE, eval=FALSE, echo=FALSE}
crazyfit <- sampling(PCAcrazy, crazy_data, pars=c('alpha', 'W', 'Z'), iter = 4000, chains = mycores)
```

```{r, eval=FALSE, include=FALSE}
saveRDS(crazyfit, 'PCAexperimentCrazy.Rds')
```

```{r, include=FALSE}
crazyfit <- readRDS('PCAexperimentCrazy.Rds')
```

```{r, echo=FALSE}
Csims <- extract(crazyfit)
orders <- apply(Csims$alpha, 1, order)
for (sim in seq_len(nsims)) {
  Csims$W[sim,,] <- Csims$W[sim,,orders[,sim]]
  Csims$Z[sim,,] <- Csims$Z[sim,,orders[,sim]]
  Csims$alpha[sim,] <- Csims$alpha[sim,orders[,sim]]
}
for (dim in seq_len(reduced_Dim)) {
  flip <- Csims$W[,anchor,dim] < 0
  Csims$W[flip,,dim] <- -Csims$W[flip,,dim]
  Csims$Z[flip,,dim] <- -Csims$Z[flip,,dim]
}
```

```{r, echo=FALSE}
Zhat <- apply(Csims$Z, 2:3, median)
What <- apply(Csims$W, 2:3, median)
preds1 <- c(Zhat[,1] %*% t(What[,1]))
preds2 <- c(Zhat[,1:2] %*% t(What[,1:2]))
preds1[(n*2+1):(n*3)] <- exp(preds1[(n*2+1):(n*3)])
preds1[(n*3+1):(n*4)] <- plogis(preds1[(n*3+1):(n*4)])
preds1[(n*4+1):(n*5)] <- plogis(preds1[(n*4+1):(n*5)])*5
preds2[(n*2+1):(n*3)] <- exp(preds2[(n*2+1):(n*3)])
preds2[(n*3+1):(n*4)] <- plogis(preds2[(n*3+1):(n*4)])
preds2[(n*4+1):(n*5)] <- plogis(preds2[(n*4+1):(n*5)])*5
PoV_Bayes <- PoVfunc(M, preds1, preds2)
```

## Chaos results

```{r, echo=FALSE}
BayesBiplot.plotly(Csims$Z, Csims$W, PoV = round(PoV_Bayes*100,1), grouping_points = Group)
```

## Moooorrreee Poowwwweeerrrr!!!! {.larger}

We're going to stop the crazy here, but some things you could try at home:

-   Bring in structure to the underlying latent dimensions
-   Test whether the dimensions are related to an external factor
-   Predict new or other things using the latent dimensions

# The uncertainty problem

## So far we just took means

The plots are done on a best guess basis; so there are many kinds of uncertainty we are ignoring, such as:

::: {.red2}
-   The latent values are not certain (in any dimension)
-   The relative orientation of the variables is not certain
:::

::: {.green2}
-   The absolute orientations of the variables are not certain
-   The distributions are not certain
-   The simulation accuracy is not certain
-   The weights are not certain
-   The relevant dimensions are not certain
-   The order of the dimensions is not certain
-   The missing values are not certain
:::

## Incorporating half of the uncertainty

If we **give each point a rectangle indicating the shortest 50% interval** in each dimension then we get a lot of uncertainty

```{r, echo=FALSE}
BayesBiplot.plotly(Csims$Z, Csims$W, PoV = round(PoV_Bayes*100,1), rectangles = 'full', grouping_points = Group)
```

## Some uncertainty doesn't matter

::: {.green2}
-   The uncertainty in the dimension weights and orders are usually small relative to those weights, so the preferred order of at least the first few dimensions is usually quite clear

-   The relevant dimensions can be chosen by practicality (say 2 for a presentation like this one)

    -   You can also do a scree plot, but with the bonus that you could add confidence limits to the plot if you wanted to

-   The distributions have little impact as long as they are chosen sensibly

    -   Like using Negative Binomial instead of Poisson if you have over-dispersion

-   The simulation accuracy is reported by Stan and if it is too large you can just do more simulations
:::

## Orientation {.smaller}

We've already looked at flipping the dimensions using a variable as an anchor.

We can go further and also rotate and scale the plot using a single anchor, to further reduce the arbitrary uncertainty that doesn't matter.

```{r}
Z <- Csims$Z[,,1:2]
W <- Csims$W[,,1:2]
anchorPoint <- colMeans(W[,1,])
anchorAngle <- atan(anchorPoint[2]/anchorPoint[1])
for (sim in 2:nsims) {
  currentAngle <- atan(W[sim, 1, 2]/W[sim, 1, 1])
  rotMat <- matrix(c(cos(anchorAngle - currentAngle), -sin(anchorAngle - currentAngle), sin(anchorAngle - currentAngle), cos(anchorAngle - currentAngle)),2)
  W[sim,,] <- W[sim,,] %*% rotMat
  Z[sim,,] <- Z[sim,,] %*% rotMat
  currentPoint <- W[sim,1,1]
  W[sim,,] <- W[sim,,] * anchorPoint[1]/currentPoint
  Z[sim,,] <- Z[sim,,] * anchorPoint[1]/currentPoint
}
```

## Uncertainty in variables {.smaller}

Let us look at the angles of variables 1 and 2 across simulations before and after this transformation:

```{r, echo=FALSE}
Angles1before <- atan(Csims$W[,1,2]/Csims$W[,1,1])
Angles2before <- atan(Csims$W[,2,2]/Csims$W[,2,1])
Angles1after <- atan(W[,1,2]/W[,1,1])
Angles2after <- atan(W[,2,2]/W[,2,1])
```

::: {.columns-2}
```{r, fig.height=2.4, fig.width=3.6, echo=FALSE}
plot_ly(x=Angles1before, type='histogram', marker = list(color = UFScolours[2], line = list(color = UFScolours[3], width = 1)))
plot_ly(x=Angles2before, type='histogram', marker = list(color = UFScolours[2], line = list(color = UFScolours[3], width = 1)))
plot_ly(x=Angles1after, type='histogram', marker = list(color = UFScolours[2], line = list(color = UFScolours[3], width = 1)))
plot_ly(x=Angles2after, type='histogram', marker = list(color = UFScolours[2], line = list(color = UFScolours[3], width = 1)))
```
:::

## Fuzzy latent points again {.smaller}

Perhaps now the remaining uncertainty in the points might be easier to visualise, or perhaps the inherent uncertainty in trying to estimate what we cannot directly observe is just that large?

```{r, echo=FALSE, fig.height=4, fig.width=7}
Zhat <- apply(Z, 2:3, median)
What <- apply(W, 2:3, median)
preds1 <- c(Zhat[,1] %*% t(What[,1]))
preds2 <- c(Zhat[,1:2] %*% t(What[,1:2]))
preds1[(n*2+1):(n*3)] <- exp(preds1[(n*2+1):(n*3)])
preds1[(n*3+1):(n*4)] <- plogis(preds1[(n*3+1):(n*4)])
preds1[(n*4+1):(n*5)] <- plogis(preds1[(n*4+1):(n*5)])*5
preds2[(n*2+1):(n*3)] <- exp(preds2[(n*2+1):(n*3)])
preds2[(n*3+1):(n*4)] <- plogis(preds2[(n*3+1):(n*4)])
preds2[(n*4+1):(n*5)] <- plogis(preds2[(n*4+1):(n*5)])*5
PoV_Final <- PoVfunc(M, preds1, preds2)
BayesBiplot.plotly(Z, W, PoV = round(PoV_Final*100,1), rectangles = 'full', grouping_points = Group)
```

::: {.red2}
-   It seems that over-analysing a biplot is a quick road to spurious results.
:::

# Conclusion

## Thank you for listening {.larger}

-   Stan is the flexible all-purpose solution for people who find themselves with unusual statistical modelling problems.

-   It can do much more advanced stuff than what was shown here. See their website for more information.

-   The help files are extensive. The reference documents are the most useful when you are trying to solve a specific problem.

-   There are many helper packages that can do a lot of the work for you, including creating basic model code.

    -   The models you saw today were created manually though.

-   Fitting models with Stan can be a lot of fun because of the freedom it provides.
