Geomorph_regression.Rmd

---
title: "geomorph_regression"
author: "SlicerMorph"
date: "4/02/2021"
output:
  html_document: default
  pdf_document: default
---
# How to perform a regression using geomorph in R from SlicerMorph outputs from the GPA module

This script is based on geomorph 3.3.2.
In this example we are using the output from the gorilla skull LMs data distributed by SlicerMorph. 
The first step is to load the geomorph library and the parser convenience function to pull all the files and analytical settings
These functions eventually will be turned into an R package. 

```{r setup}
library(geomorph)

source("https://raw.githubusercontent.com/muratmaga/SlicerMorph_Rexamples/main/log_parser.R")
```

Next we need to point to the location of the analysis.log file that was saved by SlicerMorph's GPA module either via coding the path OR interactively. 
If you are knitting the document, make sure to hard code the path. 
```{r load log file}
# coding the path to log file example 
# SM.log.file="C:\\temp\\RemoteIO\\gorilla_patches\\merged\\2021-04-05_20_51_17\\analysis.log"

# interactively choose log file
SM.log.file <- file.choose()

SM.log <- parser(SM.log.file)
```

The SM.log.file contains pointers to all relevant data files. [For details of what's provided see](https://github.com/muratmaga/SlicerMorph_Rexamples/blob/main/README.md#slicermorph_rexamples).
Example below reads the PC scores are calculated by GPA module into a data matrix. 

```{r pointers}
head(SM.log)

SM.output <- read.csv(file = paste(SM.log$output.path, 
                                SM.log$OutputData, 
                                sep = "/"))
                                
SlicerMorph.PCs <- read.table(file = paste(SM.log$output.path, 
                                        SM.log$pcScores, 
                                        sep="/"), 
                              sep = ",", header = TRUE, row.names = 1)
head(SlicerMorph.PCs)
```

### Check the number of landmarks used in the analysis and extract other variables
```{r LM number}
PD <- SM.output[,2]
if (!SM.log$skipped) {
  no.LM <- SM.log$no.LM 
  } else {
    no.LM <- SM.log$no.LM - length(SM.log$skipped.LM)
  }
PD
no.LM
```

### Reformat the GPA aligned coords into 3D LM array and apply sample names
```{r reformat}
Coords <- arrayspecs(SM.output[, -c(1:3)], 
                    p = no.LM, 
                    k = 3)

dimnames(Coords) <- list(1:no.LM, 
                        c("x","y","z"),
                        SM.log$ID)
```

### Running Stats on GPA outputs
We first construct a geomorph data frame withe data imported from SlicerMorph and fit a model to SlicerMorph's GPA aligned coordinates and centroid sizes
```{r fit model}
gdf <- geomorph.data.frame(Size = SM.output$centeroid, Coords = Coords)
fit.slicermorph <- procD.lm(Coords ~ Size, data = gdf, print.progress = FALSE)
```

This second part of the script uses the raw LM coordinates directly into R/geomorph,aligns them with gpagen(), applies PCA, and builds the same allometric regression model and then compares it to the results obtained above. 
```{r geomorph gpa comparison}
gpa <- gpagen(SM.log$LM)
pca <- gm.prcomp(gpa$coords)
geomorph.PCs <- pca$x
gdf2 <- geomorph.data.frame(size = gpa$Csize, coords = gpa$coords)
fit.rawcoords <- procD.lm(coords ~ size, data = gdf2, print.progress = FALSE)
```

Due to arbitrary rotations, we cannot compare procrustes coordinates directly, instead we look centroid sizes, procD, PC scores, and allometric regression model summary. 

geomorph does not report each sample's procD to the consensus shape, so we need to calculate PDs.
```{r}
pd <- function(M, A) sqrt(sum(rowSums((M-A)^2)))
geomorph.PD <- NULL
for (i in 1:length(SM.log$files)) 
  geomorph.PD[i] <- pd(gpa$consensus, gpa$coords[,, i])
```

### Now we can start to compare procrustes variables
*1. Centroid Size*
```{r}
plot(gpa$Csize, SM.output$centeroid, 
     pch = 20, ylab = 'SlicerMorph', 
     xlab = 'geomorph', main = "Centroid Size")
cor(gpa$Csize, SM.output$centeroid)
```

*2. Procrustes Distance of sample to their respective mean*
```{r}
plot(geomorph.PD, SM.output[,2], 
     pch = 20, ylab = 'SlicerMorph', 
     xlab = 'geomorph', main = "Procrustes Distance")
cor(geomorph.PD, SM.output[,2])
```

*3. We only plot the first two PCs but correlations reported up to 10*
* Keep in mind that PCA signs are arbitrary. 
```{r}
plot(geomorph.PCs[,1], SlicerMorph.PCs[,1], 
     pch = 20, ylab = 'SlicerMorph', 
     xlab = 'geomorph', main = "PC1 scores")

plot(geomorph.PCs[,2], SlicerMorph.PCs[,2], 
     pch = 20, ylab = 'SlicerMorph', 
     xlab = 'geomorph', main = "PC2 scores")

for (i in 1:10) print(cor(SlicerMorph.PCs[, i], 
                            geomorph.PCs[, i]))
```

*4. Finally, compare allometric regression models from SLicerMorph aligned coordinates to the one that used raw LM coordinates and gpa alignment from geomorph.*
```{r}
summary(fit.slicermorph)

summary(fit.rawcoords)
```