PCA_mixed_model.Rmd

---
title: "R Notebook"
output: html_notebook
---


```{r setup_Props_PCA, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Partial correlation analysis

Goal is to evaluate condidtional independence of electropohysiological features.

First plot partial correlations for all measured features.
```{r}
data.sc_neurons <- data.sc %>% dplyr::select(vm:fi) %>%
    na.omit
Q_neurons <- calcQ(data.sc_neurons)

Q_neurons <- as.matrix(Q_neurons)
colnames(Q_neurons) <- colnames(data.sc_neurons)
rownames(Q_neurons) <- colnames(data.sc_neurons)

order_names <- c("vm", "ir", "sag", "tau", "resf", "resmag", "rheo", "fi", "ahp", "spkmax", "spkthr", "spkhlf")

Q_neurons <- Q_neurons[match(order_names, rownames(Q_neurons)), match(order_names, colnames(Q_neurons))]

new_names <- c("Vm", "IR", "Sag", "Tm", "ResF", "ResM", "Rheo", "F-I", "AHP", "Smax", "Sthr", "SHW")

colnames(Q_neurons) <- new_names
rownames(Q_neurons) <- new_names

(Q_neurons_plot <- ggcorr(data = NULL, cor_matrix = Q_neurons, geom = "circle", min_size = 0, max_size = 5, legend.size = 7, size = 2))
```


Next plot partical correlations for intercepts of mixed model fits.
```{r}
data_summary.mm_vsris <- prep_int_slopes_PCA(data.sc_r, "property", "mm_vsris")

data_intercepts <-  spread(data_summary.mm_vsris[,1:3], property, ind_intercept)
data_slopes <-  spread(data_summary.mm_vsris[,c(1:2, 4)], property, ind_slope)


data.sc_fits <- data_intercepts %>% dplyr::select(ahp:vm) %>%
  na.omit

Q_intercepts <- calcQ(data.sc_fits)

Q_intercepts <- as.matrix(Q_intercepts)
colnames(Q_intercepts) <- colnames(data.sc_fits)
rownames(Q_intercepts) <- colnames(data.sc_fits)

Q_intercepts <- Q_intercepts[match(order_names, rownames(Q_intercepts)), match(order_names, colnames(Q_intercepts))]

colnames(Q_intercepts) <- new_names
rownames(Q_intercepts) <- new_names

(Q_intercepts_plot <- ggcorr(data = NULL, cor_matrix = Q_intercepts, geom = "circle", min_size = 0, max_size = 5, legend.size = 7, size = 2))
```




# PCA mixed model analysis

Goal is to reduce the dimensionality of the dataset and then explore how dorsoventral location and mouse identity map onto each dimension by evaluating linear models generated using the principal components.

Carry out PCA. The fi measurements are absent from a lot of cells. For a first analysis include all columns but remove cells without fi measurments. Focus analysis on ages above P32 and below P45, and on data collected by HP.
```{r}
data.pca <- dplyr::select(data.sc, vm:fi, dvlocmm, id, housing, id, mlpos, hemi, age, housing, expr, patchdir, rectime) %>%
  filter(age > 32 & age <45 & expr == "HP")

all.pca <- prcomp(drop_na(data.pca[1:12], fi),
                  retx = TRUE,
                  centre = TRUE,
                  scale = TRUE)

plot(all.pca)
summary(all.pca)

prop.var.df <- as.data.frame(summary(all.pca)$importance[2,])
colnames(prop.var.df) <- c("PropVar")
prop.var.df$components <- 1:12
(all.pca.prop.var.plot <- ggplot(prop.var.df, aes(components, PropVar)) +
  geom_bar(stat = "identity") +
  scale_x_discrete("Component") +
  ylab("Proportion of variance") +
    theme(title = element_text(size = 9), axis.text = element_text(size = 9)))

(all.pca.biplot <- pca_biplot(as.data.frame(all.pca$rotation), fontsize = 2, order_names = order_names, new_names = c("Vm", "IR", "Sag", "Tm", "Res. F", "Res. Mag.", "Rheo", "FI", "AHPmax", "Spk. max", "Spk. thr", "Spk. HW")) +
    theme(text = element_text(size = 9)) +
    xlim(-0.5, 0.55))
