---
title: "Introduction to mbend"
author: "Mohammad Ali Nilforooshan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to mbend}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

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

# Introduction to mbend

---

## Description

*Bending non-positive-definite (symmetric) matrices to positive-definite, using weighted and unweighted methods*

The `mbend` package is used for bending symmetric non-positive-definite matrices to positive-definite (PD). For a matrix to be invertible, it has to be PD. The methods of Jorjani et al. (2003) and Schaeffer (2014) are used in this package. The unweighted method of Schaeffer (2014) is also extended to weighted bending, with the possibility of choosing between unweighted and weighted bending.

## Application

Start with loading the package library:

```{r}
library(mbend)
```

Consider the following non-PD covariance matrix (Jorjani et al., 2003):

```{r}
V = matrix(nrow=5, ncol=5, c(
  100,  95,  80,  40,  40,
   95, 100,  95,  80,  40,
   80,  95, 100,  95,  80,
   40,  80,  95, 100,  95,
   40,  40,  80,  95, 100))
```

The following command can be used to bend matrix `V` to a PD matrix:

```{r}
bend(V)
```

The above command is equivalent to `bend(inmat=V)`, where `inmat` is the argument that takes the matrix to be bent, or the following command:

```{r, eval=FALSE}
bend(V, max.iter=10000, small.positive=0.0001, method="hj")
```

This runs the unweighted bending method of Jorjani et al. (2003) (`method="hj"`), with maximum 10000 number of iterations (`max.iter=10000`), and eigenvalues smaller than 0.0001 are replaced with this small positive value (`small.positive=0.0001`). Providing the default parameters, the corresponding arguments can be omitted (e.g., `bend(V)`).

The output object is a list of several items, listed below:

* bent : The bent `matrix`.
* init.ev : Eigenvalues of the initial (`inmat`) matrix.
* final.ev : Eigenvalues of the `bent` matrix.
* min.dev : `min(bent - inmat)`.
* max.dev : `max(bent - inmat)`.
* loc.min.dev : Location (indices) of `min.dev` element.
* loc.max.dev : Location (indices) of `max.dev` element.
* ave.dev : Average deviation (`bent - inmat`) of the upper triangle elements (excluding diagonal elements for correlation matrices).
* AAD : Average absolute deviation of the upper triangle elements (excluding diagonal elements for correlation matrices) of `bent` and `inmat`.
* Cor : Correlation between the upper triangle elements (excluding diagonal elements for correlation matrices) of `bent` and `inmat`.
* RMSD : Root of mean squared deviation of the upper triangle elements (excluding diagonal elements for correlation matrices) of `bent` and `inmat`.

There might be different precision involved with different elements of a non-PD matrix. In this case, a weighted bending is recommended. Jorjani et al. (2003) used the reciprocal of the number data points in common between pairs of variables, as weights. Considering the following matrix for the number of data points in common between variables (Jorjani et al., 2003):

```{r}
W = matrix(nrow=5, ncol=5, c(
  1000,  500,   20,   50,  200,
   500, 1000,  500,    5,   50,
    20,  500, 1000,   20,   20,
    50,    5,   20, 1000,  200,
   200,   50,   20,  200, 1000))
```

Matrix `V` is bent using the following command:

```{r}
bend(inmat=V, wtmat=W, reciprocal=TRUE)
```

Using `wtmat=1/W, reciprocal=FALSE`, the argument `reciprocal` could be omitted, because `FALSE` is the default parameter for the argument `reciprocal`. For the same reason, `max.iter=10000, small.positive=0.0001, method="hj"` are omitted, unless different parameters are provided to these arguments.

If there is high confidence about some elements of the non-PD matrix to remain unchanged after bending, the corresponding weights are set to zero. For example, to keep the first 2 &times; 2 block of `V` unchanged during the bending procedure:

```{r}
W2 = W; W2[1:2, 1:2] = 0
bend(V, W2, reciprocal=TRUE)
```

For weighted bending, we get extra statistics in the output:

* `w_gt_0` : Number of weight elements greater than 0, in the upper triangle of `wtmat` (for weighted bending).
* `wAAD` : Weighted `AAD` (for weighted bending).
* wCor : Weighted `Cor` (for weighted bending).
* wRMSD : Weighted `RMSD` (for weighted bending).

To bend `V` using the method of Schaeffer (2014):

```{r}
bend(inmat=V, method="lrs")
```

The method of Schaeffer (2014) does not require the argument `small.positive`, and this argument is ignored. This method is originally an unweighted bending method. However, in this package, it is extended to accommodate weighted bending (i.e., a combination of Schaeffer (2014) and Jorjani et al. (2003) methods). Weighted bending of `V` with reciprocals of `W` using the method of Schaeffer (2014):

```{r}
bend(V, W, reciprocal=TRUE, method="lrs")
```

Function `bend` automatically considers any matrix with all diagonal elements equal to one, as a correlation matrix. Considering the correlation matrix `V2` (`V` converted to a correlation matrix):

```{r}
V2 = cov2cor(V)
bend(V2, W, reciprocal=TRUE)
```

For correlation matrices, diagonal elements are not used in obtaining the statistics.
Because the argument `method` is not provided, the default parameter `"hj"` is used. To do the same using the method of Schaeffer (2014):

```{r}
bend(V2, W, reciprocal=TRUE, method="lrs")
```

Bending the same correlation matrix using the unweighted Schaeffer (2014):

```{r}
bend(V2, method="lrs")
```

---

## References

Jorjani, H., Klie. L., & Emanuelson, U. (2000). A simple method for weighted bending of genetic (co)variance matrices. *J. Dairy Sci.* 86(2): 677--679. <https://doi.org/10.3168/jds.S0022-0302(03)73646-7>

Schaeffer, L. R. (2014). Making covariance matrices positive definite. Available at: <http://animalbiosciences.uoguelph.ca/~lrs/ELARES/PDforce.pdf>
