---
title: "Analysis script for: Reproducible analysis of disease space via principal components: a brief tutorial and R package (syndRomics)"
author: "Torres-Espin A, Chou A, Huie JR, Kyritsis N, Upadhyayula PS, Ferguson AR"
date: "15/12/2020"
output: 
  html_document:
    toc: yes
    toc_float: yes
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = F, error = F, warning = F, cache = T, eval = T)
```

# Intro

This script contains the analysis of 2 case studies that demonstrate the use of syndromics analysis and the syndRomic package. Details of these can be seen in the paper: 

    Abel Torres-Espin, Austin Chou, J. Russell Huie, Nikos Kyritsis, Pavan Upadhyayula, and Adam R. Ferguson. Reproducible analysis of disease space via principal components: a brief tutorial and R package (syndRomics). 

**Libraries**

Libraries used in this notebook:

```{r}
library(syndRomics)
library(tidyverse)
library(mice)
library(MissMech)
library(naniar)
library(progress)
library(Gifi)
```


If these packages are not present, they can be installed with:

    install.packages("packageName")

syndRomics package can be installed from GitHub using the following code:

```{r, eval=FALSE}
install.packages('devtools')
devtools::install_github(repo = "ucsf-ferguson-lab/syndRomics@*release")
```

This will install the latest release. To install previous versions, you can use:

    devtools::install_github(repo = "ucsf-ferguson-lab/syndRomics@versionNumber")

For example:

```{r, eval=FALSE}
devtools::install_github(repo = "ucsf-ferguson-lab/syndRomics@0.0.1.9000")
```

# First case study

The data comes from Open Data Commons for Spinal Cord Injury (ODC-SCI) [dataset 26](https://scicrunch.org/odc-sci/about/odc-sci_26). The same data have been used in the past by Adam Ferguson team to perform SCI syndromics analysis:

    Ferguson AR, Irvine K-A, Gensel JC, Nielson JL, Lin A, Ly J, et al. Derivation of Multivariate Syndromic Outcome Metrics for Consistent Testing across Multiple Models of Cervical Spinal Cord Injury in Rats. PLOS ONE. 2013 Mar 27;8(3):e59712.

[10.1371/journal.pone.0059712](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0059712)

For details on the analysis, please see the manuscript.

## Data pre-processing

Loading the data and extracting the variables at 6 weeks after injury. To re-run the script, you will need to download the data and change the path in the code.

```{r}
# Specify the path and file name of the ODC-SCI:26 data from odc-sci.org
odc.df<-read.csv('data/odc-sci_26.csv', na.strings = '')
odc.df<-odc.df[-1,]

# select variables at week 6 and shorten the variable name
odc.df.filter<-odc.df[,str_detect(colnames(odc.df),
                                  "_wk6")]
colnames(odc.df.filter)<-str_remove(colnames(odc.df.filter), "_wk6")
colnames(odc.df.filter)<-str_remove(colnames(odc.df.filter), "AbsDev")
```

## Missing data analysis and imputation

We first visualize the missingess and the missing patterns.

```{r}
set.seed(500)

#Shadow plot
vis_miss(as.data.frame(odc.df.filter))+
  theme(plot.margin = unit(c(3,5,3,3), 'mm'),
        legend.position = 'right')

ggsave('Case1 missing shadow.png', width = 5, height = 4)

#Upset plot
upset<-naniar::gg_miss_upset(as.data.frame(odc.df.filter),nsets=20, mb.ratio=c(0.4,0.6))

upset

png('Case1 missing upset.png', units = 'in', width = 3.5, height = 3.5, res = 300)
  upset
dev.off()
```

Now we test whether the data is missing complete at random (MCAR).

```{r}
odc.df.filter<-as.matrix(apply(odc.df.filter, 2, as.numeric))
odc.df.filter%>%TestMCARNormality()
```

We impute missing data using multiple imputation (m=50) and we aggregate the imputations by averaging.

```{r}
# compute 50 imputations by PMM
odc.df.filter.mice<-mice(odc.df.filter, m=50, method='pmm', printFlag = F)
odc.df.multi.completed<-complete(odc.df.filter.mice, 'all', include = F)