```


In second analysis exclude the fi column and so include more cells.
```{r}
sub.pca <- prcomp(data.pca[1:11],
                  retx = TRUE,
                  centre = TRUE,
                  scale = TRUE)

plot(sub.pca)
summary(sub.pca)

(sub.pca.biplot <- pca_biplot(as.data.frame(sub.pca$rotation)))
```

## View relationships between principal components

```{r Prepare data for fitting model to principal components}
all.pca.x <- bind_cols(as_tibble(all.pca$x), drop_na(data.pca, fi))
sub.pca.x <- bind_cols(as_tibble(sub.pca$x), data.pca)
```


```{r Plot principal components versus location}
# For alternative analysis could replace all.pca.x with sub.pca.x
out.pca.x_g1_11 <- all.pca.x %>%
  gather("component", "value", 1:11)

pc_plot <- ggplot(data = out.pca.x_g1_11, aes(x = dvlocmm, y = value)) +
  geom_point() +
  facet_wrap(~ component)

out.pca.x_g1_5 <- all.pca.x %>%
  gather("component", "value", 1:5)
pc1to5_plot <-ggplot(data = out.pca.x_g1_5, aes(x = dvlocmm, y = value, colour = housing)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~ component, ncol = 5) +
  scale_x_continuous("DV location (mm)", c(0,1,2)) +
  theme_classic() +
  PCA_theme +
  theme(legend.position="bottom")

pc_plot
pc1to5_plot
```

## Fit mixed models to principal components

Fit mixed models to all measured properties using lmer.
```{r}
# Reform data for use with dplyr.
out.pca.x_g <- all.pca.x %>%
  gather("component", "value", 1:12) %>%
  group_by(component) %>%
  nest()

out.pca.x_g <- out.pca.x_g %>%
  mutate(mixedmodel_vsris = map(data, model_vsris)) %>%
  mutate(mixedmodel_vsris_null = map(data, model_vsris_null)) %>%
  mutate(mixedmodel_vsri = map(data, model_vsri)) %>%
  mutate(mixedmodel_vsri_null = map(data, model_vsri_null))
```

## Compare fits of mixed models to one another and to population level linear model

Extract AIC for all models
```{r Extract AIC for PCA models}
out.pca.x_g <- out.pca.x_g %>%
  mutate(vsris_glance = map(mixedmodel_vsris, broom::glance)) %>%
  mutate(vsris_null_glance = map(mixedmodel_vsris_null, broom::glance)) %>%
  mutate(vsri_glance = map(mixedmodel_vsri, broom::glance)) %>%
  mutate(AIC_vsris = map_dbl(vsris_glance, ~.$AIC)) %>%
  mutate(AIC_vsris_null = map_dbl(vsris_null_glance, ~.$AIC)) %>%
  mutate(AIC_vsri = map_dbl(vsri_glance, ~.$AIC))
```


Test whether effects of animal id are significant.
```{r Compare mixed with linear PCA models using chisq}
## linearmodel_to_fit fits: lm(value ~ dvlocmm, data = df, na.action = na.exclude)
out.pca.x_g <- out.pca.x_g %>%
  mutate(linearmodel = map(data, linearmodel_to_fit))

out.pca.x_g <- bind_cols(out.pca.x_g, mixed_vs_linear_pchisqu(out.pca.x_g, mixedmodel_vsris, linearmodel))

out.pca.x_g$mixedmodel_vsris_vslinear_pdiff_adj <- p.adjust(out.pca.x_g$mixedmodel_vsris_vslinear_pdiff, method = "BH")

table_mixedvslinear(out.pca.x_g, "mixedmodel_vsris_vslinear", "component")
```


## Extract other summary data from the model fits

Focus on the model with random intercept and slope (mixedmodel_vsris).

Store model gradient (extracted with summary / glance), marginal and conditional R2 (extracted with r.squaredGLMM) and p-value vs null model (calculated with ANOVA vs null model). Also extract model slopes.
```{r Extract model properties for PCA, warning=FALSE}
out.pca.x_g <- mixedmod_extract(out.pca.x_g, mixedmodel_vsris)

