PCA.Rmd

---
title: "PCA"
author: "Matt Nolan"
date: "02/05/2018"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

# Ensure access to libraries
library(tidyverse)
library (pls)
library(GGally)
library(Rtsne)
```

## Goals

To reduce the dimensionality of the dataset. To explore how dorsoventral location and mouse identity map onto each dimension. To evaluate linear models generated using the principal components.

Carry out PCA.
```{r}
data.pca <- dplyr::select(data.sc, vm:fi, dvlocmm, id, housing)
cols.pca <- 1:11

out.pca <- prcomp(data.pca[cols.pca],
                  retx = TRUE,
                  centre = TRUE,
                  scale = TRUE)

plot(out.pca)
summary(out.pca)
biplot(out.pca)
```



Plot components vs dosrosventral location. Colour code mouse ID.
```{r}
out.pca.x <- as_tibble(out.pca$x)
out.pca.x$dvlocmm <- data.pca$dvlocmm
out.pca.x$id <- data.pca$id
out.pca.x$housing <- data.pca$housing

out.pca.x_g1_11 <- out.pca.x %>%
  gather("component", "value", 1:11)

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

out.pca.x_g1_5 <- out.pca.x %>%
  gather("component", "value", 1:5)
pc1to5_plot <-ggplot(data = filter(out.pca.x_g1_5), aes(x = dvlocmm, y = value)) +
  geom_point() +
  facet_wrap(~ component, ncol = 5) +
  theme_classic() +
  hist_theme

pc_plot
pc1to5_plot
```


Plot components seperately for each mouse.
```{r}
ggplot(data = out.pca.x, aes(x = id, y = value, colour = housing)) +
  geom_boxplot() +
  coord_flip() +
  facet_wrap(~ component, scales = "free_x")
```

Plot components against one another.
```{r}
out.pca.x_2 <- as_tibble(out.pca$x)
out.pca.x_2$dvlocmm <- data.pca$dvlocmm
out.pca.x_2$id <- data.pca$id
pairs1 <- ggpairs(out.pca.x_2,columns = c(1:4), ggplot2::aes(colour=id, alpha = 0.1))

pairs2 <- ggpairs(out.pca.x_2,columns = c(1:2), ggplot2::aes(colour=id, alpha = 0.1))

pairs1
pairs2
```


Reform data for use with dplyr.
```{r}
out.pca.x_g <- group_by(out.pca.x, component) %>%
  nest()
```

Fit mixed model for each principal component as a function of dorsoventral position.
```{r Fit mixed model to PCA}
out.pca.x_g <- out.pca.x_g %>%
  mutate(mixedmodel = map(data, model_to_fit))
```

Fit mixed models to all measured properties using lmer. Don't fit vscri as it fails to converge.
```{r}
# model_to_fit moved to functions.
# lme4::lmer(value ~ dvlocmm +(1+dvlocmm||id), data = df, REML = FALSE, na.action = na.exclude)
# out.pca.x_g <- out.pca.x %>%
#   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))
```



Add model predictions. Using predict adds a prediction for each location inclusing random effects or at the population level (all random effects set to zero.) Using broom::augment adds predictions (.mu) and residuals (.wtres). Package broom.mixed might be useful in future (https://github.com/bbolker/broom.mixed).
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(fit_random = map(mixedmodel, predict))

out.pca.x_g <- out.pca.x_g %>%
  mutate(fit_population = map(mixedmodel, re.form = NA, predict))

out.pca.x_g <- out.pca.x_g %>%
  mutate(fits = map(mixedmodel, broom::augment))
```


Add summaries of model output
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(tidy = map(mixedmodel, broom::tidy))
```

Extract coefficients for individual mice. The function is to convert the output of coef() to a dataframe so it will work with unnest.
```{r}
out.pca.x_g <- out.pca.x_g %>%
  mutate(coefs = map(mixedmodel, coef_df))
```


Convert data from the nested format so that we can make plots of the fitted data.
```{r}
out.pca.x_g_fit <- out.pca.x_g %>%
  unnest(data, fit_random, fit_population)
```

Make plots
```{r}
ggplot(out.pca.x_g_fit, aes(x = dvlocmm, y = value, colour = housing)) +
  geom_point(alpha = 0.05) +
  geom_line(aes(y=fit_random, group = id), size=0.8, alpha = 0.5) +
  geom_line(aes(y=fit_population), colour = "black") +
  facet_wrap(~component, scales = "free")
```



Extract individual slopes and intercepts from model fits.
```{r}
out.pca.x_g_coefs <- out.pca.x_g %>%
  unnest(coefs)
