---
title: "Construction of the M-tests"
date: "`r Sys.Date()`"
author: Sermad Abbas, Barbara Brune, Roland Fried 
output: 
  rmarkdown::html_document:
    toc: true
bibliography: ../inst/REFERENCES.bib
csl: csda.csl
vignette: >
  %\VignetteIndexEntry{Construction of the M-tests}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
  
```{r, include = FALSE}
knitr::opts_chunk$set(
cache = TRUE,
collapse = TRUE,
comment = "#>"
)
```

# Introduction 
In this vignette, we will briefly describe and motivate how we constructed the test statistics used by the function `m_test` and how it derives a test decision.

# Asymptotic distribution of the M-test statistics
For a more detailed description of the asymptotic behaviour of M-estimators, we refer to @MarMarYoh19robu[p. 36ff.], which is the main reference for the following motivation.

We consider two independent samples $X_1, \ldots, X_m$ and $Y_1, \ldots, Y_n$ of i.i.d. random variables which are symmetrically distributed with variances $\sigma^2_X$ and $\sigma^2_Y$.

For M-estimators $\hat{\mu}_X$ and $\hat{\mu}_Y$ with a $\psi$-function $\psi$, it can be shown under these conditions that
\begin{align*}
\sqrt{m} \cdot \left(\hat{\mu}_X - \mu_X\right) \overset{\text{asympt.}}{\sim} \mathcal{N}\left(0, \sigma_X^2 \cdot \nu_X\right) \quad
\text{and} \quad
\sqrt{n} \cdot \left(\hat{\mu}_Y - \mu_Y\right) \overset{\text{asympt.}}{\sim} \mathcal{N}\left(0, \sigma_Y^2 \cdot \nu_Y\right),
\end{align*}

where $\mu_X \in \mathbb{R}$ and $\mu_Y \in \mathbb{R}$ are the values for which
\begin{align*}
\text{E}\left(\psi\left(\frac{X - \mu_X}{\sigma_X}\right)\right) = 0 \quad \text{and} \quad \text{E}\left(\psi\left(\frac{Y - \mu_Y}{\sigma_Y}\right)\right) = 0,
\end{align*}

and

