In harrysouthworth/margarita: Support functions for texmex

Introduction

This document illustrates some usage of the margarita package for R [@R]. The margarita package is designed to simplify extreme value modelling of clinical trial lab safety data, as described by Southworth & Heffernan [@southHeff]. Some functionality that appeared in previous versions of margarita has now been transfered to texmex.

Development of the margarita package was partially funded by AstraZeneca who have kindly made it publicly available.

The major gains in efficiency in the analysis of clinical trial safety data come from simpler simulation of return levels and upper limit of normal exceedance probabilities, as well as functions to support a standard approach to robust regression modelling.

Attempts are made at automatic threshold selection, including use of the extended generalized Pareto 3 (EGP3) distribution of Papastathopoulos & Tawn [-@egp]. However, this functionality relies on the known context of clinical trial safety data and might not work well in other contexts. In any event, automatic threshold selection appears to be a tricky matter [@scarrottMacDonald] and standard diagnostic plots should always be studied.

Previous versions of margarita provided data manipulation functions for working with CDISC SDTM data. The preferred approach is to now learn dyplr [@dplyr].

Data importation

Typically, the data will exist in SAS datasets with file extensions .sas7bdat. In what follows, it is assumed that if a SAS formats catalogue file exists, it is in a standard location (which is often a safe assumption) and is not broken (which isn't).

Supposing the SAS datasets to be held in directory user/data/ (the path being relative to the R workspace in use), we import using the following commands.

library(margarita)
library(haven)
dataPath <- "/user/data/"
dm <- read_sas(file.path(dataPath, "dm.sas7bdat"))
lb <- read_sas(file.path(dataPath, "lb.sas7bdat"))

We can examine the first few rows, and summarize the data as follows:

# Previous codeblock was not evaluated so need to attach here.
library(margarita)
library(haven)
library(knitr) # for kable

head(dm)
summary(dm)

Non-standard SAS datasets

SAS dataset can be troublesome if SAS formats catalogues (particularly broken ones) are used, or if binary compressed SAS datasets are used.

The haven package [@haven] does a good job of importing most SAS datasets, but does not support binary compressed versions. In practice, binary compressed SAS datasets appear to be uncommon. Unfortunately, haven isn't psychic, so if your formats catalogues are broken, you're going to have to put in some manual work. In any sane world, you won't be using formats catalogues anyway.

Finally, it might be the case that expected columns are missing from the data, such as the number of days on study. Whilst, in general, it would be better to construct such columns in accordance with existing standards using SAS, it is not difficult to use the dplyr package [@dplyr] to do most of the calculations.

Data aggregation

Typically, the lab dataset will contain data on a lot of lab variables. In our current case, we are interested only in ALT, so we subset down on that.

library(dplyr)
alt <- filter(lb, lbtestcd == "ALT")

We now reduce the data to the baseline values and post-baseline maxima. First we construct a dataset with baselines only, then construct one with the maxima only. Finally, we join the two.

b.alt <- filter(alt, visit == "Baseline") %>%
  select(usubjid, baseline=aval)

alt <- filter(alt, visit == "Week 6") %>%
  group_by(usubjid) %>%
  summarize(alt = max(aval, na.rm=TRUE)) %>%
  left_join(b.alt, by="usubjid")

If the lb dataset already contains the treatment decodes and population flags, the dm dataset will not likely be required. Supposing, though, that the treatment decodes need to be added to alt we do another left join:

alt <- left_join(alt, select(dm, usubjid, arm), by="usubjid")

Note that it will often be the case that additional steps are necessary to tidy up the analysis dataset, such as removing screen failures and ensuring that 'param' or 'lbtestcd' really does uniquely identify lab tests (e.g. that they don't need to be combined with a column indicating if the source was serum, urine, or something else).

Finally, note that subjects who have no baseline value recorded, or no post-baseline value recorded, cannot be used in the analysis (unless some form of multiple imputation is to be employed, but in the extreme value setting it is not obvious how this would be done). As such, all cases with missing values should be dropped.

We should now have an analysis-ready dataset, the first few rows of which look something like this:

head(alt)

Shiftplots of the data appear in Figure~\ref{fig:shiftplots}.

shiftplot(alt, aes(baseline, alt), by=~arm,
          xlab="Baseline ALT (U/L)",
          ylab="Maximum ALT (U/L)")

Robust regression

Robust regression can be performed using the lmr (linear model robust) function. This function is a wrapper to rlm in the MASS [@MASS] package, setting the method to MM, efficiency for Gaussian data to 85%, and using bisquare weight functions. It also computes the coefficient covariance matrix and attaches it to the returned object. A ggplot [@ggplot2] function for diagnostics is provided by texmex.

The use of MM-estimation with the particular settings as described above is used as default as recommended by Maronna et al [-@maronnaMartinYohai] and as described by Southworth & Heffernan [-@southHeff].

In our data, dose is ordered such that Dose $A < B < C < D$, and they are equally spaced on the log scale. We take advantage of this structure by assuming linearity in dose and recoding A -- D as 1 -- 4.

alt <- mutate(alt, ndose = as.numeric(arm)) %>%
  as.data.frame()

mm <- lmr(log(alt) ~ log(baseline) + ndose, data=alt)
ggplot(mm)

The reason for using robust regression is that we do not expect the data to be normally distributed, but to have outliers. Whilst most robust regression methods implicitly assume the distribution to be symmetric, in practice parameter estimates are very little affected by skew in the data. In the figure above, it can be seen in the QQ-plot (top left panel) that the residuals drift away from the reference line to the bottom left, but more so to the top right. This indicates outliers to the right, confirming the lack of normality in the data. The outliers can also be seen in the other panels of the above figure. The bottom left panel shows plots of residuals versus fitted values and the smoother should be near-horizontal (excusing end-effects) through the bulk of the data. The bottom right panel shows the observed values versus the fitted values and should show (as it does) a near linear relationship with little evidence of non-constant variance.

We can also look at boxplots of scaled residuals by treatment group.

boxplot(mm, by="arm")

From the boxplots, we can see a fairly clear dose-response relationship with there being few outliers in the lower dose groups, but more, and larger, outliers in the higher dose groups. In practice, it is often at this stage of the analysis that the presence or absence of treatment effects begins to become clear. Note that the residuals in the boxplots have been scaled so that values beyond about $\pm 2.5$ are inconsistent with normality, justifying the use of robust estimation.

Extreme value modelling

We now move on to fit generalized Pareto (GP) distributions to the residuals from the robust regression.

Threshold selection

Before fitting GP models, we need to select a threshold above which GP models are likely to be reasonable. We proceed firstly by producing some plots to aid threshold selection.

# Do threshold selection plots
alt$r <- resid(mm)
ggplot(gpdThresh(alt$r))

The top row of the above figure displays plots of the scale and shape parameters in the GP distribution above increasing thresholds. Following Coles [-@coles], these should be constant above a threshold suitable for GP modelling and, to make use of as much data as possible, we want to choose the lowest suitable threshold. Accounting for the uncertainty represented by the approximate pointwise confidence intervals on the graphs, a threshold close to zero, or perhaps a little higher, looks like it should be ok.

The lower left panel of the threshold selection plots shows the mean residual life (MRL) plot, again discussed by Coles [-@coles]. The MRL plot should be linear above a suitable threshold, uncertainty accounted for. In practice, interpretation of such plots requires some practice, but in this case a threshold close to zero seems like it might be reasonable.

The bottom right panel of the threshold selection plots contains a threshold selection plot suggested by Papastathoploulos and Tawn [-@egp]. The horizontal axis is a range of thresholds and the vertical axis is the value of estimated power parameter κ in the EGP3 distribution. At the point where κ = 1, the EPG3 distribution is identical to the GP distribution, suggesting that the lowest threshold for which the confidence interval for $\hat\kappa$ contains 1 as a candidate. It can be seen that the displayed pointwise confidence interval crosses or comes close to the horizontal reference line at 1 at thresholds around zero, 0.125 and then above about 0.3. Again, this suggests that a threshold of zero or a little higher might be appropriate.

From the threshold selection plots, the data are all of the residuals from the robust regression with no account being taken of dose (or any other covariate). As such, they can only be considered to be approximations of what we really want, which is a threshold above which GP modelling accounting for dose is appropriate.

A test for threshold selection

Papastathopoulos and Tawn [-@egp] suggest selecting the lowest threshold above which the confidence interval for $\hat\kappa$ contains 1. In the context of clinical trials, we typically have fewer observations than in other settings, so some care is needed in applying this approach.

Because we are modelling residuals from a robust linear model, the shape of the distribution of those residuals should be not far from symmetric and should be centered close to zero. This suggests that any threshold beneath the median of the residuals would be inappropriate because the shape above a lower threshold would not be compatible with the shape of the GP distribution.

Also, due to the sample sizes we expect in practice, if we reject all the residuals beneath about the 75^th^ percentile, then we will be left with too few to do any reasonable modelling with.

As such, when testing for a threshold, we restrict attention to values between the median and the 75^th^ percentile of the residuals (though these default values can be changed).

Further details of the implementation of the EGP3 test in margarita are as follows. The function egp3Thresh calls on egp3RangeFit from the texmex package. It then interpolates through the values of $\hat\kappa$ and its confidence intervals before looking for the lowest threshold at which the confidence interval contains 1, and returning that threshold as well as the threshold for which $\hat\kappa$ is closest to 1. If no confidence interval contains 1, the function issues a warning.

We now try the test on the residual ALT data. As with the threshold selection plots, the test can be considered only approximate because it takes no account of dose, which is what we are actually interested in.

egp3Thresh(alt$r)

None of the confidence intervals contained 1, so we got a warning and a value of Inf. The threshold for which the value of $\hat\kappa$ is closest to 1 is 0.0967. This value corresponds with what we can read off the lower right panel of the threshold selection plots. Approximately r round(mean(alt$r < .0967)*100)% of the residuals are beneath this value, suggesting a slightly lower threshold than the 70^th^ percentile used by Southworth & Heffernan [-@southHeff].

Fitting generalized Pareto distributions

Following the discussion in the previous section, we proceed to fit GP models to the residual ALT values by maximum likelihood and look at diagnostic plots. We threshold at the 60^th^ percentile.

mlmod <- evm(r, data=alt, qu=.6)
ggplot(mlmod)

In the diagnostic plots above, no account has been taken of dose, so these plots can only be considered as approximate guidance on whether GP models are likely to fit well above the chosen threshold. The pointwise simulated envelopes give no cause for concern. In general, the QQ-plot (the top right panel) is usually the one in which any lack of model fit becomes most clear.

Following Southworth & Heffernan [-@southHeff], we wish to account for dose. As with the robust regression modelling previously, we assume linearity in dose. Either one, both or none of the parameters in the GP distribution might be dependent on dose, so we fit a variety of models and compare them using the Akaike Information Criterion (AIC). The AIC provides guidance on model selection, but not a hard rule, so in some cases, we might want to consider the change in deviance due to the addition of a single parameter more carefully. The following table is reproduced from Kass & Raftery [-@kassRaftery], though versions of it have been suggested by many other authors ([@jeffreys, @raftery, @lindsey, @davison, @royall]).

evidence(kable=TRUE)

gpd0 <- evm(r, alt, qu=.6) # null model
gpd1 <- evm(r, alt, qu=.6, phi=~ndose)
gpd2 <- evm(r, alt, qu=.6, xi=~ndose)
gpd3 <- evm(r, alt, qu=.6, phi=~ndose, xi=~ndose)
AICtable(list(gpd0, gpd1, gpd2, gpd3), kable=TRUE)

The function AICtable provides a neat way of producing \LaTeX output containing the AICs, log-likelihoods and related information in a written report (like this one). Alternatively, examine any of the fitted models by printing them:

gpd2

From the table, the model with a linear term for dose in $\phi$ and $\xi$ is preferred by AIC. However, the models with dose in just $\phi$ or just $\xi$ are close contenders. Playing around with different threshods, and considering the fact that the two parameters are not independeing, we eventually choose the model with dose in $\xi$. Using either of the other models does not lead to predictions that suggest wildly different conclusions. We check the diagnostic plots for the MLE fit before refitting using Markov chain Monte Carlo (MCMC).

ggplot(gpd2)

The diagnostic plots give no cause for concern.

In this example, there is a linear dose-response relationship, with $\hat\xi$ increasing as dose increases. Such a relationship makes logical sense. In some examples, it might be the case that a decreasing or non-monotone dose-response is suggested. Such a relationship is unlikely to make sense and the analyst then needs to weigh the strength of evidence in the data (i.e considerations of AIC and deviance) with their knowledge of the data generating mechanism (i.e. the treatment regimens). Statistical inference is only one part of the larger process of scientific data modelling.

gmod <- evm(r, alt, qu=.6, xi=~ndose, method="sim")

Note that fitting using MCMC results (by default) in the function telling the user the acceptance rate. For fairly simple models with few parameters, as we have here, an acceptance rate of a little over 30% is common. If the rate is much higher or lower than this, then something might have gone wrong. Also, if more parameters are included in the model, the acceptance rate should be expected to go down. Section 12.2 of Gelman et al [-@bda] contains some discussion of acceptance rates and efficient MCMC algorithms.

Having refit, we now look at diagnostic plots for the Markov chains:

ggplot(gmod)

The top row of the above figure shows kernel density estimates of the posterior distributions as simulated by the Markov chains. The second row shows the simulated values at each step of the chain and the cumulative means. If the Markov chains have converged on their target distributions, the cumulative means should converge rapidly on stable values and the trace of the chains should look like 'fat hairy caterpillars'. The final row of the above figure shows the autocorrelation functions (ACFs) of the chains. If the algorithm has worked well, the ACFs should rapidly decay to zero.

In general, it is good practice to discard the first part of the Markov chains (the burn-in period) and to thin the chains by discarding values at regular intervals. For example, we could keep 1 in every 4 values from the chains. These steps improve the chances that the chains have converged on the target distribution, reduce autocorrelation, and make the simulated values behave more like a random sample from the target distribution. The help file for evm describes the thin and burn arguments that control these settings.

margarita

The main efficiency that the margarita package introduces to the analysis of clinical trial safety data is to provide built-in functions to aid simulation of quantities of interest such as return levels and probabilities of exceeding high values. This is achieved by first creating an object that contains the robust linear model, the MCMC simulated GP model, and other information required for the calculations.

nd <- data.frame(ndose=1:4)
mods <- margarita(mm, gmod, newdata=nd,
                  trans=log, invtrans=exp,
                  baseline="baseline")

In the call to margarita, we passed in the robust regression model, the MCMC extreme value model, and a data.frame as the newdata argument. The resulting margarita object is then interpreted by simulate functions and used to generate predicted values of interest.

The simulate functions infer from the number of draws from the posterior distributions contained in the extreme value model how many cases to simulate.

If there are no covariates in the model, simulate doesn't need the newdata argument. However, if there are covariates in either the robust linear model, the GP model or both, in order to deal with all kinds of trial designs, the simulate function requires newdata. The rows of newdata should be unique and it should not contain the baseline data. Baseline data are simulated internally by resampling from the data attached to the linear model object. The rows of newdata therefore represent combinations of covariates in which the user is interestd in.

The trans argument tells the simulation function how the response variable is transformed in the robust linear model, and the invtrans argument gives the reverse transformation. These default to log and exp, but in some cases might both be I for the identity (i.e. no) transformation, or to sqrt and function(x) x*x for the square root and square transforms. Note that trans and invtrans are functions, not character strings.

Return levels

Return levels for clinical trial data can be simulated using the simulate function with type='rl' (the default). The simulations for various return levels can be run in a single call by passing a vector argument into the function.

Note that it is common for clinical laboratory safety data to be interpreted in terms of how far the observed values are from being `normal'. As such, each lab variable tends to be issued with an upper limit of normal (ULN) and a lower limit of normal (LLN), though we are usually concerned only with one or the other.

There is a great deal of arbitrariness in the computation of the LLN and ULN. The values ought to represent a low and a high quantile of the distribution in a healthy population, but they vary from lab to lab, partly due to the machines that make the readings differing from lab to lab, partly because the local populations differ from lab to lab, and because the way that low and high quantiles are computed differs. M'Kada et al [-@mkada] describe one lab that calculated the ULN for ALT as the mean plus two standard deviations after removing the 5% most extreme values, but later recalculated it as the mean plus a single standard deviation after removing the 5% most extreme values. Given that ALT is skewed and subject to outliers, it is difficult to justify either approach.

To avoid the risk of increasing variation by using local lab normal ranges, it is recommended to use fixed ranges as published in the public domain. A little searching reveals estimated ULNs for ALT ranging from 26 U/L to 72 U/L. The Merck manual [@merck] suggests 35 U/L, and it is recommended that the Merck manual be the first place to look when attempting to decide on a reasonable reference range. Information on the ULN and LLN for a subset of the variables listed in the Merck manual appears in normalRanges for convenience. The ULN and LLN values are given in SI units and in 'conventional' units.

names(normalRanges)
normalRanges[, c("test", "uln.conventional", "uln.si")]

Finally, in clinical trials we observe patients at baseline (i.e. before treatment) and then after a period of treatment. Baseline values are almost always correlated with on-treatment values, so a more appropriate approach to interpreting the effect of drug on lab values would be to eliminate the baseline effect and study what remains. An approach that comes close to this is to look at changes from baseline. In practice, we work with the residuals from a robustly fit linear model, but transform results back to the scale of the change from baseline and multiples of ULN.

ULN <- 35

rl1000 <- simulate(mods, type="rl", M=1000)

s <- summary(rl1000)
suln <- s / ULN
sFold <- summary(rl1000, scale="proportional")
sDiff <- summary(rl1000, scale="difference")

# s and the others are lists with an element for each element of M
rownames(s[[1]]) <- paste("Dose", LETTERS[1:4])
rownames(sFold[[1]]) <- rownames(sDiff[[1]]) <-
  rownames(suln[[1]]) <- rownames(s[[1]])
gs <- ggplot(s, xlab="U/L")
gsuln <- ggplot(suln, xlab="Multiples of ULN")
gsFold <- ggplot(sFold, xlab="Fold-change from baseline")
gsDiff <- ggplot(sDiff, xlab="Change in U/L from baseline")
grid.arrange(gs, gsuln, gsFold, gsDiff)

Notice that for return levels, because of the way they are calculated, it is straightforward to change the scale after the calculations have been done. It is achieved by providing optional arguments to summary or by manipulating the object returned by summary directly.

Here we have included the absolute change from baseline in the output, though for ALT it is usually the proportional change that is of more interest.

Threshold exceedance probabilities

When simulating threshold exceedance probabilities, the scale of the predicted values needs to be specified in the call to simulate rather than summary.

First we find probabilities of threshold exceedance on the scale of the raw data (i.e. multiples of ULN), then in terms of multiples of baseline.

pULNs <- simulate(mods, type="prob",
                  M=ULN * c(1, 3, 10, 20),
                  Mlabels=c("P(ALT > ULN)", "P(ALT > 3xULN)",
                            "P(ALT > 10xULN)", "P(ALT > 20xULN)"))
ps <- summary(pULNs)

# Manually set the treatment group names
names(ps) <- LETTERS[1:4]
ggplot(ps, ncol=1, pointest="median")

# Simulate fold-changes from baseline
pBase <- simulate(mods, type="prob", M=c(2, 5, 10, 20),
                  Mlabels=c("2-fold", "5-fold",
                            "10-fold", "20-fold"),
                  scale="proportional")
pbs <- summary(pBase)

# Manually change treatment group names
names(pbs) <- LETTERS[1:4]
ggplot(pbs, ncol=1,
       xlab="Proportion of patients with exceedance")

Note that in the preceding Figures, the scales of the horizontal axes vary. The reason is that we are primarily interested in the large increases from baseline, but if we use the same scale for these as for the smaller increases, the useful information is compressed into a tidy part of the graph, making it difficult to read.

Note some of the complexity caused by the weird geometry of having a fixed threshold (multiples of the ULN) against which to compare changes from baseline. It is possible for any particular patient to have a baseline value and regression coefficients that cause his expected on-treatment value to be above threshold M. So we need to simulate exceedance proabilities in 2 steps.

We are interested in P(y = r + u > M | u > M). If the simulated patient is in the group whose residuals fall above the threshold, the probability is equal to the probability of exceeding the modelling threshold in the original GP model. Otherwise, the probability depends on the distribution of residuals beneath the modelling threshold. We thus have a mixture distribution:

$f(y) = \theta g(y) + (1 - \theta) h(y)$

in which $\theta$ is our modelling threshold, $g(y)$ is the GP distribution and $h(y) = P(y > M | u > M)$.

$g(y)$ is what most of the focus of this article has been on. To account for the second component of the mixture distribution, we resample residuals below the threshold and count how many of them breach the return level M.

Note that the distribution of probabilities of threshold exceedance will usually be highly skewed with a mode at 0. As such, the expected value (mean) will often be greater than the median and can fall outside of the interval estimates (especially the 50% interval estimate). To avoid confusion, the point estimate that is computed by default is the median, which will underestimate the expected proportion.

lapply(unclass(pULNs)[1:4], function(x) apply(x, 2, summary))

References

harrysouthworth/margarita documentation built on Aug. 19, 2021, 5 a.m.