# Aggregate the imputations by average
odc.df.multi.completed<-lapply(odc.df.multi.completed, as.matrix)
odc.df.completed<-apply(simplify2array(odc.df.multi.completed), 1:2, mean)

naniar::vis_miss(as.data.frame(odc.df.completed))
```

## PCA

Perform linear PCA with centered and unit variance scaled data
```{r}
pca_data<-odc.df.completed
pca<-prcomp(pca_data, center = T, scale. = T)

# calculate the standardizes loadings from the PCA solution using the package
original_loadings<-syndRomics::stand_loadings(pca, pca_data)

write.csv(original_loadings,'Case1 original loadings.csv')
```

## Component selection

We perform a permutation test with P 10000 to test for significance of VAF on the first 5 PCs. As a result, we select 3 components for subsequent analysis (see paper).

```{r}
per.VAF<-syndRomics::permut_pc_test(pca, pca_data, P=3000, ndim=5, statistic = 'VAF')

# Plot the original VAF vs. the permuted VAF
plot(per.VAF, plot_type = "VAF", ndim=1:5, plot_resample = T)

ggsave('Case1 permut VAF plot.png', width = 2.5, height = 3)

# Results of the permutation test
per.VAF$results

write.csv(per.VAF$results,'Case1 permut VAF results.csv')
```

## Component interpretation

Here we calculate the permutation test of the loadings using the permV methods (see paper for details on interpretation)
```{r}
# Loading permutation test on 3000 resamples of the selected first 3 PCs
per.load<-syndRomics::permut_pc_test(pca, pca_data, P=3000, ndim=3, 
                                     statistic = 's.loadings', perm.method = 'permV')

# Plot the results
plot(per.load, ndim=1:3, plot_resample = T, colors=c('grey','grey','grey'), gradient_color=F, plot_cutoff=F)+
  ggtitle('s.loadings. permut. test (3000 samples)')+
  xlab(NULL)+ylab(NULL)+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10))

ggsave('Case1 permuted loadings.png', width = 6.5, height = 3.5)

# Results table
per.load$results
write.csv(per.load$results, "Case1 permutation load p3000.csv")
```

We do the same with the communalities
```{r}
# Communalities permuation test with 3000 resamples
per_com<-permut_pc_test(pca, pca_data, P=3000, ndim=3, 
                                     statistic = 'commun', perm.method = 'permV')

# Plot the results
plot(per_com, plot_resample = T, ndim=1:3, plot_list_center=T, plot_legend=F, colors=c("grey","grey"), gradient_color=F)+
  ggtitle('Commun. permut. test (3000 samples) on 3 PCs')+
  xlab(NULL)+ylab(NULL)+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10))

ggsave('Case1 permuted communalities.png', width = 3.5, height = 3)

# Results table
per_com$results
write.csv(per_com$results, "Case1 permutation comm p3000.csv")

```

## Component stability

Bootstrap the loadings and calculate similarity metrics for the PCs using 3000 resamples with the balanced bootstrap.
```{r}
lb<-pc_stability(pca, pca_data, B = 3000, ndim=3, s_cut_off = 0.1, 
                              sim = 'balanced', barmap_plot = T)
#table of PC similarities
lb$PC_similarity

#Plot loadings
plot(lb, ndim = 1:3, plot_resample = T, gradient_color=F, colors=c("grey", "grey","grey"), plot_cutoff=F, plot_list_center=T)+
  ggtitle('Bootstrap loadings (3000 samples)')+
  xlab(NULL)+ylab(NULL)+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10))

ggsave('Case1 boot_loads.png', width = 5.2)

# Plot communalities
plot(lb, communalities = T,ndim = 1:3, plot_list_center=T,plot_resample = T, gradient_color=F, colors=c("grey", "grey","grey"))+
  ggtitle('Bootstrap communalities (3000 samples) on 3PCs')+
  xlab(NULL)+ylab(NULL)+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10))

ggsave('Case1 boot_commun.png', width = 3)

#Save similarity metrics
write.csv(lb$PC_similarity, 'Case1 boot similarity.csv')
```

### Stability of missing data

Run PCA and extract the loadings for each imputed dataset (m=50)
```{r}
missing_loadings<-lapply(odc.df.multi.completed, function(x){
  temp_data<-x
  temp_pca<-prcomp(temp_data, scale. = T, center = T)
  
  syndRomics:::extract_loadings(temp_pca, temp_data)
})
```

We use Procrustes rotation to rotate each loading matrix obtained from each of the 50 imputed datasets towards the parent matrix to deal with indeterminacies (see paper)
```{r}
parent<-syndRomics:::extract_loadings(pca, pca_data)
parent$Variables<-NULL