out.pca.x_g <- out.pca.x_g %>%
  mutate(anova = map2(mixedmodel_vsris, mixedmodel_vsris_null, ~anova(.x,.y))) %>%
  mutate(tidy_anova = map(anova, broom::tidy)) %>% 
  mutate(anova_p_val = map_dbl(tidy_anova, ~.$p.value[2])) %>%
  mutate(anova_p_val_adj = p.adjust(anova_p_val, method = "BH"))
```


## Generate summary table for PCA

Show model fitting results as a table.
```{r Make and save table with PCA mixed model properties}
props_for_table_PCA <- c("component", "mixedmodel_vsris_gradient_slopes", "anova_p_val_adj", "mixedmodel_vsris_marginal.r2", "mixedmodel_vsris_conditional.r2", "mixedmodel_vsris_slope_min", "mixedmodel_vsris_slope_max", "mixedmodel_vsris_vslinear_pdiff_adj")
props_table_PCA <- as.tibble(out.pca.x_g[props_for_table_PCA])
props_table_PCA$anova_p_val_adj <- format(props_table_PCA$anova_p_val_adj, digits = 3)
props_table_PCA$mixedmodel_vsris_vslinear_pdiff_adj <- format(props_table_PCA$mixedmodel_vsris_vslinear_pdiff_adj, digits = 3)
(props_table_unnest_PCA <- unnest(props_table_PCA) %>%
  knitr::kable(
  digits = 5,
  col.names = c("Feature", "Slope", "p (slope)", "Marginal R2", "Conditional R2", "Slope (min)", "Slope (max)", "p (vs linear)")
) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover")))

table_path <- paste0(getwd(), "/Tables/props_PCA")
if (stajpeg == 1) {
  save_kable_jpg_html(props_table_unnest_PCA, table_path)
} else  {props_table_unnest_PCA %>% kableExtra::save_kable(file = paste0(table_path, ".html"), self_contained = T)}
```


## Plot fits of mixed models of PCA data

As above, we want to plot for each model the prediction at location = 0 for each animal (I), the model prediction for location = 1 mm (I + S) and a line indicating the slope with start centred at the value of the population level model at location = 0.


Extract model predictions ready to generate plots.
```{r Format PCA mixed model fits ready for plotting, warning=FALSE}
combined_intercepts_slopes_PCA <- prep_int_slopes(out.pca.x_g, "component", "mixedmodel_vsris")

id_housing_PCA <-  distinct(all.pca.x, id, housing)

combined_intercepts_slopes_PCA <- left_join(combined_intercepts_slopes_PCA, id_housing_PCA, by = "id")

combined_intercepts_slopes_PCA$component_factors <- as.factor(combined_intercepts_slopes_PCA$component)

combined_intercepts_slopes_PCA$component_factors = factor(combined_intercepts_slopes_PCA$component_factors, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "PC11"))

```


Now generate the plot.
```{r Make facetted plot of PCA model fits}
IS_figure_PCA_1_11 <- ggplot(combined_intercepts_slopes_PCA, aes(x = measure, y = value_1, colour = housing)) +
  geom_line(aes(group = id)) +
  geom_jitter(aes(y = value_2), width = 0.2) +
  scale_x_discrete(limits = c("ind_intercept", "ind_intercept_slope", "global_intercept", "global_intercept_slope"), label = c("I", "I + S", "", "")) +
  facet_wrap(~component_factors) +
  theme_classic() +
  hist_theme +
  theme(axis.line.x = element_blank(), axis.ticks.x = element_blank())

IS_figure_PCA_1_5 <- ggplot(subset(combined_intercepts_slopes_PCA, component %in% c("PC1", "PC2", "PC3", "PC4", "PC5")), aes(x = measure, y = value_1, colour = housing)) +
  geom_line(aes(group = id)) +
  geom_jitter(aes(y = value_2), width = 0.2) +
  scale_x_discrete(limits = c("ind_intercept", "ind_intercept_slope", "global_intercept", "global_intercept_slope"), label = c("I", "I + S", "", "")) +
  facet_wrap(~component_factors, ncol = 5) +
  theme_classic() +
  hist_theme +
  theme(axis.line.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_text(angle = 90)) +
  theme(legend.position="none")

