ZICRM-ASCA_plus.Rmd

---
title: "Zero-inflated GLMM ASCA Simulation"
author: "Auke Haver"
date: "15/05/2022"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
source("R/ZICRM-ASCA_plus_supplementary_functions.R",local=TRUE)
```

# Introduction
In this simulation, we will apply the ASCA framework on unbalanced, zero-inflated, Poisson-distributed data. Furthermore, we will consider repeated measurements on unbalanced groups. This will require the use of Generalized Linear Mixed Models (GLMM) for estimating coefficients which can be used to construct effect matrices, which can be decomposed into scores and loadings through Singular Value Decomposition (SVD). We will utilize a 7-fold Jackknife approach for validation of our method.  


# Data Simulation
## Design
First, the parameters of the study need to be determined. Using the function gen_sim_info(), we simulate a design info file which constitutes a list with information on the design. Then, we construct a design matrix from the design info with sum-coded variables. Lastly, we construct a fixed effect model matrix for later multiplication with the fixed effect parameter matrix.
```{r Simulation Design}
# design size
design_info <- gen_sim_info(n_treatments       = 5,
                            n_timepoints       = 3,
                            timepoints         = c("0","6","12") %>% 
                              factor(x      =.,
                                     levels =.),
                            treatments         = c("1","2","3","4","5") %>% 
                              factor(x      =.,
                                     levels =.),
                            n_PC               = 2 ,
                            n_variables        = 16,
                            n_subjects         = 208,
                            n_dummy_replicates = 300,
                            family             = "nbinom", # choose pois, zipois, nbinom, zibinom
                            prob_zero          = 0.50,
                            sigma               = 1,
                            jackknife_fold     = 7) 

# Generate Design Matrix
design_matrix <- gen_sim_design(design_info)

# Construct a fixed effect model matrix
X <- model.matrix(~timepoint+timepoint:treatment,
                  design_matrix)


output_dir <- paste0(
  "Result/",
  "ZICRM-ASCA_plus_",
  design_info$family,
  "_",
  Sys.Date(),
  "_",
  design_info$prob_zero,
  "_",
  design_info$sigma,
  "/"
)
# Prepare output directory
for(dir in c("Result", output_dir)){
  while(!file.exists(dir)){
    dir.create(dir)
  }
}

```

## Scores and Loadings
Having decided the design of the simulation, we will continue with setting the scores ($T$) and loadings ($P$) for the time ($\alpha$) and time-treatment interaction ($\alpha\beta$) effects, in order to create effect matrices. The requirement for the scores and loadings are that the former must all be orthogonal and the latter must all be orthonormal.
```{r Simulate scores and loadings}
# SCORES
## Time Effect
Ta <- c(-4.00, 0.00, 4.00,      # PC 1
        1.00,-2.00, 1.00) %>%  # PC 2
  rep(each = design_info$n_treatments*design_info$n_dummy_replicates) %>% 
  matrix(ncol=2) 
## Time:Treatment Interaction Effect
Tab<-c(0.00, 0.02,-0.01,-0.02, 0.01, 3.00, 1.00,-0.50,-1.50,-2.00, 4.00,-0.50,-1.00,-1.50,-1.00,    # PC 1
       4.00,-0.50,-1.50,-0.50,-1.50, 0.40,-0.70, 0.00, 0.20, 0.10, 0.01,-0.09, 0.06, 0.01, 0.01)%>% # PC 2
  rep(each = design_info$n_dummy_replicates)%>% 
  matrix(ncol=design_info$n_PC) 

# LOADINGS
# Time Effect
Pa <- cbind(
  unit_vector(
    c(-1,-1,-1, 2, 0, 0, 0, 0, 1, 1,-2, 1, 0, 0, 0, 0)),
  unit_vector(
    c(-2, 2,-2,-1, 0, 0, 0, 0, 1,-1, 1, 2, 0, 0, 0, 0)))
# Time:Treatment Interaction Effect
Pab <- cbind(
  unit_vector(
    c( 0, 0, 0, 0,-1, 1, 2, 1,-2,-1,-1, 1, 0, 0, 0, 0)),
  unit_vector(
    c( 0, 0, 0, 0,-2, 1,-1,-1,-1, 2, 1, 1, 0, 0, 0, 0)))