rotated_missing_loadings<-lapply(missing_loadings, function(x){
  x$Variables<-NULL
  r<-pracma::procrustes(as.matrix(parent)[,1:3],as.matrix(x)[,1:3])
  r$P
})
```

Now we extract the similarity metrics between the loadings of the missing datasets.
```{r}
sim_missing<-syndRomics::component_similarity(rotated_missing_loadings,
                                              ndim = 3, s_cut_off = 0.2)

sim_missing$index_mean
sim_missing$index_sd

write.csv(sim_missing, 'Case1 boot similarity missing.csv')
```

We can also plot the average loadings and their 95% confidence intervals
```{r}

#Build the loading dataset
load_list<-rotated_missing_loadings
load_list<-lapply(load_list, function(x){
  x[,1:3]
})

load_df<-as.data.frame(do.call(rbind, load_list))
colnames(load_df)<-paste0("PC", 1:ncol(load_df))
load_df$Variables<-rep(rownames(original_loadings),50)
load_df$id<-rep(1:50, each=nrow(original_loadings))

load_df<-load_df%>%pivot_longer(-c(id,Variables),
                              names_to = "component", values_to = "loading")

# Summarize the loadings by average and 95%CI
results<-load_df%>%
  group_by(Variables, component)%>%
  dplyr::summarise(mean=mean(loading), ci_low=quantile(loading, (1-0.95)/2, na.rm = T),
            ci_high=quantile(loading, 1-(1-0.95)/2, na.rm = T))

# Plot the loadings and the CI using the barmap_loading() function from the package
syndRomics::barmap_loading(pca, pca_data, ndim=1:3, resample_ci = results,plot_cutoff = F,
               plot_list_center = T, gradient_color=F, colors=c("grey", "grey","grey"))+
  ggtitle("Multiple imputation (m=50) loading stability")+
  xlab(NULL)+ylab(NULL)+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10))

ggsave('barmap_missing.png', height = 4.01, width = 7.5)
```


## Component visualization

Generating and saving syndRomic plots with cutoff 0.45 and heatmaps and barmap plots with cutoff 0.21, 0.25 and 0.4 for PC1, PC2 and PC3, respectively (see paper for more details)

```{r, fig.height=9, fig.width=9}
# Extract syndromic plots
splot<-syndromic_plot(pca, pca_data, cutoff = 0.45)
splot$PC1
splot$PC2
splot$PC3

# Save the plots in PDF
ggsave(plot = splot$PC1,filename = 'Case1 splot PC1.pdf', width = 9, height = 9)
ggsave(plot = splot$PC2,filename = 'Case1 splot PC2.pdf', width = 9, height = 9)
ggsave(plot = splot$PC3,filename = 'Case1 splot PC3.pdf', width = 9, height = 9)
```

```{r }
#Heatmap plot
heatmap_loading(pca, pca_data, ndim=1:3, cutoff = c(0.21,0.25,0.4), star_values = T, text_values = F)+
  ggtitle("Standardized loadings")+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10))

