---
title: "R Notebook"
output: html_notebook
---
Packages
```{r}
library(tidyverse)
library(cowplot)
library(RRPP)
library(MVN)
```
Data upload - will need to change data location
```{r}
pta_data <- read_csv("./data/PTA_data.csv",
col_names = TRUE,
col_types = list(pathogen = col_factor(levels = c("CTRL", "C1", "C14", "C20")),
mat_temp = col_factor(levels = c("20", "25")),
acc_temp = col_factor(levels = c("20", "25")),
temp_treat = col_factor(levels = c("20-20", "20-25", "25-20", "25-25")),
exposed = col_factor(NULL))) %>%
na.omit %>%
unite("path_acc_temp", c(pathogen, acc_temp), sep = "_", remove = FALSE) %>%
mutate(path_acc_temp = as.factor(path_acc_temp))
pta_data
```
Predictor variables:
pathogen
mat_temp
acc_temp
Phenotype variables:
size_at_mat
size_at_20
growth_rate
age_at_maturity
age_of_first_clutch
size_first_clutch
total_offspring
total_clutches
mean_clutch_size
little_r
lifespan
sample sizes
```{r}
pta_data %>%
group_by(pathogen, mat_temp, acc_temp) %>%
summarise(n = sum(!is.na(total_clutches)))
```
Transforming total_offspring, total_clutches & little_r
```{r}
pta_data2 <- pta_data %>%
mutate(log_offspring = log(total_offspring),
log_clutches = log(total_clutches),
sq_little_r = little_r^2)
pta_data2
```
Scaleing and centering variables
```{r}
scaled_pta_data <- pta_data2 %>%
mutate_at(vars(size_at_mat, size_at_20, growth_rate, age_at_maturity, age_of_first_clutch, size_first_clutch, log_offspring, log_clutches, mean_clutch_size, sq_little_r, lifespan),
funs(z = scale(., center = TRUE, scale = TRUE)))
```
Checking distributions
```{r}
#run through each predictor
hist(scaled_pta_data$log_clutches_z)
```
Selecting scaled variables and creating matrix
```{r}
Y <- scaled_pta_data %>%
select(size_at_mat_z, size_at_20_z, growth_rate_z, age_at_maturity_z, age_of_first_clutch_z, size_first_clutch_z, log_offspring_z, log_clutches_z, mean_clutch_size_z, sq_little_r_z, lifespan_z) %>%
as.matrix()
ta_temp_data <- rrpp.data.frame(Y = Y,
pathogen = scaled_pta_data$pathogen,
mat_temp = scaled_pta_data$mat_temp,
acc_temp = scaled_pta_data$acc_temp,
path_acc_temp = scaled_pta_data$path_acc_temp)
glimpse(ta_temp_data)
```
Fit multivaraite model
```{r}
fit <- lm.rrpp(Y ~ mat_temp + acc_temp * pathogen, SS.type = "III",
data = ta_temp_data,
print.progress = FALSE)
anova(fit, effect.type = "F")
```
Convert model to MANOVA
```{r}
fitm <- manova.update(fit, print.progress = FALSE, tol = 0)
summary(fitm, test = "Pillai")
```
Conduct the PTA - trajectory analysis
```{r}
TA <- trajectory.analysis(fitm, groups = ta_temp_data$pathogen, traj.pts = ta_temp_data$acc_temp,
print.progress = FALSE)
# Magnitude difference (absolute difference between path distances)
TA_summary <- summary(TA, attribute = "MD")
TA_summary$summary.table
```
```{r}
# Correlations (angles) between trajectories
TA_summary_tc <- summary(TA, attribute = "TC", angle.type = "deg")
TA_summary_tc$summary.table
```
PCA rotations
```{r}
TA$pca$rotation
```
principle components
```{r}
TA$pca$sdev %>% as_tibble() %>%
rename(SD = value) %>%
mutate(Variance = SD^2) %>%
mutate(Percent_var = round((Variance/sum(Variance))*100, 3))
```
Basic plot
```{r}
# Plot results
TP <- plot(TA, pch = as.numeric(ta_temp_data$pathogen) + 20, bg = as.numeric(ta_temp_data$acc_temp),
cex = 0.7, col = "gray")
add.trajectories(TP)
```
Collecting data & PCs for plotting
```{r}
gg_tp_pca <- TP$pc.points %>% as_tibble() %>%
bind_cols(pathogen = ta_temp_data$pathogen) %>%
bind_cols(acc_temp = ta_temp_data$acc_temp) %>%
bind_cols(path_acc_temp = ta_temp_data$path_acc_temp)
gg_tp_pca_sum <- gg_tp_pca %>% group_by(pathogen, acc_temp) %>%
summarise(.groups = "drop",
PC1_mn = mean(PC1, na.rm = TRUE),
PC1_SE = sd(PC1, na.rm = TRUE)/sqrt(sum(!is.na(PC1))),
PC2_mn = mean(PC2, na.rm = TRUE),
PC2_SE = sd(PC2, na.rm = TRUE)/sqrt(sum(!is.na(PC2))))
ggplot(gg_tp_pca_sum, aes(PC1_mn, PC2_mn, group = pathogen, colour = pathogen)) +
geom_point() +
geom_line(colour = "black") +
geom_point(gg_tp_pca, mapping = aes(PC1, PC2),
inherit.aes = FALSE) +
theme_bw() +
theme(aspect.ratio=1)
```
Phenotypic Trajectory plot
```{r}
(PCA_plot_1 <- ggplot(gg_tp_pca_sum, aes(PC1_mn, PC2_mn, group = pathogen, shape = acc_temp, colour = pathogen)) +
geom_errorbar(aes(ymin = PC2_mn - PC2_SE, ymax = PC2_mn + PC2_SE), width = 0) +
geom_errorbarh(aes(xmin = PC1_mn - PC1_SE, xmax = PC1_mn + PC1_SE)) +
geom_line() +
geom_point(size = 3) +
scale_shape_manual(values = c(15, 17), labels = c("20°C", "25°C"), name = "Temperature") +
scale_colour_manual(values = c("#008000", "#e55326","#5e3096","#e72a8a"),
labels = c("Uninfected", "C1", "C14", "C20"),
name = "Pathogen") +
xlab("PC1 (62.02% variance explained)") + ylab("PC2 (35.38% variance explained)") +
theme_bw() +
coord_fixed(ratio = 1) +
ylim(c(-3, 3)) +
xlim(c(-3, 3)))
```
Principle components plot
```{r}
test_pca <- TA$pca$rotation %>% as_tibble(rownames = "trait")
test_pca.arrows = data.frame(x1 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
y1 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
x2 = test_pca$PC1,
y2 = test_pca$PC2,
trait = test_pca$trait) %>%
mutate(trait = recode(trait,
size_at_mat_z = "Size at maturity",
size_at_20_z = "Size at day 20",
growth_rate_z = "Growth rate",
age_at_maturity_z = "Age at maturity",
age_of_first_clutch_z = "Age at first clutch",
size_first_clutch_z = "Size of first clutch",
log_offspring_z = "Lifetime fecundity",
log_clutches_z = "Number of clutches",
mean_clutch_size_z = "Mean clutch size",
sq_little_r_z = "Rate of increase (r)",
lifespan_z = "Lifespan"))
(PCA_plot_2 <- ggplot(gg_tp_pca, aes(x = PC1, y = PC2, colour = pathogen, shape = pathogen)) +
coord_fixed(ratio = 1) +
geom_point(size = 2, alpha = 0.5) +
stat_ellipse(aes(x = PC1, y = PC2, group = path_acc_temp, linetype = acc_temp), level = 0.95) +
scale_colour_manual(values = c("#008000", "#e55326","#5e3096","#e72a8a"),
labels = c("Uninfected", "C1", "C14", "C20"),
name = "Pathogen") +
scale_shape_manual(values = c(6,1,2,0),
labels = c("Uninfected", "C1", "C14", "C20"),
name = "Pathogen") +
scale_linetype_manual(values = c(1, 2),
labels = c("20°C", "25°C"),
name = "Temperature") +
geom_segment(data = test_pca.arrows, aes(x = x1, y = y1, xend = x2*6, yend = y2*6),
colour = "black", arrow = arrow(length = unit(1/2, 'picas')), inherit.aes = FALSE) +
# ylim(c(-4,4)) + xlim(c(-4,4)) +
xlab("PC1 (62.02% variance explained)") + ylab("PC2 (35.38% variance explained)") +
theme_bw() +
geom_text(data=test_pca.arrows, mapping = aes(x = x2*6.5, y = y2*6.5, label = trait),
nudge_x = c(-1.5,0,0.8,0.2,0.6,-1.5,-0.5,0.5,-0.5,0,0),
nudge_y = c(0,0.2,-0.2,-0.1,-0.6,0.2,0,0.3,0.2,0.5,0.5),
inherit.aes = FALSE))
```
Assembling final figure
```{r}
plots_aligned <- align_plots(PCA_plot_1, PCA_plot_2, align = "hv")
PCA_full_plot <- plot_grid(plots_aligned[[1]], plots_aligned[[2]], ncol = 1, labels = c("A", "B"))
#will need to change save location
ggsave2("./figures/PCA_full_plot.pdf", PCA_full_plot, width = 6, height = 8)
```