ASR019. Obtaining predictions for specific values of an explanatory variable - Synthetic weight loss data

The complete script for this example can be downloaded here:

Dataset

The model that we will fit here is based on the D015 synthetic dataset, and the first few rows are presented below:

trt	patient	weight	loss
trtA	1	198	2
trtA	2	207	3
trtA	3	262	8
trtA	4	189	1
trtA	5	230	3
trtA	6	278	5
trtA	7	222	3
trtA	8	240	4

Model

In this example we will fit a linear model with the following specification:

\[ y = \mu + trt + weight + trt:weight + e\ \]

where,

\(y\) is the yield of oats,

\(\mu\) is the overall mean,

\(trt\) is the fixed effect of treatment,

\(weight\) is the fixed effect of the initial weight,

\(trt:weight\) is the fixed effect of the interaction between treatment and initial weight,

\(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 the model with ASReml-R. Note that before fitting the model, trt needs to be set as factor.

asr019 <- asreml(
  fixed = loss ~ trt + weight + trt:weight, 
  residual = ~units,
  data = d015
)

The above model has only fixed effects, and a simple residual term.

Exploring output

Let’s take a look at the ANOVA table:

wald(asr019, denDF = 'numeric', ssType = 'conditional')$Wald

            Df denDF   F.inc  F.con Margin           Pr
(Intercept)  1    18 292.800 27.600        5.381724e-05
trt          2    18  61.920 61.160      A 9.410860e-09
weight       1    18  60.280 60.280      A 3.743320e-07
trt:weight   2    18   6.178  6.178      B 9.061740e-03

Note that in this case, we are using ssType = 'conditional', which allows for testing the significance of each model term given that all other terms are included in the model. This is in contrast to the sequential options that is used by default that tests the terms in a sequential order.

We can calculate the predictions such as:

predict(asr019, classify = 'trt:weight')$pvals

      trt   weight predicted.value std.error    status
1 placebo 235.5833        1.914169 0.5374119 Estimable
2    trtA 235.5833        4.050402 0.5481806 Estimable
3    trtB 235.5833        9.751260 0.5324321 Estimable

Here we can see the predictions of the different treatments based on the average initial weight of 235.6 pounds.

However, if we are interested to estimate predictions based on specific initial weight values, it is possible to add another component to the function predict() as shown in the following code:

predict(asr019, classify = 'trt:weight', levels = list(weight = c(190, 230, 270)))$pvals

      trt weight predicted.value std.error    status
1 placebo    190      0.05409661 0.9505139 Estimable
2 placebo    230      1.68633574 0.5569940 Estimable
3 placebo    270      3.31857486 0.6990672 Estimable
4    trtA    190      1.40614042 0.8806557 Estimable
5    trtA    230      3.72651645 0.5323593 Estimable
6    trtA    270      6.04689249 0.9327052 Estimable
7    trtB    190      4.51815708 0.9063103 Estimable
8    trtB    230      9.11027680 0.5445855 Estimable
9    trtB    270     13.70239653 0.7268101 Estimable

Now, we are able to obtain predictions on how effective the different treatments are based on a given initial weight of a patient.