ggsave('Case1 heatmap.png',height = 4.01, width = 3.2)
```

# Second case study

The data is a selected subset of variables of the TRACK-TBI PILOT dataset that was used and shared by the publication: 

    Nielson JL, Cooper SR, Yue JK, Sorani MD, Inoue T, Yuh EL, Mukherjee P, Petrossian TC, Paquette J, Lum PY, Carlsson GE, Vassar MJ, Lingsma HF, Gordon WA, Valadka AB, Okonkwo DO, Manley GT, Ferguson AR; TRACK-TBI Investigators. Uncovering precision phenotype-biomarker associations in traumatic brain injury using topological data analysis. PLoS One. 2017 Mar 3;12(3):e0169490. doi: 10.1371/journal.pone.0169490. PMID: 28257413; PMCID: PMC5336356.
    
[10.1371/journal.pone.0169490](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0169490)
  
For details on the analysis, please see the manuscript.

```{r}
# Specify the path and file name of the data from 10.1371/journal.pone.0169490.s010
plos_df<-read.csv("data/journal.pone.0169490.s010.csv")
```

## Data pre-processing

We change the variable names for clarity. The data dictionary containing information for this dataset can be found at:

[10.1371/journal.pone.0169490](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0169490)

```{r}
#Rename the variables for easy use
plos_df<-plos_df%>%
  rename(CT_Marshall=Marshall, CT_Rotterdam=Rotterdam,
         CT_brain_pathology=CT_INTRACRANIAL_FINAL,
         CT_skull_FX=SKULLFX, CT_skullbase_FX=SKULLBASEFX,
         CT_facial_FX=FACIALFX, CT_EDH=EDH_FINAL, CT_SDH=SDH_FINAL,
         CT_SAH=SAH_FINAL, CT_contusion=CONTUSION_FINAL,
         CT_midlineshift=MIDLINESHIFT_FINAL,
         CT_cisterncomp=CISTERNCOMP_FINAL,
         PTSD_diagnosis_6mo=PTSD_6mo,
         GOSE_3mo=GOSE_OverallScore3M, GOSE_6mo=GOSE_OverallScore6M,
         WAIS_PSI_6mo=WAIS_PSI_Composite_6mo,
         CVLT_short_6mo=CVLTShortDelayCuedRecallStandardScore_6mo,
         CVLT_long_6mo=CVLTLongDelayCuedRecallStandardScore_6mo,
         SNP_COMT=rs4680, SNP_DRD2=rs6277, SNP_PARP1=rs3219119, 
         SNP_ANKK1_Gly318Arg=rs11604671, SNP_ANKK1_Gly442Arg=rs4938016,
         SNP_ANKK1_Glu713Lys=rs1800497)%>%
  relocate(PTSD_diagnosis_6mo, .before=GOSE_3mo)
```

## Missing data analysis
```{r}
set.seed(500)

# Shadow missing plot
vis_miss(plos_df[,-1])+
  theme(plot.margin = unit(c(3,5,3,3), 'mm'),
        legend.position = 'right')

ggsave('Case2 missing shadow.png', width = 7, height = 4)

# Upset plot
png('Case2 missing upset.png', units = 'in', width = 7, height = 3.5, res = 300)
  naniar::gg_miss_upset(plos_df,nsets=20, mb.ratio=c(0.4,0.6))
dev.off()
```

The MR_result variable has a lot of missigness, and we will drop it for this analysis. We also subtract the patient identifier.

```{r}
plos_df<-plos_df%>%select(-PatientNum, -MR_result)
```

We transform the variables into their type for further analysis. Specifically, all CT and SNP, as well as PTSD and GOSE are categorical and we represent them as a factor in R

```{r}
factor_vars<-colnames(plos_df)[!colnames(plos_df)%in%c("PatientNum", "CVLT_short_6mo","CVLT_long_6mo","WAIS_PSI_6mo")]

plos_df<-plos_df%>%
  mutate_at(.vars = factor_vars, .funs=factor)%>%
  mutate(WAIS_PSI_6mo=as.numeric(WAIS_PSI_6mo))
```

Now we test whether the data is missing complete at random (MCAR).

```{r}
plos_df<-as.matrix(apply(plos_df, 2, as.numeric))
plos_df%>% 
  TestMCARNormality()
```

We impute missing data using multiple imputation (m=50) and we aggregate the imputations by averaging the numeric variables or taking the mode of the categorical variables.

```{r}
plos_df.mice<-mice(plos_df, m=50, method='pmm', printFlag = F)
plos_df.completed<-complete(plos_df.mice, 'all', include = F)

# Aggregate the dataset to the closest integer to the median
nlpca_data<-lapply(plos_df.completed, as.matrix)
nlpca_data<-apply(simplify2array(nlpca_data), 1:2, function(x){
  round(median(x))
})
```

## nonlinear PCA

Perform nonlinear PCA using *princals()*. Variables are not restricted to ordinal. A three degree polynomial is specified for the continuous transformation function for all variables.

For details on using *princals()*, we refer to:

    De Leeuw, J., Mair, P., Groenen, P. J. F. (2017). Multivariate Analysis with Optimal Scaling.

and the function documentation *?princals*

```{r}
# NLPCA with princals
nlpca<-princals(nlpca_data, ndim = ncol(nlpca_data),
                ordinal = FALSE,
                degrees = 3, knots = knotsGifi(nlpca_data, type="E"))