```


Add slope to intercept to predict values at 1 mm for each mouse.
```{r}
out.pca.x_g_coefs$is <- out.pca.x_g_coefs$intercept + out.pca.x_g_coefs$slope
```

To enable plotting of slopes on the same graph, but seperately from the intercepts, make columns containing population intercepts, and population intercept + slope for each mouse.
```{r}
pop_intercepts_PCA <- unnest(out.pca.x_g, tidy) %>% filter(term == "(Intercept)") %>% dplyr::select(component, estimate)
out.pca.x_g_coefs_regather <- out.pca.x_g_coefs %>%
  dplyr::select(component, id, slope) %>%
  spread(key = id, value = slope) %>%
  left_join(pop_intercepts_PCA, by = "component") %>%
  gather("id", "slope", -estimate, -component) %>%
  mutate(intercept_slope = estimate + slope) %>%
  dplyr::select(component, id, estimate, intercept_slope) %>%
  gather(measure, value, estimate, intercept_slope)
```



Plot predicted values at 0 mm and at 1mm for each mouse.
```{r}
out.pca.x_g_coefs_01 <- out.pca.x_g_coefs %>%
  dplyr::select(id, component, intercept, is) %>%
  gather(measure, value, intercept, is)

intercept_slope_plot_a <- ggplot(out.pca.x_g_coefs_01, aes(x = measure, y = value)) +
  geom_line(aes(group = id)) +
  facet_wrap(~component, scales = "free") +
  theme_classic() +
  hist_theme

intercept_slope_plot_b <- ggplot(out.pca.x_g_coefs_01, aes(x = measure, y = value)) +
  geom_jitter(width = 0.2, height = 0) +
  facet_wrap(~component, scales = "free") +
  theme_classic() +
  hist_theme
intercept_slope_plot_a
intercept_slope_plot_b
```

Combine plots of intercept and slopes. Need to combine out.pca.x_g_coefs_01 and out.pca.x_g_coefs_regather and add housing. Make property into factors to enable ordering of plots. Then make appropriately faceted plots.

```{r}
out.pca.x_g_coefs_regather <- mutate(out.pca.x_g_coefs_regather, value_1 = value)
out.pca.x_g_coefs_01 <- mutate(out.pca.x_g_coefs_01, value_2 = value)
combined_intercepts_slopes_PCA <- bind_rows(out.pca.x_g_coefs_regather, out.pca.x_g_coefs_01)

id_housing <-  distinct(out.pca.x, id, housing)
combined_intercepts_slopes_PCA <- left_join(combined_intercepts_slopes_PCA, id_housing, 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"))

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("intercept", "is",  "estimate", "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_11
```

---------------------- Compare mixed with linear models -------------

To test whether effects of animal id are significant compare mixed model fits with linear model fits. Modified from: https://web.stanford.edu/class/psych252/section/Mixed_models_tutorial.html.
```{r Compare mixed with linear model 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))
```

Compare the linear and mixed models using ANOVA and mixed model fit with nlme. 
```{r}
# nlmemodel_to_fit: implements lme(value ~ dvlocmm, random = ~1|id, data = df, method = "ML", na.action = na.exclude)

ctrl <- lmeControl(opt='optim')
out.pca.x_g <- out.pca.x_g %>%
  mutate(nlmemodel = map(data, nlmemodel_to_fit))

# quickly check first 4 principal components
anova(out.pca.x_g$nlmemodel[[1]],out.pca.x_g$linearmodel[[1]])
anova(out.pca.x_g$nlmemodel[[2]],out.pca.x_g$linearmodel[[2]])
anova(out.pca.x_g$nlmemodel[[3]],out.pca.x_g$linearmodel[[3]])
anova(out.pca.x_g$nlmemodel[[4]],out.pca.x_g$linearmodel[[4]])
```

-------------- t-SNE ----------------------

```{r}
p.x <- as_tibble(out.pca$x)
p.x_unique <- unique(p.x)
p.x_mat <- as.matrix(p.x_unique[,1:11])
set.seed(1)
tSNE_out <- Rtsne(p.x_mat, perplexity = 50, max_iter = 10000, pca = FALSE)
```


Plot tSNE output.
```{r}
t.o <- as.tibble(tSNE_out$Y)
t.o$dvlocmm <- data.pca$dvlocmm
t.o$housing <- data.pca$housing
t.o$id <- data.pca$id
ggplot(t.o, aes(V1,V2, colour = dvlocmm)) +
  geom_point()
```

----------------- Model analysis ---------------------


Compare different model with null models.
```{r}
# Store model gradient (extracted with summary), marginal and conditional R2 (extracted with r.squaredGLMM) and p-value vs null model (calculated with ANOVA vs null model).
out.pca.x_g <- out.pca.x_g %>%
  mutate(vsris_summary = map(mixedmodel_vsris, summary)) %>%
  mutate(vsris_null_summary = map(mixedmodel_vsris_null, summary)) %>%
  mutate(vsri_summary = map(mixedmodel_vsri, summary)) %>%
  mutate(gradient_slopes = map_dbl(vsris_summary, ~.$coefficients[[2]])) %>%
  mutate(extractR2 = map(mixedmodel_vsris, r.squaredGLMM)) %>%
  mutate(marginal.r2 = map_dbl(extractR2, ~.[[1]])) %>%
  mutate(conditional.r2 = map_dbl(extractR2, ~.[[2]])) %>%
  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]))

