Posterior_test_comparative_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Posterior test for phylogenetic comparative analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "library(\"tidyverse\")\n",
    "library(\"ggpubr\")\n",
    "library(\"rstatix\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing arguments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "args<-c(\"Path_to_replicates\", \"pattern\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Listing files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_files<-list.files(args[1], pattern = as.character(args[2]), full.names = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing meta dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df<-data.frame()\n",
    "Model<-c()\n",
    "AICc<-c()\n",
    "Replicate<-c()\n",
    "for( i in list_files){\n",
    "        iteration<-gsub(paste(args[1], \"/df_\", as.character(args), sep =\"\"),'', i)\n",
    "        iteration<-gsub('.tsv', '',  iteration)\n",
    "        file<-read.csv(i, sep =\"\\t\")\n",
    "        for (j in 1:nrow(file)){\n",
    "            Model<-c(Model, rownames(file[j,]))\n",
    "            AICc<-c(AICc, file[j,2])   \n",
    "            Replicate<-c(Replicate, as.numeric(iteration))\n",
    "        }\n",
    "}\n",
    "df<-as.data.frame(cbind(Model, AICc, Replicate))\n",
    "df$AICc<-as.numeric(df$AICc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shapiro test for normality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab_test_shap <- df %>%\n",
    "  group_by(Model) %>%\n",
    "  shapiro_test(AICc)\n",
    "\n",
    "tab_test_shap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AICc summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "tab_test_mean <-df %>%\n",
    "  group_by(Model) %>%\n",
    "  get_summary_stats(AICc, type = \"mean_sd\")\n",
    "tab_test_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove duplicates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_dups <- df[c(\"Model\", \"Replicate\")]\n",
    "\n",
    "df<-df[!duplicated(df_dups),]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform Friedmann test & pairwise Wilcoxon test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.aov <- df %>% friedman_test(AICc ~ Model |Replicate)\n",
    "    pwc <- df %>%\n",
    "  wilcox_test(AICc ~ Model, paired = TRUE, p.adjust.method = \"bonferroni\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "res.aov\n",
    "pwc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwc2 <- pwc %>% add_xy_position(x = \"time\")\n",
    "    ggboxplot(df, x = \"Model\", y = \"AICc\", add = \"point\") +\n",
    "      stat_pvalue_manual(pwc2, hide.ns = TRUE) +\n",
    "      labs(\n",
    "        subtitle = get_test_label(res.aov,  detailed = TRUE),\n",
    "        caption = get_pwc_label(pwc2)\n",
    "          )\n",
    "ggsave(\"result_AICc_boxplot.pdf\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.3.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}