barmap_loading(nlpca, nlpca_data, ndim=1:5)

par(mfrow=c(1,1))
plot(nlpca, plot.type = 'screeplot')


# calculate the standardizes loadings from the PCA solution using the package
plos_df_original_loadings<-syndRomics::stand_loadings(nlpca, nlpca_data)

write.csv(plos_df_original_loadings,'Case2 original loadings.csv')
```

## Component selection

We perform a permutation test with P 500 to test for significance of VAF on the first 10 PCs. As a result, we select 3 components for subsequent analysis (see paper).

```{r}
nlpca_per_VAF<-permut_pc_test(nlpca, nlpca_data, P=500, ndim=10, statistic = 'VAF')

# Plot the original VAF vs. the permuted VAF
plot(nlpca_per_VAF, ndim=1:8, plot_resample = T)
ggsave('Case2 permut VAF plot plos.png', width = 3, height = 4)

nlpca_per_VAF$results
write.csv(nlpca_per_VAF$results,'Case2 permut VAF results plos.csv')

save(nlpca_per_VAF,file="nlpca_per_VAF_plos.Rda")
```

## Component interpretation

Here we calculate the permutation test of the loadings using the permV methods (see paper for details on interpretation)

This might take several minutes to compute
```{r}
nlpca_per_load<-syndRomics::permut_pc_test(nlpca, nlpca_data, P=500, ndim=3, 
                                     statistic = 's.loadings', perm.method = 'permV')

plot(nlpca_per_load, ndim=1:3, plot_resample = T, colors=c("grey","grey","grey"), gradient=F, cutoff=0)+
  xlab(NULL)+ylab(NULL)+
  ggtitle('Permutation loadings test (500 samples)')+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10, hjust = 0.5))

ggsave('Case2 permuted loadings, plos.png', width = 4.9)

nlpca_per_load$results
write.csv(nlpca_per_load$results, "Case2 permutation load nlpca.csv")
```

## Component stability

We bootstrap the loadings and calculate similarity metrics for the PCs using 1000 resamples with the balanced bootstrap.
```{r}
lb_nlpca_plos<-syndRomics::pc_stability(nlpca, nlpca_data, B = 1000, 
                                        ndim=3, s_cut_off = 0.1, sim = 'balanced')
lb_nlpca_plos$PC_similarity

#Plot the loadings
plot(lb_nlpca_plos, ndim = 1:3, plot_resample = T, gradient_color=F, colors=c("grey", "grey","grey"),plot_cutoff=F)+
  ggtitle('Bootstrap loadings (1000 samples)')+
  xlab(NULL)+ylab(NULL)+
  theme(plot.title = element_text(size = 10,hjust = 0.5), plot.title.position = "plot")

ggsave('Case2 boot_loads_plos.png', width = 4.9)

#Plot the communalities
plot(lb_nlpca_plos, communalities = T,ndim = 1:3, plot_list_center=T,plot_resample = T, gradient_color=F, colors=c("grey", "grey","grey"))+
  ggtitle('Bootstrap communalities (1000 samples) on 3 PCs')+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10))

ggsave('Case2 boot_commun_plos.png', width = 3.5)
write.csv(lb_nlpca_plos$PC_similarity, 'Case2 boot similarity_plos.csv')
```

### Stability of missing data

Helper function to run PCA and extract the loadings for each imputed dataset (m=50)
```{r}
missing_loadings_plos<-lapply(plos_df.completed, function(x){
  temp_pca<-princals(x, ndim = ncol(x),
                ordinal = FALSE,
                degrees = 3, knots = knotsGifi(x,"E"))
  
  syndRomics:::extract_loadings(temp_pca, x)
})
```

We have to rotate each loading matrix towards the parent matrix to deal with indeterminacies (see paper)
```{r}
parent<-syndRomics:::stand_loadings(nlpca, nlpca_data)