# Explained variance per component
expl_var_a <-  perform_pca(Ta %*% t(Pa))$singular_values %>% 
  expl_var_from_svd()
expl_var_ab<-  perform_pca(Tab%*%t(Pab))$singular_values %>% 
  expl_var_from_svd()
```

## Plot Scores and Loadings
We are going to create 2 plots with each 4 subplots. For the Time and Time-Treatment Interaction effect we will make a scoreplot and loadings plot for principal component 1 and 2. 
```{r}
### TIME EFFECT ###
ggarrange(
  # SCORES PC1
  design_matrix %>% 
    # Modify design matrix
    mutate(score = Ta[,1])%>%
    # Initialize Subplot
    ggplot(aes(x=timepoint,
               y=score))+
    # Geometric objects
    geom_point(size=5)+
    geom_line(aes(x=timepoint,
                  y=score))+
    # Theme
    theme_classic(base_family="Times")+
    theme(plot.title = element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Labels
    labs(x="Timepoint",
         y=paste0("Scores PC1 (",
                  expl_var_a[1],
                  "%)"),
         color="Treatment Group"),
  # LOADINGS PC1
  tibble(var = LETTERS[1:nrow(Pa)],
         loading = Pa[,1])%>%
    # Initialize Subplot
    ggplot(aes(x=var,
               y=loading))+
    # Geometric objects
    geom_bar(stat="identity",
             fill = c(rep(dark_purple,4),
                      rep(teal,4),
                      rep(dark_purple,4),
                      rep(teal,4)))+
    geom_point()+
    # Theme
    theme_classic(base_family="Times")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Labels
    labs(x = "Variable",
         y = paste0("Loadings PC1 (",
                    expl_var_a[1],
                    "%)")),
  # SCORES PC2
  design_matrix %>% 
    # Modify design matrix
    mutate(score = Ta[,2])%>%
    # Initialize Subplot
    ggplot(aes(x=timepoint,
               y=score))+
    # Geometric objects
    geom_point(size=5)+
    geom_line(aes(x=timepoint,
                  y=score))+
    # Theme
    theme_classic(base_family="Times")+
    theme(plot.title = element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Labels
    labs(x="Timepoint",
         y=paste0("Scores PC2 (",
                  expl_var_a[2],
                  "%)")),
  # LOADINGS PC2
  tibble(var = LETTERS[1:nrow(Pa)],
         loading = Pa[,2])%>%
    # Initialize Subplot
    ggplot(aes(x=var,
               y=loading))+
    # Geometric objects
    geom_bar(stat="identity",
             fill = c(rep(dark_purple,4),
                      rep(teal,4),
                      rep(dark_purple,4),
                      rep(teal,4)))+
    geom_point()+
    # Theme
    theme_classic(base_family="Times")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Labels
    labs(x = "Variable",
         y = paste0("Loadings PC2 (",
                    expl_var_a[2],
                    "%)")),
  # Plot settings
  ncol=2,
  nrow=2,
  align="hv") %>%
  # Add title
  annotate_figure(top=text_grob("Time Effect (A)",
                                size=30,
                                family="Times")) %>%
  # Export to pdf
  ggexport(
    filename=paste0(
      output_dir,
      "Time_Simulated.pdf"),
    plot=.,
    width = 8,
    height= 8)

# Time-Treatment Interaction
ggarrange(
  # SCORES PC1
  design_matrix %>% 
    # Modify design matrix
    mutate(score = Tab[,1])%>%
    # Initialize plot
    ggplot(aes(x=timepoint,
               y=score,
               color=treatment))+
    # Geometric objects
    geom_point(size=5)+
    geom_line(aes(group=treatment),
              size=1)+
    # Theme
    theme_classic(base_family="Times")+
    theme(plot.title = element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    scale_color_viridis_d()+
    # Labels
    labs(x="Timepoint",
         y=paste0("Scores PC1 (",
                  expl_var_ab[1] %>%
                    round(1)%>%
                    format(nsmall=1),
                  "%)"),
         color="Treatment Group"),
  # LOADINGS PC1
  tibble(var = LETTERS[1:nrow(Pab)],
         loading = Pab[,1])%>%
    # Initialize plot
    ggplot(aes(x=var,
               y=loading))+
    # Geometric objects
    geom_bar(stat="identity",
             fill = c(rep(teal,4),
                      rep(dark_purple,8),
                      rep(teal,4)))+
    geom_point()+
    # Theme
    theme_classic(base_family="Times")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Labels
    # Labels
    labs(x = "Variable",
         y = paste0("Loadings PC1 (",
                    expl_var_ab[1] %>%
                      round(1)%>%
                      format(nsmall=1),
                    "%)")),
  # SCORES PC2
  design_matrix %>% 
    # Modify design matrix
    mutate(score = Tab[,2])%>%
    # Initialize plot
    ggplot(aes(x=timepoint,
               y=score,
               color=treatment))+
    # Geometric objects
    geom_point(size=5)+
    geom_line(aes(group=treatment),
              size=1)+
    # Theme
    theme_classic(base_family="Times")+
    theme(plot.title = element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    scale_color_viridis_d()+
    # Labels
    labs(x="Timepoint",
         y=paste0("Scores PC2 (",
                  expl_var_ab[2] %>%
                    round(1)%>%
                    format(nsmall=1),
                  "%)"),
         color="Treatment Group"),
  # LOADINGS PC2
  tibble(var = LETTERS[1:nrow(Pab)],
         loading = Pab[,2])%>%
    # Initialize plot
    ggplot(aes(x=var,
               y=loading))+
    # Geometric objects
    geom_bar(stat="identity",
             fill = c(rep(teal,4),
                      rep(dark_purple,8),
                      rep(teal,4)))+
    geom_point()+
    # Theme
    theme_classic(base_family="Times")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Labels
    labs(x="Timepoint",
         y=paste0("Loadings PC2 (",
                  expl_var_ab[2] %>%
                    round(1)%>%
                    format(nsmall=1),
                  "%)")),
  # Plot settings
  ncol=2,
  nrow=2,
  align="hv",
  common.legend = TRUE) %>%
  # Add title
  annotate_figure(top=text_grob("Interaction Effect (AB)",
                                size=30,
                                family = "Times")) %>%
  # Export to pdf
  ggexport(filename=paste0(output_dir, "Interaction_Simulated.pdf"),
           plot=.,
           width = 8,
           height= 8)
```

## Compose Matrices and unbalance dataset
Next, we will create balanced simulation matrices from the design info and the fixed effect components.
```{r Compose Matrices and Unbalance}
# Compose (Effect) Matrices
balanced_simulation_matrices <- compose_matrices(design_info,
                                                 function_seed=1,
                                                 mu=5,
                                                 effects = c("time","interaction"),
                                                 scores_time = Ta,
                                                 scores_interaction = Tab,
                                                 loadings_time = Pa,
                                                 loadings_interaction = Pab,
                                                 include_subject = TRUE,
                                                 subject_mu = 0,
                                                 subject_Sigma = .05)

# Sample from distribution with certain mean
balanced_response_matrix <- 
  generate_response_matrix(balanced_simulation_matrices$logYhat,
                           type=design_info$family,
                           prob_str_0 = design_info$prob_zero,
                           sigma = design_info$sigma) 

# Generate unbalance in dataset and merge into one list
unbalanced_simulation_matrices <- 
  gen_sim_sub_matrices(design_info,
                       design_matrix,
                       balanced_response_matrix,
                       function_seed=1,
                       balanced_simulation_matrices,
                       fixed_model_parameter_matrix=X)

# Modify Design list
design_info$X <- unbalanced_simulation_matrices$X
design_info$design <- unbalanced_simulation_matrices$design
```

## Plot Zero-inflation per zOTU
```{r}
gen_long_response_df(unbalanced_simulation_matrices) %>%
  mutate(var=zOTU)%>%
  group_by(zOTU)%>%
  group_nest()%>%
  mutate(plot = map(data,
                    function(x) x %>%
                      ggplot(aes(x=count))+
                      geom_histogram(bins=15)+
                      theme_classic(base_family="Times",
                                    base_size = 14)+
                      theme(axis.text.x = element_text(angle=45,
                                                       hjust=1))+
                      labs(x=element_blank(),
                           y=element_blank()))) %>%
  pull("plot") %>%
  ggarrange(plotlist=.,
            ncol=4,
            nrow=4,
            labels="AUTO",
            label.x=.75)%>%
  annotate_figure(bottom = text_grob("zOTU Count",family = "Times",size=20),
                  left = text_grob("Frequency (counts)",
                                   size=20,
                                   family="Times",
                                   rot=90)) %>%
  ggexport(filename=paste0(
    output_dir,
    "Zero-inflation_per_zOTU.pdf"),
    plot=.,
    width = 6,
    height= 6)
```



# Perform ZIGLMM fit on complete dataset
```{r GLMM ASCA with Jackknife}
# Fit GLMMs and Extract coefficients
glmm_fit_results <- gen_long_response_df(unbalanced_simulation_matrices) %>%
  fit_glmm(long_response_df = .,
           dist_family = design_info$family,
           cores=8) %>%
  extr_glmm_coef(model_df = .,
                 family = design_info$family)

# Perform 7-fold Jackknife
glmm_jackknife_models <- fit_jackknife_models(glmm_fit_results,
                                              design_info)
glmm_jackknife_results <- process_jackknife_models(glmm_jackknife_models,
                                                   design_info)


# Modify Loadings into factors for later plotting
glmm_jackknife_results$Pa$var  <- factor(glmm_jackknife_results$Pa$var,
                                         levels=seq(1:16))
glmm_jackknife_results$Pa      <- arrange(glmm_jackknife_results$Pa,
                                          var,
                                          PC)
glmm_jackknife_results$Pab$var <- factor(glmm_jackknife_results$Pab$var,
                                         levels=seq(1:16))
glmm_jackknife_results$Pab     <- arrange(glmm_jackknife_results$Pab,
                                          var,
                                          PC)
```



# Plot Jackknife results
## Plot Time Effect
```{r}
# Time EFFECT 
# SCORES PC1
ggarrange(
  glmm_jackknife_results$Ta %>% 
    filter(PC==1) %>%
    group_by(timepoint,treatment)%>%
    mutate(quantile_025 = quantile(score,.025),
           quantile_500 = quantile(score,.500),
           quantile_975 = quantile(score,.975)) %>%
    # Initialize Plot
    ggplot(aes(x=timepoint,
               y=quantile_500))+
    # Add Geometric Objects
    geom_point(size=5)+
    geom_errorbar(aes(ymin=quantile_025,
                      ymax=quantile_975))+
    # Modify Theme
    theme_classic(base_family="Times")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Add Labels
    labs(x="Timepoint",
         y=paste0("Scores PC1 (",
                  filter(glmm_jackknife_results$perc_a,
                         PC==1)$perc_expl %>%
                    quantile(.025) %>% 
                    round(1) %>% 
                    format(nsmall=1),
                  "-",
                  filter(glmm_jackknife_results$perc_a,
                         PC==1)$perc_expl %>%
                    quantile(.975) %>% 
                    round(1) %>% 
                    format(nsmall=1),
                  "%)")),
  # LOADINGS PC1
  glmm_jackknife_results$Pa %>% 
    filter(PC==1)%>%
    group_by(var)%>%
    mutate(quantile_025 = quantile(loading,.025),
           quantile_500 = quantile(loading,.500),
           quantile_975 = quantile(loading,.975)) %>%
    select(-loading)%>%
    distinct()%>%
    # Initialize Plot
    ggplot(aes(x=var,
               y=quantile_500))+
    # Add Geometric Objects
    geom_bar(stat="identity",
             fill = c(rep(dark_purple,4),
                      rep(teal,4),
                      rep(dark_purple,4),
                      rep(teal,4)))+
    geom_point()+
    geom_errorbar(aes(ymin=quantile_025,
                      ymax=quantile_975))+
    # Modify Theme
    theme_classic(base_family="Times")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Add Labels
    labs(x = "Variable",
         y = paste0("Loadings PC1 (",
                    filter(glmm_jackknife_results$perc_a,
                           PC==1)$perc_expl %>%
                      quantile(.025) %>% 
                      round(1) %>% 
                      format(nsmall=1),
                    "-",
                    filter(glmm_jackknife_results$perc_a,
                           PC==1)$perc_expl %>%
                      quantile(.975) %>%  
                      round(1) %>% 
                      format(nsmall=1),
                    "%)")),
  # SCORES PC2
  glmm_jackknife_results$Ta %>% 
    filter(PC==2) %>%
    group_by(timepoint,treatment)%>%
    mutate(quantile_025 = quantile(score,.025),
           quantile_500 = quantile(score,.500),
           quantile_975 = quantile(score,.975)) %>%
    # Initialize Plot
    ggplot(aes(x=timepoint,
               y=quantile_500))+
    # Add Geometric Objects
    geom_point(size=5)+
    geom_errorbar(aes(ymin=quantile_025,
                      ymax=quantile_975))+
    # Modify Theme
    theme_classic(base_family="Times")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Add Labels
    labs(x="Timepoint",
         y=paste0("Scores PC2 (",
                  filter(glmm_jackknife_results$perc_a,
                         PC==2)$perc_expl %>%
                    quantile(.025) %>% 
                    round(1) %>% 
                    format(nsmall=1),
                  "-",
                  filter(glmm_jackknife_results$perc_a,
                         PC==2)$perc_expl %>%
                    quantile(.975)%>% 
                    round(1) %>% 
                    format(nsmall=1),
                  "%)")),
  # LOADINGS PC2
  glmm_jackknife_results$Pa %>% 
    filter(PC==2)%>%
    group_by(var)%>%
    mutate(quantile_025 = quantile(loading,.025),
           quantile_500 = quantile(loading,.500),
           quantile_975 = quantile(loading,.975)) %>%
    select(-loading)%>%
    distinct()%>%
    # Initialize Plot
    ggplot(aes(x=var,
               y=quantile_500))+
    # Add Geometric Objects
    geom_bar(stat="identity",
             fill = c(rep(dark_purple,4),
                      rep(teal,4),
                      rep(dark_purple,4),
                      rep(teal,4)))+
    geom_point()+
    geom_errorbar(aes(ymin=quantile_025,
                      ymax=quantile_975))+
    # Modify Theme
    theme_classic(base_family="Times")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    # Add Labels
    labs(x = "Variable",
         y=paste0("Loadings PC2 (",
                  filter(glmm_jackknife_results$perc_a,
                         PC==2)$perc_expl %>%
                    quantile(.025) %>% 
                    round(1) %>% 
                    format(nsmall=1),
                  "-",
                  filter(glmm_jackknife_results$perc_a,
                         PC==2)$perc_expl %>%
                    quantile(.975) %>% 
                    round(1) %>% 
                    format(nsmall=1),
                  "%)")),
  # Plot settings
  ncol=2,
  nrow=2,
  align="hv") %>%
  annotate_figure(top=text_grob("Time Effect (A)",
                                family = "Times",
                                size=30)) %>%
  ggexport(filename=paste0(output_dir,
                           "Time_Estimated.pdf"),
           plot=.,
           width = 8,
           height= 8)

# Time-Treatment Interaction
ggarrange(
  # SCORES PC1
  glmm_jackknife_results$Tab %>% 
    filter(PC==1) %>%
    group_by(timepoint,treatment)%>%
    mutate(quantile_025 = quantile(score,.025),
           quantile_500 = quantile(score,.500),
           quantile_975 = quantile(score,.975)) %>%
    ggplot(aes(x=timepoint,
               y=-quantile_500,
               color=treatment))+
    geom_point(size=5)+
    geom_line(aes(group=treatment))+
    geom_errorbar(aes(ymin=-quantile_975,
                      ymax=-quantile_025))+
    theme_classic(base_family="Times")+
    labs(x="Timepoint",
         y=paste0("Scores PC1 (",
                  filter(glmm_jackknife_results$perc_ab,PC==1)$perc_expl %>%                     quantile(.025) %>%  round(1) %>% format(nsmall=1),
                  "-",
                  filter(glmm_jackknife_results$perc_ab,PC==1)$perc_expl %>%                     quantile(.975) %>%  round(1) %>% format(nsmall=1),
                  "%)"),
         color="Treatment Group")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    scale_color_viridis_d(),
  # LOADINGS PC1
  glmm_jackknife_results$Pab %>% 
    filter(PC==1)%>%
    group_by(var)%>%
    mutate(quantile_025 = quantile(loading,.025),
           quantile_500 = quantile(loading,.500),
           quantile_975 = quantile(loading,.975)) %>%
    select(-loading)%>%
    distinct()%>%
    ggplot(aes(x=var,
               y=-quantile_500))+
    geom_bar(stat="identity",
             fill = c(rep(teal,4),rep(dark_purple,8),rep(teal,4)))+
    geom_point()+
    geom_errorbar(aes(ymin=-quantile_975,
                      ymax=-quantile_025))+
    theme_classic(base_family="Times")+
    labs(x = "Variable",
         y=paste0("Loadings PC1 (",
                  filter(glmm_jackknife_results$perc_ab,PC==1)$perc_expl %>%                     quantile(.025) %>%  
                    round(1) %>% 
                    format(nsmall=1),
                  "-",
                  filter(glmm_jackknife_results$perc_ab,PC==1)$perc_expl %>%                     quantile(.975) %>%  
                    round(1) %>% 
                    format(nsmall=1),
                  "%)"))+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1),
  # SCORES PC2
  glmm_jackknife_results$Tab %>% 
    
    filter(PC==2) %>%
    group_by(timepoint,treatment)%>%
    mutate(quantile_025 = -quantile(score,.025),
           quantile_500 = -quantile(score,.500),
           quantile_975 = -quantile(score,.975)) %>%
    ggplot(aes(x=timepoint,
               y=quantile_500,
               color=treatment))+
    geom_point(size=5)+
    geom_line(aes(group=treatment))+
    geom_errorbar(aes(ymin=quantile_025,
                      ymax=quantile_975))+
    theme_classic(base_family="Times")+
    labs(x="Timepoint",
         y=paste0("Scores PC2 (",
                  filter(glmm_jackknife_results$perc_ab,PC==2)$perc_expl %>%                     quantile(.025) %>%  round(1) %>% format(nsmall=1),
                  "-",
                  filter(glmm_jackknife_results$perc_ab,PC==2)$perc_expl %>%                     quantile(.975) %>%  round(1) %>% format(nsmall=1),
                  "%)"),
         color = "Treatment Group")+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1)+
    scale_color_viridis_d(),
  # LOADINGS PC2
  glmm_jackknife_results$Pab %>% filter(PC==2)%>%
    group_by(var)%>%
    mutate(quantile_025 = -quantile(loading,.025),
           quantile_500 = -quantile(loading,.500),
           quantile_975 = -quantile(loading,.975)) %>%
    select(-loading)%>%
    distinct()%>%
    ggplot(aes(x=var,
               y=quantile_500))+
    geom_bar(stat="identity",
             fill = c(rep(teal,4),rep(dark_purple,8),rep(teal,4)))+
    geom_point()+
    geom_errorbar(aes(ymin=quantile_025,
                      ymax=quantile_975))+
    theme_classic(base_family="Times")+
    labs(x = "Variable",
         y=paste0("Loadings PC2 (",
                  filter(glmm_jackknife_results$perc_ab,PC==2)$perc_expl %>%
                    quantile(.025) %>%  
                    round(1) %>% 
                    format(nsmall=1),
                  "-",
                  filter(glmm_jackknife_results$perc_ab,PC==2)$perc_expl %>%
                    quantile(.975) %>%  
                    round(1) %>% 
                    format(nsmall=1),
                  "%)"))+
    theme(plot.title=element_text(hjust=.5),
          text=element_text(size=20),
          aspect.ratio = 1),
  # Plot settings
  ncol=2,
  nrow=2,
  align="hv",
  common.legend = TRUE)%>%
  annotate_figure(top=text_grob("Interaction Effect (AB)",
                                family="Times",
                                size=30)) %>%
  ggexport(filename=paste0(
    output_dir,
    "Interaction_Estimated.pdf"
  ),
  plot=.,
  width = 8,
  height= 8)
```

End of Simulation