---
title: "Details of Hypothesis Testing"
package: mmrm
bibliography: '`r system.file("REFERENCES.bib", package = "mmrm")`'
csl: '`r system.file("jss.csl", package = "mmrm")`'
output:
  rmarkdown::html_vignette:
          toc: true
vignette: |
  %\VignetteIndexEntry{Details of Hypothesis Testing}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options:
  chunk_output_type: console
  markdown:
    wrap: 72
---

# Introduction to Type I/II/III Hypothesis Testing

Type I/II/III hypothesis testing originates from analysis of variance, however, 
@goodnight1980tests suggests viewing this issue as hypothesis testing of fixed effects
in linear models, which avoid the "usual assumptions" and defined three types of estimable
functions.
Given that type I estimable functions are dependent on the specific order of used covariates,
we focus on the type II and III testing here.

For MMRM, obtaining estimable functions is equivalent to obtaining the contrasts for linear hypothesis tests.
For more details, see @goodnight1980tests.

# Contained Effect

Before we discuss the testing, it is important that we understand the concept of "contained effect".
The definition of "contained effect" is:

For effect $E_2$, it is said to contain $E_1$ if

1. the involved numerical covariate(s) are identical for these two effects
1. all the categorical covariate(s) associated with $E_1$ are also associated with $E_2$

For example, for the following model formula

$Y = A + B + A*B$

using $E_A$, $E_B$ and $E_{A*B}$ to denote the effect for $A$, $B$ and $A*B$ respectively,
we have

1. If $A$ and $B$ are both categorical, then $E_{A*B}$ contains $E_{A}$ and $E_{B}$.
1. If $A$ and $B$ are both numerical, then $E_{A*B}$ do not contain $E_{A}$ or $E_{B}$.
1. If $A$ is categorical and $B$ is numerical, then $E_{A*B}$ contains $E_{B}$ but not $E_{A}$.

# Type II Hypothesis Testing

For type II hypothesis testing, the contrasts can be obtained through the following steps:

Create a contrast matrix $L$, with rows equal to the number of parameters associated with the effect,
columns equal to the number of all parameters, *including aliased parameters*.

1. All columns of $L$ associated with effect not containing $E_1$ (except $E_1$) are set to 0.
1. The submatrix of $L$ associated with $E_1$ is set to $(X_1^\top M X_1)^{-}X_1^\top M X_1$
1. Each of the remaining submatrices of L associated with an effect $E_2$ that contains $E_1$ is $(X_1^\top M X_1)^{-}X_1^\top M X_2$
1. Columns associated with aliased parameters are dropped.

