--- 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) ```