---
title: "EDMA form matrix"
author: "Peter Solymos"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{EDMA form 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 use the **EDMAinR** package
to fit a nonparametric model to landmark data and estimate
the mean form.

Let's load the package and read in an example data set:

```{r}
library(EDMAinR)
file <- system.file("extdata/crouzon/Crouzon_P0_Global_MUT.xyz", 
    package="EDMAinR")
x <- read_xyz(file)
```

# Nonparametric fit

We use the `edma_fit` function to run the nonparametric estimation.
The `B` argument is the number of bootstrap replicates to take
(resampling the $n$ specimens from object `x` with replacement).
We chose a small number here, but `B` is ideally much higher
(i.e. `B=99`):

```{r}
fit <- edma_fit(x, B=9)
fit
```

The nonparametric estimator gives the mean form matrix ($\hat{M}$)
and $\hat{\Sigma}_{K}^\ast$. We can extract these from the fitted
model object (`fit`) using the `Meanform` and `SigmaKstar` functions:

```{r}
str(Meanform(fit))
str(SigmaKstar(fit))
```

> *Note*: we are using the `str` function to peek into the - otherwise
> huge - objects, just printing the `str`ucture.

# Form matrix

We can extract stacked Euclidean distances (the stacked form matrix $FM$)
with the `get_fm` function. The functions allows to sort
the table based on decreasing or increasing order by distance:

```{r stackeddist}
head(get_fm(fit))
head(get_fm(fit, sort=TRUE, decreasing=TRUE))
head(get_fm(fit, sort=TRUE, decreasing=FALSE))
```

The columns `row` and `col` refer to the landmark pairs.
`dist` gives the Euclidean distance between the mean form corrdinates
for the two landmarks. `lower` and `upper` indicates
the 95% confidence limits based on the bootstrap distribution
of the pairwise distances (from `B` replicates).
The confidence `level` can be changed as
`get_fm(fit, level=0.9)`.

# Visualizing the mean form

The 2D plot gives the projection of the mean form
using multidimensional scaling (it uses $FM$ based on $\hat{M}$).
The dot size is proportional to `SigmaKstar(fit)` diagonal elements
(scaling can be chosen by the `cex` argument).
The `plot_2d` function returns the plotting coordinates,
which can be used to add the landmark names:

```{r fig.width=7,fig.height=5,out.width='60%'}
xy <- plot_2d(fit, cex=2)
text(xy, labels=rownames(xy), pos=1)
```

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

See `?rgl::plot3d` for options to modify the 3D plot.

# Changing the default colors

It is possible to supply your own color values to the
`plot_2d` and `plot_3d` functions to allow full control:

```{r fig.width=7,fig.height=5,out.width='60%'}
col <- rainbow(nrow(xy))
plot_2d(fit, col=col, cex=2)
```

```{r fig.width=7,fig.height=5,out.width='100%'}
plot_3d(fit, col=col)
rglwidget(width = 600, height = 600, reuse = FALSE)
```

An easier way, however, is to provide a color palette
from the choices listed on the help page of the function
`?hcl.colors`. The available diverging color palettes are:


```{r}
hcl.pals(type = "diverging")
```

The **EDMAinR** default palettes are stored as 'options'.

```{r}
getOption("edma_options")
```

The default palettes produces these colors:

```{r fig.width=4,fig.height=2,out.width='50%'}
plot_edma_colors(101)
```

Here is how to change the default palette:

```{r}
op <- options("edma_options" = list(
    diverging = "Green-Orange",
    qualitative = "Dark 2"))
op
```

This produces the following palettes:

```{r echo=FALSE,fig.width=4,fig.height=2,out.width='50%'}
plot_edma_colors(101)
```

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

```{r fig.width=7,fig.height=5,out.width='100%'}
plot_3d(fit)
rglwidget(width = 600, height = 600, reuse = FALSE)
```

> *Note*: we are only using the 'high' end of the palette
> because the dots are scaled by standard deviation like measure
> based on `sqrt(diag(SigmaKstar(fit)))`.

We can reset the defaults using the `op` list that we saved before:

```{r}
options(op)
```