Raw_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparing morphological measurments for comparative analyses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cleaning and loading data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Le chargement a nécessité le package : RRPP\n",
      "\n",
      "Le chargement a nécessité le package : rgl\n",
      "\n",
      "This build of rgl does not include OpenGL functions.  Use\n",
      " rglwidget() to display results, e.g. via options(rgl.printRglwidget = TRUE).\n",
      "\n",
      "Le chargement a nécessité le package : Matrix\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(\"geomorph\")\n",
    "library(\"stringr\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Landmarks file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "No curves detected; all points appear to be fixed landmarks.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in .readland.tps(file, specID, negNA, readcurves, warnmsg):\n",
      "“Not all specimens have scale adjustment (perhaps because they are already scaled); \n",
      "no rescaling will be performed in these cases\n",
      "”\n"
     ]
    }
   ],
   "source": [
    "landmarks <- readland.tps(\"data/HYDRO_LM_all_withNAs.tps\",\n",
    "                          specID = \"imageID\",\n",
    "                          negNA = FALSE,\n",
    "                          readcurves = FALSE,\n",
    "                          warnmsg = TRUE\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Verterbrae count file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vert_Count <- read.csv(\"data/Vertebral_Count.csv\")\n",
    "Vert_Count$Short_Code <- sapply(strsplit(Vert_Count$Full_Code, \"\"), function(x) paste(x[1:7], collapse = \"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ecological data file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_eco<-read.table(\"data/Cat_Var.csv\", header = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loading vertebral count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "vertnumb <- read.csv(\"data/Caudal_Count_All.csv\", header = FALSE, row.names = 1, sep = \",\")\n",
    "rnames <- rownames(vertnumb)\n",
    "colnames(vertnumb) <- \"Vert_numb\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cleaning data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Remove missidentify specimen (Protochromys fellowsi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "landmarks <- landmarks[,,-74]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Substetting datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "genus<-as.factor(substr(dimnames(landmarks)[[3]], 1,3))\n",
    "species <- as.factor(substr(dimnames(landmarks)[[3]],1,7))\n",
    "locomotion<-as.factor(substr(dimnames(landmarks)[[3]], 8,10))\n",
    "climbing<-as.factor(substr(dimnames(landmarks)[[3]], 11,13))\n",
    "prehensility_TER<-as.factor(substr(dimnames(landmarks)[[3]], 14,15))\n",
    "prehensility_BIN<-as.factor(substr(dimnames(landmarks)[[3]], 16,16))\n",
    "prehensility_POST<-as.factor(substr(dimnames(landmarks)[[3]], 17,18))\n",
    "tribe<-as.factor(substr(dimnames(landmarks)[[3]], 19,21))\n",
    "coll<-as.factor(substr(dimnames(landmarks)[[3]], 22,25))\n",
    "sex<-as.factor(substr(dimnames(landmarks)[[3]], 26,27))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Attributing species names to landmark table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "dimnames(landmarks)[[3]]<-substr(dimnames(landmarks)[[3]],1,7)\n",
    "namesp<-dimnames(landmarks)[[3]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Linear extraction from landmarks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "vec<-c(\"FMRMID\")\n",
    "for(i in 1:46){\n",
    "       if(i < 10){\n",
    "            vec<-c(vec,c(paste(\"V0\", i, \"TPW\", sep = \"\"), paste(\"V0\", i, \"MDW\", sep = \"\"), paste(\"V0\", i, \"CTL\", sep = \"\"), paste(\"V0\", i, \"CTH\", sep = \"\")))\n",
    "       }\n",
    "        else{\n",
    "            vec<-c(vec,c(paste(\"V\", i, \"TPW\", sep = \"\"), paste(\"V\", i, \"MDW\", sep = \"\"), paste(\"V\", i, \"CTL\", sep = \"\"), paste(\"V\", i, \"CTH\", sep = \"\")))\n",
    "        }\n",
    "}\n",
    "\n",
    "lmks <- matrix(c(1:370), ncol=2, byrow=TRUE, \n",
    "               dimnames = list(vec,c(\"start\", \"end\")))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inter-landmark distance computing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "lineardists <- interlmkdist(landmarks, lmks)\n",
    "write.csv(lineardists, \"data/Linear_Measurements_all_with_NAs.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Removing outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "lineardists_out <- lineardists\n",
    "lineardists_out[8,130:137] <- NA #remove 2 last vertebrae from CHIglir BOR014\n",
    "lineardists_out[9,126:129] <- NA #remove last vertebrae from CHImajo\n",
    "lineardists_out[40,170:173] <- NA #remove last vertebrae from MALroth\n",
    "lineardists_out[48,130:133] <- NA #remove last vertebrae from MICspce"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Preparing subsets per vertebrae type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "lineardists_out <- as.data.frame(lineardists_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Removing and isolating femur measure (size proxy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizef <- lineardists_out[,\"FMRMID\"]\n",
    "sizef <- as.matrix(sizef)\n",
    "rownames(sizef) <- c(namesp)\n",
    "colnames(sizef) <- \"Proxy_size_raw\"\n",
    "TBLsizef <- as.matrix(sizef)\n",
    "lineardists_out <- lineardists_out[,-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing other size proxy (geometric mean of other measurements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "size_prox <- apply(lineardists_out, 1, function(x) prod(x, na.rm = T)^(1/(length(x)-sum(is.na(x)))))\n",
    "size_prox <- as.matrix(size_prox)\n",
    "rownames(size_prox) <- c(namesp)\n",
    "colnames(size_prox) <- \"Proxy_size_GM\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extracting specific measurments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableTPW <- lineardists_out[,seq(1, 184, 4)]\n",
    "tableMDW <- lineardists_out[,seq(2, 184, 4)]\n",
    "tableCTL <- lineardists_out[,seq(3, 184, 4)]\n",
    "tableCTH <- lineardists_out[,seq(4, 184, 4)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tail length proxy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "TailLength <- rowSums(tableCTL, na.rm = TRUE)\n",
    "TBLtaillength <- as.matrix(TailLength)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Attributing the longest vertebrae number for each specimen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "LongestVertebra <- apply(tableCTL,1,function(row) colnames(tableCTL)[which.max(row)])\n",
    "LongestVertebra <-sub(\"V\", \"\" , LongestVertebra)\n",
    "LongestVertebra <-sub(\"CTL\", \"\" , LongestVertebra)\n",
    "tableLV <- as.matrix(LongestVertebra)\n",
    "tableLV <- as.data.frame(tableLV)\n",
    "tableLV[,1]<-as.numeric(tableLV[,1])\n",
    "tableLV <- as.matrix(tableLV)\n",
    "rownames(tableLV) <- c(namesp)\n",
    "write.csv(tableLV, \"data/Longest_Vertebra_All.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Robusticity index (RI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Robusticity index for all specimens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "RI <- ((tableMDW)/(tableCTL))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mean robusticity index per species (multiple specimens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "RImean <- rowMeans(RI, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relative expansion of the transverse process index (TPEI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transverse process index for all specimens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEI <- ((tableTPW)/(tableCTL))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mean transverse process index per species (multiple specimens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEImean <- rowMeans(TPEI, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Per vertebral region species analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "lineardists_out$TV <- Vert_Count$Transition_Vertebra\n",
    "lineardists_out$LV <- Vert_Count$Longest_Vertebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subsetting per region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proximal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "proximal <- apply(lineardists_out, 1, function(x){\n",
    "  limit_prox <- as.numeric(x[\"TV\"])\n",
    "  vert_index <- sapply(strsplit(names(x), \"\"), function(x) paste(x[2:3], collapse = \"\"))\n",
    "  vert_index <- as.numeric(vert_index[-c(length(vert_index)-1, length(vert_index))])\n",
    "  limit_prox <- which(vert_index == limit_prox)\n",
    "  limit_prox <- limit_prox[length(limit_prox)]\n",
    "  proximal <- x[1:limit_prox]\n",
    "  names(proximal) <- NULL\n",
    "  return(proximal)\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transitional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "transitional <- apply(lineardists_out, 1, function(x){\n",
    "  limit_prox <- as.numeric(x[\"TV\"])\n",
    "  vert_index <- sapply(strsplit(names(x), \"\"), function(x) paste(x[2:3], collapse = \"\"))\n",
    "  vert_index <- as.numeric(vert_index[-c(length(vert_index)-1, length(vert_index))])\n",
    "  limit_prox <- which(vert_index == limit_prox)\n",
    "  limit_prox <- limit_prox[length(limit_prox)]\n",
    "  limit_trans <- as.numeric(x[\"LV\"])\n",
    "  limit_trans <- which(vert_index == limit_trans)\n",
    "  limit_trans <- limit_trans[length(limit_trans)]\n",
    "  transitional <- x[(limit_prox+1):limit_trans]\n",
    "  names(transitional) <- NULL\n",
    "  return(transitional)\n",
    "  \n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "distal <- apply(lineardists_out, 1, function(x){\n",
    "  vert_index <- sapply(strsplit(names(x), \"\"), function(x) paste(x[2:3], collapse = \"\"))\n",
    "  vert_index <- as.numeric(vert_index[-c(length(vert_index)-1, length(vert_index))])\n",
    "  limit_trans <- as.numeric(x[\"LV\"])\n",
    "  limit_trans <- which(vert_index == limit_trans)\n",
    "  limit_trans <- limit_trans[length(limit_trans)]\n",
    "  distal <- x[(limit_trans+1):(length(x)-2)]\n",
    "  names(distal) <- NULL\n",
    "  return(distal)\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assemble table in a list and fill gaps with NA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "regions <- list(proximal = proximal,\n",
    "                transitional = transitional,\n",
    "                distal = distal)\n",
    "regions_mat <- lapply(regions, function(x){\n",
    "  max_length <- max(sapply(x, length))\n",
    "  NA_rep <- max_length - sapply(x, length)\n",
    "  x <- lapply(x, function(y) c(y, rep(NA, max_length - length(y))))\n",
    "  x <- do.call(rbind, x)\n",
    "  return(x)\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a table with last 25% from the distal region per specimen (\"Last Quarter\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Isolating distal data of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "distL <- rowSums(!is.na(regions_mat$distal)) \n",
    "distLvrt <- (distL/4)\n",
    "distLQ <- (25*distLvrt)/100\n",
    "distLQ <- ceiling(distLQ)\n",
    "distLQm <- distLQ*4\n",
    "distLQv <- as.vector(distLQm)\n",
    "distalDF <- as.data.frame(regions_mat$distal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Creating a function extracting the last 25% from the distal region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "tail_func<-function(df, vec, species){\n",
    "  max_tail<-max(vec)\n",
    "  data_tail<-c(1:max_tail)\n",
    "  for (i in 1: nrow(df)){\n",
    "    vec_temp<-df[i,]\n",
    "    vec_temp<-vec_temp[!is.na(vec_temp)]\n",
    "    vec_temp<-rev(rev(vec_temp)[1:vec[i]])\n",
    "    if (length(vec_temp) < max_tail){\n",
    "      vec_temp<-c(vec_temp, rep(NA, (max_tail - length(vec_temp))))\n",
    "    }\n",
    "    data_tail<-cbind(data_tail, vec_temp)\n",
    "  }\n",
    "  data_tail<-t(data_tail)\n",
    "  data_tail<-data_tail[-1,]\n",
    "  rownames(data_tail)<-species\n",
    "return(as.data.frame(data_tail))  \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Apply the function to the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "LastQ <- tail_func(distalDF, distLQv, species)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Per region measure extraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Transverse process width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableTPWreg <- lapply(regions_mat, function(x){\n",
    "  x[, seq(1, ncol(x), by = 4)]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Mid-distal width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableMDWreg <- lapply(regions_mat, function(x){\n",
    "  x[, seq(1, ncol(x), by = 4)+1]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Centrum length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableCTLreg <- lapply(regions_mat, function(x){\n",
    "  x[, seq(1, ncol(x), by = 4)+2]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Centrum height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableCTHreg <- lapply(regions_mat, function(x){\n",
    "  x[, seq(1, ncol(x), by = 4)+3]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Last Quarter measure extraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Transverse process width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableTPWq <- LastQ[,seq(1,40,by=4)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Mid-distal width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableMDWq <- LastQ[,seq(1,40,by=4)+1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Centrum length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableCTLq <- LastQ[,seq(1,40,by=4)+2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Centrum height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "tableCTHq <- LastQ[,seq(1,40,by=4)+3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Per vertebrae region RI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proximal (All)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "RIprox <- ((tableMDWreg$proximal)/(tableCTLreg$proximal))*100 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proximal (Species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "RIproxmean <- rowMeans(RIprox, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transitional (All)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "RItran <- ((tableMDWreg$transitional)/(tableCTLreg$transitional))*100 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transitional (Species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "RItranmean <- rowMeans(RItran, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distal (All)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "RIdist <- ((tableMDWreg$distal)/(tableCTLreg$distal))*100 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distal (Species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "RIdistmean <- rowMeans(RIdist, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Last Quarter (All)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "RIlq <- ((tableMDWq)/(tableCTLq))*100 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Last Quarter (Species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "RIlqmean <- rowMeans(RIlq, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Per vertebrae region TPEI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proximal (All)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEIprox <- ((tableTPWreg$proximal)/(tableCTLreg$proximal))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proximal (Species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEIproxmean <- rowMeans(TPEIprox, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transitional (All)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEItran <- ((tableTPWreg$transitional)/(tableCTLreg$transitional))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transitional (Species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEItranmean <- rowMeans(TPEItran, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distal (All)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEIdist <- ((tableTPWreg$distal)/(tableCTLreg$distal))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Distal (Species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEIdistmean <- rowMeans(TPEIdist, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Last Quarter (All)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEIlq <- ((tableTPWq)/(tableCTLq))*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Last Quarter (Species)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "TPEIlqmean <- rowMeans(TPEIlq, na.rm = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing synthetic output table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merging all measurments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_synth <- cbind(vertnumb, TailLength, sizef, size_prox, RImean, RIproxmean, RItranmean, RIdistmean, RIlqmean, TPEImean, TPEIproxmean, TPEItranmean, TPEIdistmean, TPEIlqmean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Removing specimen not included in the phylogeny (Chiropodomys pusillus & Haeromys pusillus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "table_synth <- table_synth[!rownames(table_synth)==c(\"CHIpusi\"), ]\n",
    "table_synth <- table_synth[!rownames(table_synth)==c(\"HAEpusi\"), ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Isolating Nomenclature code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_names <- as.data.frame(as.matrix(str_split_fixed(rownames(table_synth), \"[.]\", n = 2)[,1]))\n",
    "colnames(data_names)<-\"Nomenclature\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merging with ecological data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nomenclature <- str_split_fixed(string = rownames(table_synth), pattern = \"[.]\", 2)[,1]\n",
    "\n",
    "table_synth <- aggregate(table_synth, list(Nomenclature), mean)\n",
    "\n",
    "colnames(table_synth)[1]<-\"Nomenclature\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_eco<-merge(data_names, data_eco, by = \"Nomenclature\")\n",
    "data_eco<-unique(data_eco)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_synth <- merge(data_eco, table_synth, by = \"Nomenclature\")\n",
    "rownames(table_synth)<-NULL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "write.table(table_synth, \"data/Complete_Table.tsv\", sep=\"\\t\")"
   ]
  }
 ],
 "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
}