ASR021. Three-way interaction linear models - Treating hypertension

The complete script for this example can be downloaded here:

Dataset

The models that we will fit here are based on the D017 dataset, and the first few rows are presented below:

drug biofeed diet pressure
X Absent No 173
X Absent No 194
X Absent No 197
X Absent No 190
X Absent No 176
X Absent No 198
X Absent Yes 164
X Absent Yes 190


Models

In this example, we will fit two three-way interaction models to analyze the effect of three types of treatments on hypertension. The first model will include all the combinations of the three variables while the second one will only include the interaction term between the three types of treatments. The specifications of the full model is:

\[ y = \mu + drug * biofeed * diet + e\ \]

where,

    \(y\) is the blood pressure,

    \(\mu\) is the overall mean,

    \(drug*biofeed*diet\) corresponds to the following fixed effects:

      \(drug\) is the fixed effect of the medication,

      \(biofeed\) is the fixed effect of psychological feedback,

      \(diet\) is the fixed effect of diet,

      \(drug:biofeed\) is the fixed interaction term between drug and biofeed,

      \(drug:diet\) is the fixed interaction term between drug and diet,

      \(biofeed:diet\) is the fixed interaction term between biofeed and diet,

      \(drug:biofeed:diet\) is the fixed interaction term between drug, biofeed, and diet,

    \(e\) is the random residual effect, with \(e \sim \mathcal{N}(0,\,\sigma^{2}_{e})\).


The specification of the partial model is:

\[ y = \mu + drug:biofeed:diet + e\ \] where,

    \(y\) is the blood pressure,

    \(\mu\) is the overall mean,

    \(drug:biofeed:diet\) is the fixed term defined by the combination (concatenation) between the three treatment factors, and

    \(e\) is the random residual effect, with \(e \sim \mathcal{N}(0,\,\sigma^{2}_{e})\).


Now, let’s take a look at how to write these models with ASReml-R. Note that before fitting them, drug, biofeed, and diet need to be set as factors.

asr021_full <- asreml(
  fixed = pressure ~ drug*biofeed*diet, 
  residual = ~units,
  data = d017
)

Note that drug*biofeed*diet with the use of * will be expanded to its main effects, two-way interactions and three-way interaction.

asr021_partial <- asreml(
  fixed = pressure ~ drug:biofeed:diet, 
  residual = ~units,
  data = d017
)

In this case, there will not be expansion of the term drug:biofeed:diet with the use of :, here a single concatenated factor will be assumed.


Exploring output

Let’s see what we obtain when estimating the predictions for the full model:

preds_full <- predict(asr021_full, classify = 'drug*biofeed*diet')$pvals
preds_full
   drug biofeed diet predicted.value std.error    status
1     X  Absent   No             188  5.109903 Estimable
2     X  Absent  Yes             173  5.109903 Estimable
3     X Present   No             168  5.109903 Estimable
4     X Present  Yes             169  5.109903 Estimable
5     Y  Absent   No             200  5.109903 Estimable
6     Y  Absent  Yes             187  5.109903 Estimable
7     Y Present   No             204  5.109903 Estimable
8     Y Present  Yes             172  5.109903 Estimable
9     Z  Absent   No             209  5.109903 Estimable
10    Z  Absent  Yes             182  5.109903 Estimable
11    Z Present   No             189  5.109903 Estimable
12    Z Present  Yes             173  5.109903 Estimable

And then for the partial model:

preds_partial <- predict(asr021_partial, classify = 'drug:biofeed:diet')$pvals
preds_partial
   drug biofeed diet predicted.value std.error    status
1     X  Absent   No             188  5.109903 Estimable
2     X  Absent  Yes             173  5.109903 Estimable
3     X Present   No             168  5.109903 Estimable
4     X Present  Yes             169  5.109903 Estimable
5     Y  Absent   No             200  5.109903 Estimable
6     Y  Absent  Yes             187  5.109903 Estimable
7     Y Present   No             204  5.109903 Estimable
8     Y Present  Yes             172  5.109903 Estimable
9     Z  Absent   No             209  5.109903 Estimable
10    Z  Absent  Yes             182  5.109903 Estimable
11    Z Present   No             189  5.109903 Estimable
12    Z Present  Yes             173  5.109903 Estimable

We can observe that we get the same predictions with the two different models. This is because in both cases we are requesting predictions that correspond to the combination of the levels of those 3 factors.


However, the big difference will be on the ANOVA tables for the two models:

wald(asr021_full, denDF = 'numeric')$Wald
                  Df denDF     F.inc           Pr
(Intercept)        1    60 1.564e+04 0.000000e+00
drug               2    60 1.173e+01 5.018624e-05
biofeed            1    60 1.307e+01 6.150698e-04
diet               1    60 3.320e+01 3.053311e-07
drug:biofeed       2    60 8.266e-01 4.424566e-01
drug:diet          2    60 2.882e+00 6.381527e-02
biofeed:diet       1    60 2.043e-01 6.529374e-01
drug:biofeed:diet  2    60 3.431e+00 3.883423e-02
wald(asr021_partial, denDF = 'numeric')$Wald
                  Df denDF     F.inc           Pr
(Intercept)        1    60 15640.000 0.000000e+00
drug:biofeed:diet 11    60     7.656 4.644875e-08

With the full model, the ANOVA tests the levels of each of the factors for main effects, two-way and three-way interactions. In contrast, the partial model, the ANOVA table only indicates if there are any significant differences between any of the 12 means (df = 11). Hence, with the full model, we obtain a clearer picture on where the possible interactions are happening.