# Extract model slopes
out.pca.x_g <- out.pca.x_g %>%
  mutate(vsris_simcoefs = map(mixedmodel_vsris, ~summary(coef(.x)[[1]][[2]]))) %>%
  mutate(modelslope_min = map_dbl(vsris_simcoefs, ~.[[1]])) %>%
  mutate(modelslope_median = map_dbl(vsris_simcoefs, ~.[[3]])) %>%
  mutate(modelslope_max = map_dbl(vsris_simcoefs, ~.[[5]]))
  
# Extract AIC for all models
out.pca.x_g <- out.pca.x_g %>% 
  mutate(AIC_vsris = map_dbl(vsris_summary, ~.$AIC[[1]])) %>%
  mutate(AIC_vsris_null = map_dbl(vsris_null_summary, ~.$AIC[[1]])) %>%
  mutate(AIC_vsri = map_dbl(vsri_summary, ~.$AIC[[1]]))

```

Show model fitting results as a table.
```{r}
props_for_table <- c("component", "gradient_slopes", "modelslope_min", "modelslope_max", "anova_p_val", "marginal.r2", "conditional.r2", "pdiff")
props_table <- as.tibble(out.pca.x_g[props_for_table])
props_table_unnest <- unnest(props_table)

  knitr::kable(
  props_table_unnest,
  digits = 5,
  caption = "Fit of measured membrane properties as a function of location"
)

write_csv(props_table_unnest, "PCA_results_model_table.csv")
```


-----Directly compare inter-animal variance with population variance---------

Predictions from model coefficients from each mouse are in fit_random and for the global model are in fit_population. Subtract coefficients from datavalues to obtain the residuals.
```{r}
out.pca.x_g_fit <- out.pca.x_g_fit %>%
  mutate(random_resid = fit_random - value) %>%
  mutate(population_resid = fit_population - value)
```

Plot residuals for each principal component
```{r}
out.pca.x_g_fit$component_factors <- as.factor(out.pca.x_g_fit$component)
out.pca.x_g_fit$component_factors = factor(out.pca.x_g_fit$component_factors, c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "PC11"))

ggplot(out.pca.x_g_fit, aes(dvlocmm, random_resid)) +
  geom_point() +
  facet_wrap(~component_factors)
```


Make plots of distribution of residuals for each principal component.
```{r Distribution of residuals}
ggplot(out.pca.x_g_fit, aes(random_resid)) +
  geom_density(colour = "red") +
  geom_density(aes(population_resid), colour = "black") +
  facet_wrap(~component_factors, scales = "free")
```

Calculate squared errors for each component.
```{r}
out.pca.x_g_fit <- out.pca.x_g_fit %>%
  mutate(random_square = random_resid^2) %>%
  mutate(population_square = population_resid^2)
```

Plot distributions of squared errors for each component.
```{r}
ggplot(out.pca.x_g_fit, aes(random_square)) +
  geom_density(colour = "red") +
  geom_density(aes(population_square), colour = "black") +
  xlim(0, 10) +
  facet_wrap(~component_factors, scales = "free")
```



Statistical comparison? Manual ANOVA? What's a correct approach for predictions from a mixed model?

Calculate variance as the mean of the squares for each principal component.
```{r}
grouped_pca <- group_by(out.pca.x_g_fit, component)
summarise(grouped_pca, var_PCA = mean(random_square))
summarise(grouped_pca, var_PCA = mean(population_square))
```


Calculate f from the ratio of between group variance to within group variance.

Calcluate p from f.


A slightly different approach, directly fitting model to PCR output.
```{r}
data.pca <- dplyr::select(data.sc.norm, vm:ahp, dvlocmm)
pcr.fit=pcr(dvlocmm~., data=data.pca,scale=TRUE , validation ="CV")
summary(pcr.fit)
```

Evaluate mean squared error of prediction
```{r}
validationplot(pcr.fit ,val.type="MSEP")
```