rotated_missing_loadings_plos<-lapply(missing_loadings_plos, function(x){
  x$Variables<-NULL
  r<-pracma::procrustes(as.matrix(parent)[,1:3],as.matrix(x)[,1:3])
  r$P
})
```

Now we extract the similarity metrics between the loadings of the missing datasets.
```{r}
sim_missing_nlpca<-component_similarity(rotated_missing_loadings_plos,
                                        ndim = 3, s_cut_off = 0.2)

sim_missing_nlpca$index_mean
sim_missing_nlpca$index_sd

write.csv(sim_missing_nlpca, 'Case2 boot similarity missing.csv')
```

We can also plot the average loadings and their 95% confidence intervals
```{r}

# Create the loadings data for ploting
load_list_plos<-rotated_missing_loadings_plos
load_list_plos<-lapply(load_list_plos, function(x){
  x[,1:3]
})

load_df_plos<-as.data.frame(do.call(rbind, load_list_plos))
colnames(load_df_plos)<-paste0("PC", 1:ncol(load_df_plos))
load_df_plos$Variables<-rep(colnames(nlpca_data),50)
load_df_plos$id<-rep(1:50, each=ncol(nlpca_data))

load_df_plos<-load_df_plos%>%pivot_longer(-c(id,Variables),
                              names_to = "component", values_to = "loading")

# Summarize the loadings with the average and 95%CI
results<-load_df_plos%>%
  group_by(Variables, component)%>%
  dplyr::summarise(mean=mean(loading), 
                   ci_low=quantile(loading, (1-0.95)/2, na.rm = T),
                   ci_high=quantile(loading, 1-(1-0.95)/2, na.rm = T))

# Plot the loadings and the 95%CI
barmap_loading(nlpca, nlpca_data, ndim=1:3, resample_ci = results,plot_list_center = T,colors=c("grey","grey","grey"), gradient=F, cutoff=0)+
  ggtitle("Multiple imputation (m=50) loading stability")+
  xlab(NULL)+ylab(NULL)+
  theme(plot.title.position = 'plot', plot.title = element_text(size = 10, hjust = 0.5))

ggsave('Case2 barmap_missing_plos.png', height = 4.01, width = 7.5)
```

## Component visualization
```{r}
#Extract and save the syndromic plots
nlpca_syndromics<-syndromic_plot(nlpca, nlpca_data, cutoff = 0.4)

nlpca_syndromics$PC1
nlpca_syndromics$PC2
nlpca_syndromics$PC3

ggsave(plot = nlpca_syndromics$PC1,filename = 'Case2 splot PC1 plos.pdf', 
       width = 9, height = 9)
ggsave(plot = nlpca_syndromics$PC2,filename = 'Case2 splot PC2 plos.pdf', 
       width = 9, height = 9)
ggsave(plot = nlpca_syndromics$PC3,filename = 'Case2 splot PC3 plos.pdf', 
       width = 9, height = 9)
```

# Apendix

## Speed and precision

This script measures the computation time and the precision on the estimates of the permutation and the bootstrapping algorithms with 2 different versions of the first use case datasets. The analyzed version is considered a "small" size dataset.  The "small" dataset is repeated several times to artificially increase the volume of the data and is regarded as the "large" dataset.

THIS CODE CAN TAKE A WHILE TO RUN. The script is not run at the time of rendering this file. Please make sure to set eval=TRUE in the following code chunks to run the script. Alternatively, the script can be directly evaluated in RStudio.

```{r, eval=F}
#Helper function to run the test
test_speed<-function(pca, pca_data, test, samples, k, var_id, perm="permD"){
  
  # Function taking the results of a pca and the data used on the pca to run k
  # permutation test through permut_pc_test().
  # 
  # pca: pca result from prcomp function
  # pca_data: data passed to the prcomp function to run the pca
  # samples:  numeric vector of number of permutation samples to run
  # test: Either "permutation" or "bootstrap"
  # k:  number of replicates per permutation sample
  # var_id: variable number to return the results from
  # perm: permutation method to test. Either "permD" or "permV"
  
  cat(test, "k:",k, "var_id:", var_id, "/n")
  
  #vectors to save the results on
  results_time<-vector(mode='numeric',length = k)
  results_mean<-vector(mode='numeric',length = k)
  results_m_list<-list()
  results_t_list<-list()
  pb<-progress_bar$new(total = k*10)
  
  
  #' Iterate k times for each sample to produce computation time and confident interval on the estimates
  
  for (j in 1:k){
    set.seed(j)
    for (i in 1:length(samples)){
      
      if (test=="permutation"){
        suppressMessages(s<-system.time(
        p<-permut_pc_test(pca, pca_data, ndim = 10, P=samples[i],
                          statistic="s.loading", perm.method = perm)))
      }else if (test=="bootstrap"){
        suppressMessages(s<-system.time(
        p<-pc_stability(pca, pca_data, ndim = 3, communalities = F,
                        test_similarity = F,ci_type = "perc",B = samples[i])))
      }
      
      results_time[i]<-s[3]
      results_mean[i]<-p$results$mean[var_id]
      
      pb$tick()
    }
    results_t_list[[j]]<-results_time
    results_m_list[[j]]<-results_mean
  }
  
  return(list(results_t_list, results_m_list))
}