\begin{align*}
\nu_X = \frac{\text{E}\left(\psi\left(\frac{X - \mu_X}{\sigma_X}\right)^2\right)}{\left(\text{E}\left(\psi'\left(\frac{X - \mu_X}{\sigma_X}\right)\right)\right)^2} \quad \text{and} \quad
\nu_Y = \frac{\text{E}\left(\psi\left(\frac{Y - \mu_Y}{\sigma_Y}\right)^2\right)}{\left(\text{E}\left(\psi'\left(\frac{Y - \mu_Y}{\sigma_Y}\right)\right)\right)^2}.
\end{align*}

From this, it follows that
\begin{align*}
\hat{\mu}_X \overset{\text{asympt.}}{\sim} \mathcal{N}\left(\mu_X, \frac{\sigma^2_X \cdot \nu_X}{m}\right) \quad \text{and} \quad
\hat{\mu}_y \overset{\text{asympt.}}{\sim} \mathcal{N}\left(\mu_Y, \frac{\sigma^2_Y \cdot \nu_Y}{n}\right),
\end{align*}

implying

\begin{align*}
\frac{\hat{\mu}_X - \hat{\mu}_Y - \left(\mu_X - \mu_Y\right)}{\sqrt{\frac{n \cdot \sigma^2_X \cdot \nu_X + m \cdot \sigma^2_Y \cdot \nu_Y}{m \cdot n}}} \overset{\text{asympt.}}{\sim} \mathcal{N}\left(0, 1\right).
\end{align*}

In order to use this statistic as a test statistic for our M-tests, we need to estimate $\sigma_X$, $\sigma_Y$, $\nu_X$, and $\nu_Y$.
We use the $\tau$-scale estimator [@MarZam02robu] to estimate $\sigma^2_X$ and $\sigma^2_Y$ by $\hat{\sigma}_X^2$ and $\hat{\sigma}_Y^2$ robustly and estimate $\nu_X$ and $\nu_Y$ by

\begin{align*}
\hat{\nu}_X = \frac{\frac{1}{m} \sum_{i = 1}^m \psi\left(\frac{X_i - \hat{\mu}_X}{\hat{\sigma}_X}\right)^2}{\left(\frac{1}{m} \sum_{i = 1}^m \psi'\left(\frac{X_i - \hat{\mu}_X}{\hat{\sigma}_X}\right)\right)^2} \quad
\text{and} \quad
\hat{\nu}_Y = \frac{\frac{1}{n} \sum_{j = 1}^n \psi\left(\frac{Y_j - \hat{\mu}_Y}{\hat{\sigma}_Y}\right)^2}{\left(\frac{1}{n} \sum_{j = 1}^n \psi'\left(\frac{Y_j - \hat{\mu}_Y}{\hat{\sigma}_Y}\right)\right)^2}.
\end{align*}

Under the previous considerations, the test statistic of the M-tests we implemented in the package is given by

\begin{equation*}
\frac{\hat{\mu}_X - \hat{\mu}_Y - \Delta}{\sqrt{\frac{n \cdot \hat{\sigma}^2_X \cdot \hat{\nu}_X + m \cdot \hat{\sigma}^2_Y \cdot \hat{\nu}_Y}{m \cdot n}}} \overset{\text{asympt.}}{\sim} \mathcal{N}\left(0, 1\right),
\end{equation*}

where $\Delta = \mu_X - \mu_Y$ is the location difference between both distributions.

The M-tests are implemented in the function `m_test`. 
More details on the usage of the function can be found in the vignette [Getting started with `robnptests`](robnptests.html).
Inside `m_test`, we use the function `scaleTau2` from the `R` package [robustbase](https://cran.r-project.org/package=robustbase) [@MaeRouCro22robu] to compute the $\tau$-scale estimates for the samples.

# Simulation results
The following figure shows the simulated test sizes from a small simulation study with 1000 replications, where we applied the M-tests with different $\psi$-functions to samples from the $\mathcal{N}(0, 1)$-distribution, the $t_2$-distribution, and the $\chi^2_3$-distribution. 
We chose the significance level $\alpha = 0.05$. 
The results are shown in the following figure.

Under the $\mathcal{N}(0, 1)$- and the $t_2$-distribution we make similar observations: 
For equal sample sizes $m = n \geq 30$, the simulated test size is quite close to the the specified value of $\alpha$. When $m \neq n$, it seems to be important that both values are rather large and do not deviate too much from each other.
Otherwise, the tests may become very anti-conservative. 
In general, the three test statistics lead to similar results for the considered sample sizes.

Under the $\chi^2_3$-distribution, all tests are anti-conservative. 
While there seems to be some improvement when the sample sizes become larger, the estimated sizes are still rather far away from 0.05. 
A reason might be that the asymptotic variance we use is only a good approximation for symmetric distributions [@MarMarYoh19robu, p. 37].

```{r, echo = FALSE, fig.height = 10, fig.width = 8}
load(file.path("..", "R", "sysdata.rda"))

# M-estimators to plot ----
estimators <- unique(results.normal$estimator)

# Plot results ----
par(mfrow = c(3, 1))

# Normal distribution ----

# Equal sample sizes
for (i in seq_along(estimators)) {
  inds <- which(results.normal[, "estimator"] == estimators[i])
  p.value <- results.normal[inds, "size"]
  std.err <- results.normal[inds, "std.error"]
  x <- 1:length(inds)
  m <- results.normal[inds, "m"]
  n <- results.normal[inds, "n"]
  
  cols <- c("black", "red", "green")
  
  if (i == 1) {
    plot(x, p.value,
         ylim = c(0, 0.18),
         xlab = "Sample sizes",
         ylab = "Size of the test",
         main = "N(0, 1)-distribution",
         col = cols[i],
         xaxt = "n",
         yaxt = "n"
    )
    axis(side = 1, at = 1:10, labels = unique(paste0("(", m, ", ", n, ")")))
    axis(side = 2, at = seq(0.025, 0.175, by = 0.025), labels = c("", "0.05", "", "0.1", "", "0.15", ""))
  } else {
    points(x, p.value,
           col = cols[i],
    )
  }
  
  abline(h = 0.05, lty = "dashed")
}

# t2-distribution ----
for (i in seq_along(estimators)) {
  inds <- which(results.t[, "estimator"] == estimators[i])
  p.value <- results.t[inds, "size"]
  std.err <- results.t[inds, "std.error"]
  x <- 1:length(inds)
  m <- results.t[inds, "m"]
  n <- results.t[inds, "n"]
  
  cols <- c("black", "red", "green")
  
  if (i == 1) {
    plot(x, p.value,
         ylim = c(0, 0.18),
         xlab = "Sample sizes",
         ylab = "Size of the test",
         main = "t(2)-distribution",
         col = cols[i],
         xaxt = "n",
         yaxt = "n"
    )
    axis(side = 1, at = 1:10, labels = unique(paste0("(", m, ", ", n, ")")))
    axis(side = 2, at = seq(0.025, 0.175, by = 0.025), labels = c("", "0.05", "", "0.1", "", "0.15", ""))
  } else {
    points(x, p.value,
           col = cols[i],
    )
  }
  
    abline(h = 0.05, lty = "dashed")
}

# Chi-square distribution ----
for (i in seq_along(estimators)) {
  inds <- which(results.chi[, "estimator"] == estimators[i])
  p.value <- results.chi[inds, "size"]
  std.err <- results.chi[inds, "std.error"]
  x <- 1:length(inds)
  m <- results.chi[inds, "m"]
  n <- results.chi[inds, "n"]
  
  cols <- c("black", "red", "green")
  
  if (i == 1) {
    plot(x, p.value,
         ylim = c(0, 0.18),
         xlab = "Sample sizes",
         ylab = "Size of the test",
         main = "Chi^2(3)-distribution",
         col = cols[i],
         xaxt = "n",
         yaxt = "n"
    )
    axis(side = 1, at = 1:10, labels = unique(paste0("(", m, ", ", n, ")")))
    axis(side = 2, at = seq(0.025, 0.175, by = 0.025), labels = c("", "0.05", "", "0.1", "", "0.15", ""))
  } else {
    points(x, p.value,
           col = cols[i],
    )
  }
  
  legend("bottomleft", legend = c("Huber", "Hampel", "Bisquare"), col = 1:3, pch = 1)
  
  abline(h = 0.05, lty = "dashed")
}
```

Based on these results, we discourage using the tests for asymmetric distributions.
For symmetric distributions, the asymptotic test should only be used for large samples. In all other cases, the randomization or permutation test might be preferable.

# Session Info
```{r}
library(robnptests)

sessionInfo()
```


# References