where $X_0$ stands for columns of $X$ whose associated effect do not contain $E_1$,
$X_1$ stands for columns of $X$ associated with $E_1$,
$X_2$ stands for columns of $X$ whose associated effect contains $E_1$,
$M = I - X_0(X_0^\top X_0)^{-}X_0^\top$,
and $Z^{-}$ stands for the [g2 inverse](https://blogs.sas.com/content/iml/2018/11/21/generalized-inverses-for-matrices.html) of $Z$.

Note: Here we do not allow for singularity in general, so the g2 inverse is just the usual inverse.
Thus we can use an identity matrix as the submatrix in step 2.

Using `fev_data` to create our example, we have

```{r}
library(mmrm)
fit <- mmrm(FEV1 ~ ARMCD + RACE + ARMCD * RACE + ar1(AVISIT | USUBJID), data = fev_data)
```

For this given example, we would like to test the effect of `RACE`, $E_{RACE}$.
We initialize the contrast matrix with

$$
\begin{matrix} 
\mu & ARMCD & RACE & RACE & ARMCD:RACE & ARMCD:RACE\\
0 & 0 & 0 & 0 & 0 & 0\\ 
0 & 0 & 0 & 0 & 0 & 0\\ 
\end{matrix}
$$

Please note that this is a $2\times 6$ matrix, the rows is the number of coefficients that is associated with `RACE`,
and the columns is the total number of coefficients `r length(coef(fit))`.

Then following step 2, we filled in the identity matrix in the corresponding submatrix

$$
\begin{matrix} 
\mu & ARMCD & RACE & RACE & ARMCD:RACE & ARMCD:RACE\\
0 & 0 & 1 & 0 & 0 & 0\\ 
0 & 0 & 0 & 1 & 0 & 0\\ 
\end{matrix}
$$

In the last step, for the last 2 columns, we fill the values with the calculated result.

```{r}
x <- component(fit, "x_matrix")
x0 <- x[, c(1, 2)]
x1 <- x[, c(3, 4)]
x2 <- x[, c(5, 6)]
m <- diag(rep(1, nrow(x))) - x0 %*% solve(t(x0) %*% x0) %*% t(x0) # solve is used because the matrix is inversible
sub_mat <- solve(t(x1) %*% m %*% x1) %*% t(x1) %*% m %*% x2
sub_mat
```

$$
\begin{matrix} 
\mu & ARMCD & RACE & RACE & ARMCD:RACE & ARMCD:RACE\\
0 & 0 & 1 & 0 & 0.42618692 & 0.02751985\\ 
0 & 0 & 0 & 1 & -0.04372702 & 0.58570969\\ 
\end{matrix}
$$

# Type III Hypothesis Testing

For Type III hypothesis testing, we assume that the hypothesis tested should be the same for all
designs with the same general form of estimable functions.

An important consideration for Type III testing is that it should be performed with orthogonal contrasts for the categorical covariates. In R this can be achieved via using e.g. the `contr.sum` contrasts for the categorical covariates. 
However, the default in R are treatment contrasts, via `contr.treatment`. Therefore, to ensure that Type III tests are performed correctly, independently of which contrasts happen to be used in the model fitting, we follow this algorithm:

1. Based on the user specification in the data frame we construct a design matrix, say $X_t$ (this could be using all `contr.treatment` or a mixture of contrasts for factor variables).
1. We also construct the design matrix with `contr.sum` for all factor variables, say $X_c$.
1. We construct the Type III contrast matrices $L_{c,j}$ for all terms $j$ for the coefficient estimate based on the $X_c$ design matrix:
  - All columns of $L_{c,j}$ associated with effect not being $E_j$ are set to 0.
  - The submatrix of $L_{c,j}$ associated with $E_j$ is set to the identity matrix.
1. We calculate a correction matrix $C = (X_c^\top X_c)^{-1} X_c^\top X_t$ and obtain corrected Type III contrast matrices subsequently as $L_{t,j} = L_{c,j} C$.
1. We perform the linear hypothesis tests with the contrast matrices $L_{t,j}$ for our coefficient estimate.

Note that an intercept term should be included in a model when performing Type III tests: otherwise one of the categorical covariates will have an arbitrary high test statistic value. This behavior is different from SAS, but consistent with other `car::Anova()` results e.g. from linear or generalized linear models. In order to guide the user, a warning message is issued when Type III tests are requested for models without an intercept.

# Hypothesis Testing in SAS

In `PROC MIXED` we can specify `HTYPE=1,2,3` to enable these hypothesis testings, and we can use option `e1 e2 e3` to
view the estimable functions.

## Special Considerations

### Reference Levels

For `PROC MIXED`, it is important that the categorical values have the correct order, especially
for structured covariance (otherwise the correlation between visits can be messed up).
We usually use `CLASS variable(REF=)` to define the reference level.
However, SAS will put the reference level in the last, so if the visit variable is included in fixed effect,
there can be issues with the comparison between SAS and R, as in R we still use the first level as our reference.

When we conduct the hypothesis testing, we will also face this issue.
If reference level is different, the testing we have can be different, e.g. if we have 3 levels, $l_1$, $l_2$ and $l_3$.
In theory we can test either $l_2 - l_1$, $l_3 - l_1$ (using $l_1$ as reference), or test $l_1 - l_3$, $l_2 - l_3$ (using $l_3$ as reference).
But the result will be slightly different.
However, if we have identical tests the result will be closer.
The following example illustrate this.

#### Example of Reference Levels

Assume in SAS we have the following model

```
PROC MIXED DATA = fev cl method=reml;
  CLASS AVISIT(ref = 'VIS4') ARMCD(ref = 'PBO') USUBJID;
  MODEL FEV1 = ARMCD AVISIT ARMCD*AVISIT / ddfm=satterthwaite htype=3;
  REPEATED AVISIT / subject=USUBJID type=ar(1);
RUN;
```

And in R we have the following model
```{r eval=FALSE}
fit <- mmrm(FEV1 ~ ARMCD + AVISIT + ARMCD * AVISIT + ar1(AVISIT | USUBJID), data = fev_data)
Anova(fit)
```

Note that for `AVISIT` there are 4 levels, `VIS1`, `VIS2`, `VIS3` and `VIS4`.
In SAS `VIS4` is the reference, while in R `VIS1` is the reference.

They both use the following matrix to test the effect of `AVISIT`
(ignoring other part of contrasts that is not associated with `AVISIT`)

$$
\begin{matrix} 
1 & 0 & 0\\ 
0 & 1 & 0\\
0 & 0 & 1\\
\end{matrix}
$$

But given that the reference level is different, we are testing different things.
In SAS we are testing `VIS1 - VIS4`, `VIS2 - VIS4` and `VIS3 - VIS4`,
while in R we are testing `VIS2 - VIS1`, `VIS3 - VIS1` and `VIS4 - VIS1`.

To correctly test `VIS1 - VIS4`, `VIS2 - VIS4` and `VIS3 - VIS4` in R, we need to update the contrast
to the following matrix (ignoring the other part) to get a closer result.

$$
\begin{matrix} 
0 & 1 & 0\\ 
0 & 0 & 1\\
-1 & -1 & -1\\
\end{matrix}
$$


### Why Model Covariate Order Changes My Testing in SAS?

For type II/III testing, it is true that the order of covariate should not affect the result.
However, if there is an interaction term, the reference level can be changed.
With different reference levels the result will be slightly different.

### Why `mmrm` Gives More Covariates Than SAS?

For `PROC MIXED`, if a model is specified as follows

```
MODEL FEV1 = ARMCD SEX ARMCD * SEX ARMCD * FEV1_BL
```

Then it is equivalent to the following `mmrm` model

```
FEV1 ~ ARMCD * SEX + ARMCD * FEV1_BL - FEV1_BL
```

or

```
FEV1 ~ ARMCD * SEX + ARMCD:FEV1_BL
```

Because SAS will not include the covariate `FEV1_BL` by default, unless manually added.
However in R we will add `FEV1_BL` by default.
Please note that in this example we exclude the covariance structure part.

### Excluding columns due to collinearity

`mmrm::mmrm()` (which employs `qr()`) and `PROC MIXED` have similar strategies to handle collinearity of model parameters. Both prioritize earlier columns in the design matrix, excluding a column if it is a linear combination of earlier ones. Whereas `mmrm::mmrm()` completely excludes such parameters from the model, `PROC MIXED` retains them but assigns a coefficient of 0. The presence of these columns is necessary for the accurate calculation of type-II contrast matrices. 

Therefore, an `mmrm` object contains the full design matrix including aliased columns in a component called `x_matrix_complete`. This matrix is utilized to calculate the contrast matrices when running a type-II ANOVA. Once all values are correctly calculated, the final step in contrast matrix creation is to drop the columns corresponding to the aliased parameters.

### Intercept-free models

`stats::model.matrix()` and `PROC MIXED` handle intercept-free models in different ways. This has important implications for the creation of the contrast matrices used for type-II ANOVAs.

#### Intercept-free Models with `stats::model.matrix()`

When the default method of `stats::model.matrix()` is creating a design matrix for an intercept-free model, it considers the terms sequentially. When it encounters the first term to contain a categorical variable, it will include the reference level when creating this term's design matrix columns. Therefore, this term de facto contains the model's intercept.

#### Intercept-free Models with `PROC MIXED`

`PROC MIXED`, on the other hand, treats all terms the same regardless of the presence or absence of an intercept: it always considers all factor levels for all categorical variables when creating design matrix columns (see the **Main Effects** and **Implications of the Non-Full-Rank Parameterization** sections [here](https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect023.htm)). Therefore, all terms involving categorical variables effectively share the model's intercept. 

#### Type-II Contrast Matrices in Intercept-free Models

To accurately run a type-II ANOVA when all levels of a categorical variable are represented in a design matrix, one of the levels of each categorical variable must contain subtractive (i.e., negative) values its corresponding contrast matrix columns. Since `PROC MIXED` includes all factor levels of all categorical variables, the incorporation of negative values is a commonplace occurrence. Negative values are invariably assigned to the last columns, which are always the reference levels' columns. 

However, in `mmrm::mmrm()`, reference level columns only exist in intercept-free models for the first term containing a categorical variable. Thus, `mmrm`'s helper functions incorporate a special handling step for such first-categorical-terms in which their last contrast matrix columns are assigned straight -1s.

