Type I/II/III hypothesis testing originates from analysis of variance, however, Goodnight (1980) 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 Goodnight (1980).
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
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
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.
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 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
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 6.
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.
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
## ARMCDTRT:RACEBlack or African American
## RACEBlack or African American 0.42618692
## RACEWhite -0.04372702
## ARMCDTRT:RACEWhite
## RACEBlack or African American 0.02751985
## RACEWhite 0.58570969
\[ \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} \]
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. It can be defined with 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.
Here equate the values means that, given that the coefficients in step 2 are all 1 (diagonal), in step 3 this 1 is divided by the total number of new levels.
Using the same example for type II testing, we also 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} \]
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 to
equate the coefficients in the sub matrix. We use the coefficients for
RACE
, and divide it by the additional levels (there is
interaction between ARMCD
and RACE
, the
additional covariate ARMCD
has 2 levels) and we get 0.5 for
the coefficients.
\[ \begin{matrix} \mu & ARMCD & RACE & RACE & ARMCD:RACE & ARMCD:RACE\\ 0 & 0 & 1 & 0 & 0.5 & 0\\ 0 & 0 & 0 & 1 & 0 & 0.5\\ \end{matrix} \]
Similarly, if we are testing the effect of ARMCD
we will
have the following contrast because RACE
has 3 levels.
\[ \begin{matrix} \mu & ARMCD & RACE & RACE & ARMCD:RACE & ARMCD:RACE\\ 0 & 1 & 0 & 0 & \frac{1}{3} & \frac{1}{3}\\ \end{matrix} \]
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.
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.
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
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} \]
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.
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.