#number of resamples to test
samples<-c(10,25,50,100,250,500,750,1000,1500, 2000)

# construct the large dataset
pca_data_large<-matrix(rep(odc.df.completed,30), ncol = 54)
colnames(pca_data_large)<-rep(colnames(odc.df.filter),3)

# Run the PCA on the large dataset
pca_large<-prcomp(pca_data_large, center = T, scale. = T)
```

### Permutation

```{r, eval=FALSE}
# Run the different speed and precision test
small_permut_test_7<-test_speed(pca, pca_data, test="permutation", samples, 10, 7)
small_permut_test_5<-test_speed(pca, pca_data, test="permutation", samples, 10, 5)
large_permut_test_7<-test_speed(pca_large, pca_data_large, test="permutation",
                                samples, 10, 7)
large_permut_test_5<-test_speed(pca_large, pca_data_large, test="permutation",
                                samples, 10, 5)
```

**Plot speed and precision test**
```{r, eval=FALSE}
# Create the datafrem with the results of the test for ploting
df.plot_perm<-data.frame(mean=c(unlist(small_permut_test_7[[2]]),
                                  unlist(small_permut_test_5[[2]]),
                                  unlist(large_permut_test_7[[2]]),
                                  unlist(large_permut_test_5[[2]])),
                    n_sample=rep(samples, 40),
                    time=c(unlist(small_permut_test_7[[1]]),
                                  unlist(small_permut_test_5[[1]]),
                                  unlist(large_permut_test_7[[1]]),
                                  unlist(large_permut_test_5[[1]])),
                    df_size=c(rep('159x18', 200),rep('1590x54',200)),
                    loading_size=c(rep("High loading", 100), rep("Low loading", 100),
                                   rep("High loading", 100), rep("Low loading", 100)))

df.plot_perm$df_size<-factor(df.plot_perm$df_size, levels=c("159x18", "1590x54"))

df.plot_perm<-df.plot_perm%>%
  group_by(loading_size)%>%
  mutate(center_mean=scale(mean, center = T, scale = F))

ggplot(df.plot_perm, aes(log10(n_sample),mean))+
  stat_summary(fun = 'mean', geom='line', color='red')+
  stat_summary(fun = 'mean', geom='point', color='red')+
  stat_summary(fun.data = mean_cl_normal, geom='errorbar', color='red')+
  theme_minimal(base_size = 9)+
  xlab('P size (log10)')+
  ylab('Var 1 loading in PC1')+
  geom_hline(yintercept = 0, linetype="dashed")+
  facet_grid(df_size~loading_size)

ggsave('sample size p perm.png',height = 3, width = 3)

ggplot(df.plot_perm, aes(n_sample, time, color=df_size, group=df_size))+
  stat_summary(fun= 'mean', geom='line')+
  geom_jitter(alpha=0.3)+
  scale_color_manual(values=c('steelblue', 'orange'))+
  theme_minimal(base_size = 9)+
  xlab('P size')+
  ylab('Computing time (s)')+
  theme(legend.position = c(0.3,0.85))+
  labs(color='Data size')

ggsave('sample size P time perm.png', height = 3, width = 1.7)

# Create the dataframe with the confidence interval and its width for ploting
df.plot.CI_perm<-df.plot_perm%>%
  group_by(n_sample, df_size, loading_size)%>%
  summarise(sd=sd(mean),
            mean=mean(abs(mean)),
            ci_low=mean-2.26*sd,
            ci_high=mean+2.26*sd,
            per_CIwidth=abs(ci_high-ci_low))