IS_figure_PCA_1_11
IS_figure_PCA_1_5
```

## Make figure with conditional independence and PCA analyses.
```{r}
ggdraw() +
  draw_plot(IS_figure_PCA_1_5, x = 0.57, y = 0.05, width = 0.4, height = 0.4) +
  draw_plot(pc1to5_plot, x = 0.57, y = 0.45, width = 0.4, height = 0.5) +
  draw_plot(Q_neurons_plot, x = 0, y = 0.5, width = 0.35, height = 0.45) +
  draw_plot(Q_intercepts_plot, x = 0, y = 0, width = 0.35, height = 0.45) +
  draw_plot(all.pca.prop.var.plot, x = 0.35, y = 0.55, width = 0.2, height = 0.35) +
  draw_plot(all.pca.biplot, x = 0.35, y = 0, width = 0.2, height = 0.4) +
  draw_plot_label(label = c("A", "B", "C", "D", "E", "F"), size = 20,
                  x = c(0, 0, 0.35, 0.35, 0.55, 0.55), y = c(1, 0.5, 1, 0.5, 1, 0.5))
  
```

Save the figure.
```{r Save plot of PCA model fits}
ggsave("Figures/CondInd_PCA_figure.png", width = 220, height = 120, units = "mm")
```



## Evaluate influence of individual fixed effects on PCA components

Housing
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(mm_vsris_housing = map(data, model_vsris_housing)) %>%
  summary_2fixedeffects(mm_vsris_housing, "housing") %>%
  FE2_p.adjust("housing")

(housing_table_PCA <- FE_table(out.pca.x_g, "housing", "component"))

table_path <- paste0(getwd(), "/Tables/PCA_FE_housing")

if (stajpeg == 1) {
  save_kable_jpg_html(housing_table_PCA, table_path)
} else  {housing_table_PCA %>% kableExtra::save_kable(file = paste0(table_path, ".html"), self_contained = T)}
```

Evaluate housing in datset containing only male mice.

Prepare data
```{r}
data.pca_male <- dplyr::select(data.sc, vm:fi, dvlocmm, id, housing, id, mlpos, hemi, age, sex, housing, expr, patchdir, rectime) %>%
  filter(age > 32 & age <45 & expr == "HP" & sex == "M")

all_male.pca <- prcomp(drop_na(data.pca_male[1:12], fi),
                  retx = TRUE,
                  centre = TRUE,
                  scale = TRUE)

plot(all_male.pca)
summary(all_male.pca)

all.pca_male.x <- bind_cols(as_tibble(all_male.pca$x), drop_na(data.pca_male, fi))

out.pca.x_male_g <- all.pca_male.x %>%
  gather("component", "value", 1:12) %>%
  group_by(component) %>%
  nest()
```

Do analysis
```{r}
out.pca.x_male_g <- out.pca.x_male_g %>%
  mutate(mm_vsris_housing = map(data, model_vsris_housing)) %>%
  summary_2fixedeffects(mm_vsris_housing, "housing") %>%
  FE2_p.adjust("housing")

(housing_male_table_PCA <- FE_table(out.pca.x_male_g, "housing", "component"))

table_path <- paste0(getwd(), "/Tables/PCA_FE_housing_male")

if (stajpeg == 1) {
  save_kable_jpg_html(housing_male_table_PCA, table_path)
} else  {housing_male_table_PCA %>% kableExtra::save_kable(file = paste0(table_path, ".html"), self_contained = T)}
```


Evaluate whether interanimal differences remain in models that account for housing.

```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(lm_vsris_housing = map(data, linearmodel_housing))
out.pca.x_g <- mixed_vs_linear_pchisqu(out.pca.x_g,
                                         mm_vsris_housing,
                                         lm_vsris_housing)

out.pca.x_g$mm_vsris_housing_vslinear_pdiff_adj <- p.adjust(out.pca.x_g$mm_vsris_housing_vslinear_pdiff, method = "BH")
(PCA_mix_vs_linear_table <- table_mixedvslinear(out.pca.x_g, "mm_vsris_housing_vslinear", "component"))

if (stajpeg == 1) {
  table_path <- paste0(getwd(), "/Tables/PCA_mm_vsris_housing_vslinear.jpg")
   PCA_mix_vs_linear_table %>% kableExtra::save_kable(file = table_path, self_contained = T)
} 

table_path <- paste0(getwd(), "/Tables/PCA_mm_vsris_housing_vslinear.html")
PCA_mix_vs_linear_table %>% kableExtra::save_kable(file = table_path, self_contained = T)
```

