---
title: "EDMA growth matrix"
author: "Peter Solymos"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{EDMA growth matrix}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup,include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
set.seed(429)
```

# Introduction

This tutorial explains how to calculate growth matrix based
on 4 EDMA data objects with homologous landmarks.

We will use data sets where landmarks were measured on specimens
of 2 different ages (E17.5 embryonic and newborn mice), 
taking a subset (`l`) of 10 landmarks:

```{r}
library(EDMAinR)

file_a1 <- system.file("extdata/growth/CZEM_wt_global.xyz",
    package="EDMAinR")
file_a2 <- system.file("extdata/growth/CZP0_wt_global.xyz",
    package="EDMAinR")

l <- c("amsph", "bas", "loci", "lpto", "lsqu",
        "lsyn", "roci", "rpto", "rsqu", "rsyn")

a1 <- read_xyz(file_a1)[l,,]
a2 <- read_xyz(file_a2)[l,,]
a1
a2
```

# Estimating the growth matrices

We first estimate the mean forms (no bootstrap replicates are necessary).

```{r}
fit_a1 <- edma_fit(a1, B=25)
fit_a2 <- edma_fit(a2, B=25)
```

Growth matrices ($GM$) are formed as pairwise Euclidean distances
between landmarks from EDMA fit objects using the estimated mean forms
from objects `a1` and `a2`.

The growth matrix is calculated as the ratio of form matrices
($FM$) from the numerator and denominator objects: $FDM(a1,a2) = FM(a2)/FM(a1)$.
`a2` is taken as the numerator, `a1` as the denominator.
We put the older sample (newborn) in the numerator spot 
and the younger sample (embryonic) in the denominator spot

```{r}
gm <- edma_gm(a1=fit_a1, a2=fit_a2, B=25)
gm
```

## Global T-test

The global testing is explained on the form difference page:

```{r fig.width=7,fig.height=5,out.width='100%'}
global_test(gm)
plot_test(gm)
```

## Local testing

The local testing is explained on the form difference page:

```{r}
head(confint(gm))
head(get_gm(gm))
head(get_gm(gm, sort=TRUE, decreasing=TRUE))
head(get_gm(gm, sort=TRUE, decreasing=FALSE))
```

The `plot_ci` function shows the pairwise differences and confidence intervals:

```{r fig.width=7,fig.height=7,out.width='100%'}
plot_ci(gm)
```


# Influential landmarks

Influence is calculated similarly to $FDM$:

```{r}
get_influence(gm)
plot(get_influence(gm))
```

# Ordination and clustering for specimens


```{r fig.width=7,fig.height=7,out.width='100%'}
plot_ord(gm)
```

The dendrogram leaves (specimen labels) are also colored by groups:

```{r fig.width=7,fig.height=7,out.width='100%'}
plot_clust(gm)
```


# Visualizing landmarks

The 2D and 3D plots produce a plot of the mean form
from the reference object ('prototype').
The color intensity for the landmarks (dots) is associated with the `Tdrop`
influence value (larger the difference, the more intensive the color;
red by default).
Lines between the landmarks represent distances.
We use the diverging palettes:
<1 differences are colored blue (1st half of the palette), >1 differences 
are colored red (2st half of the palette).

```{r fig.width=7,fig.height=5,out.width='100%'}
plot_2d(gm, cex=2)
```



```{r fig.width=7,fig.height=5,out.width='100%'}
library(rgl)
xyz <- plot_3d(gm)
text3d(xyz, texts=rownames(xyz), pos=1) # this adds names
decorate3d() # this adds the axes
rglwidget(width = 600, height = 600, reuse = FALSE)
```

# Growth difference matrix

Growth difference matrix ($GDM$) is calculated as 
$GDM(A1,A2,B1,B2) = GM(B1,B2) / GM(A1,A2)$.

We will use two Crouzon mutant samples, same age groups as 
for the unaffected samples (embryonic and newborn):

```{r}
file_b1 <- system.file("extdata/growth/CZEM_mut_global.xyz",
    package="EDMAinR")
file_b2 <- system.file("extdata/growth/CZP0_mut_global.xyz",
    package="EDMAinR")

b1 <- read_xyz(file_b1)[l,,]
b2 <- read_xyz(file_b2)[l,,]
b1
b2

fit_b1 <- edma_fit(b1, B=25)
fit_b2 <- edma_fit(b2, B=25)
```

Growth matrices ($GM$) are formed as pairwise Euclidean distances
between landmarks from EDMA fit objects using the estimated mean forms
from objects `a1` and `a2`.

The growth matrix is calculated as the ratio of form matrices
($FM$) from the numerator and denominator objects: $FDM(a1,a2) = FM(a2)/FM(a1)$.
`a2` is taken as the numerator, `a1` as the denominator.
We put the older sample (newborn) in the numerator spot 
and the younger sample (embryonic) in the denominator spot

```{r}
gdm <- edma_gdm(a1=fit_a1, a2=fit_a2, b1=fit_b1, b2=fit_b2, B=25)
gdm

global_test(gdm)
plot_test(gdm)
plot_ci(gdm)

plot_ord(gdm)
plot_clust(gdm)

plot_2d(gdm)
```

```{r fig.width=7,fig.height=5,out.width='100%'}
xyz <- plot_3d(gdm)
text3d(xyz, texts=rownames(xyz), pos=1) # this adds names
decorate3d() # this adds the axes
rglwidget(width = 600, height = 600, reuse = FALSE)
```