ggplot(df.plot.CI_perm, aes(n_sample,per_CIwidth, color=loading_size))+
  geom_line()+
  geom_point()+
  theme_minimal(base_size = 9)+
  scale_color_manual(values=c('green3', 'purple1'))+
  xlab('P size')+
  ylab('95% Confidence Interval width')+
  facet_grid(~df_size)+
  theme(legend.position = c(0.8,0.8))+
  labs(color = "Loading size")

ggsave('sample CI width perm.png',height = 3, width = 3)
```

### Bootstrapping

```{r, eval=FALSE}
small_boot_test_7<-test_speed(pca, pca_data, test="bootstrap", samples, 10, 7)
small_boot_test_5<-test_speed(pca, pca_data, test="bootstrap", samples, 10, 5)
large_boot_test_7<-test_speed(pca_large, pca_data_large, test="bootstrap",
                              samples, 10, 7)
large_boot_test_5<-test_speed(pca_large, pca_data_large, test="bootstrap",
                              samples, 10, 5)
```

**Plot speed and precision test**
```{r, eval=FALSE}

# Create the datafrem with the results of the test for ploting
df.plot<-data.frame(mean=c(unlist(small_boot_test_7[[2]]),
                                  unlist(small_boot_test_5[[2]]),
                                  unlist(large_boot_test_7[[2]]),
                                  unlist(large_boot_test_5[[2]])),
                    n_sample=rep(samples, 40),
                    time=c(unlist(small_boot_test_7[[1]]),
                                  unlist(small_boot_test_5[[1]]),
                                  unlist(large_boot_test_7[[1]]),
                                  unlist(large_boot_test_5[[1]])),
                    df_size=c(rep('159x18', 200),rep('1590x54',200)),
                    loading_size=c(rep("High loading", 100), rep("Low loading", 100),
                                   rep("High loading", 100), rep("Low loading", 100)))

df.plot$df_size<-factor(df.plot$df_size, levels=c("159x18", "1590x54"))

df.plot<-df.plot%>%
  group_by(loading_size)%>%
  mutate(center_mean=scale(mean, center = T, scale = F))

ggplot(df.plot, aes(log10(n_sample),center_mean))+
  stat_summary(fun = 'mean', geom='line', color='red')+
  stat_summary(fun = 'mean', geom='point', color='red')+
  stat_summary(fun.data = mean_cl_normal, geom='errorbar', color='red')+
  theme_minimal(base_size = 9)+
  xlab('B size (log10)')+
  ylab('Var 1 loading in PC1 (centered to mean)')+
  facet_grid(df_size~loading_size)+
  geom_hline(yintercept = 0, linetype="dashed")

ggsave('sample size B.png',height = 3, width = 3)

ggplot(df.plot, aes(n_sample, time, color=df_size, group=df_size))+
  stat_summary(fun= 'mean', geom='line')+
  geom_jitter(alpha=0.3)+
  scale_color_manual(values=c('steelblue', 'orange'))+
  theme_minimal(base_size = 9)+
  xlab('B size')+
  ylab('Computing time (s)')+
  theme(legend.position = c(0.32,0.85))+
  labs(color='Data size')

ggsave('sample size B time.png', height = 3, width = 1.7)

# Create the dataframe with the confidence interval and its width for ploting
df.plot.CI<-df.plot%>%
  group_by(n_sample, df_size, loading_size)%>%
  summarise(sd=sd(mean),
            mean=mean(abs(mean)),
            ci_low=mean-2.26*sd,
            ci_high=mean+2.26*sd,
            per_CIwidth=abs(ci_high-ci_low))

ggplot(df.plot.CI, aes(n_sample,per_CIwidth, color=loading_size))+
  geom_line()+
  geom_point()+
  theme_minimal(base_size = 9)+
  scale_color_manual(values=c('green3', 'purple1'))+
  xlab('B size')+
  ylab('95% Confidence Interval width')+
  facet_grid(~df_size)+
  theme(legend.position = c(0.8,0.8))+
  labs(color = "Loading size")

ggsave('sample CI width B.png',height = 3, width = 3)
```

## Session information
```{r}
sessionInfo()
```