Sex
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(mm_vsris_sex = map(data, model_vsris_sex)) %>%
  summary_2fixedeffects(mm_vsris_sex, "sex") %>%
  FE2_p.adjust("sex")

(sex_table_PCA <- FE_table(out.pca.x_g, "sex", "component"))

table_path <- paste0(getwd(), "/Tables/PCA_FE_sex.html")
sex_table_PCA %>% kableExtra::save_kable(file = table_path, self_contained = T)
```



Mediolateral position
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(mm_vsris_mlpos = map(data, model_vsris_ml)) %>%
  summary_2fixedeffects(mm_vsris_mlpos, "mlpos") %>%
  FE2_p.adjust("mlpos")

(ml_table_PCA <- FE_table(out.pca.x_g, "mlpos", "component"))

table_path <- paste0(getwd(), "/Tables/PCA_FE_mlpos.html")
ml_table_PCA %>% kableExtra::save_kable(file = table_path, self_contained = T)
```

Hemisphere.
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(mm_vsris_hemi = map(data, model_vsris_hemi)) %>%
  summary_2fixedeffects(mm_vsris_hemi, "hemi") %>%
  FE2_p.adjust("hemi")

(hemi_table_PCA <- FE_table(out.pca.x_g, "hemi", "component"))

table_path <- paste0(getwd(), "/Tables/PCA_FE_hemi.html")
hemi_table_PCA %>% kableExtra::save_kable(file = table_path, self_contained = T)
```

Direction.
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(mm_vsris_dir = map(data, model_vsris_dir)) %>%
  summary_2fixedeffects(mm_vsris_dir, "dir") %>%
  FE2_p.adjust("dir")

(dir_table_PCA <- FE_table(out.pca.x_g, "dir", "component"))

table_path <- paste0(getwd(), "/Tables/PCA_FE_dir.html")
dir_table_PCA %>% kableExtra::save_kable(file = table_path, self_contained = T)
```


Recording time.
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(mm_vsris_rect = map(data, model_vsris_rect)) %>%
  summary_2fixedeffects(mm_vsris_rect, "rect") %>%
  FE2_p.adjust("rect")

(rect_table_PCA <- FE_table(out.pca.x_g, "rect", "component"))

table_path <- paste0(getwd(), "/Tables/PCA_FE_rect.html")
rect_table_PCA %>% kableExtra::save_kable(file = table_path, self_contained = T)
```


## Evaluate PCA using intercepts and slopes from models fits as inputs.

Obtain intercepts and slopes from previously fit models.
```{r, warning=FALSE}
data_summary.pca <- prep_int_slopes_PCA(data.sc_r, "property", "mm_vsris")

data_intercepts_PCA <-  spread(data_summary.pca[,1:3], property, ind_intercept)
data_slopes_PCA <-  spread(data_summary.pca[,c(1:2, 4)], property, ind_slope)
```

Carry out PCA for intercepts
```{r}
intercepts.pca <- prcomp(data_intercepts_PCA[2:13],
                  retx = TRUE,
                  centre = TRUE,
                  scale = TRUE)

plot(intercepts.pca)
summary(intercepts.pca)

(intercepts.pca.biplot <- pca_biplot(as.data.frame(intercepts.pca$rotation)))
```

Carry out PCA for slopes.
```{r}
slopes.pca <- prcomp(data_slopes_PCA[2:13],
                  retx = TRUE,
                  centre = TRUE,
                  scale = TRUE)

plot(slopes.pca)
summary(slopes.pca)

(slopes.pca.biplot <- pca_biplot(as.data.frame(slopes.pca$rotation)))
```

Examine distribution of points
```{r}
id_lookup <- unique(data.sc[c("id", "age", "housing")])

slopes.pca_x <- bind_cols(as_tibble(slopes.pca$x), as.tibble(data_intercepts_PCA$id))

slopes.pca_x$age <- id_lookup$age[match(slopes.pca_x$value, id_lookup$id)]
slopes.pca_x$housing <- id_lookup$housing[match(slopes.pca_x$value, id_lookup$id)]

(pc_plot_12 <- ggplot(data = slopes.pca_x, aes(x = PC1, y = PC2)) +
  geom_point(aes(colour = housing)))
```