Bayesian survival analysis with BUGS

AFT models are the survival counterpart of linear models. Survival time $T$ in logarithmic scale is expressed in terms of a linear combination of covariates ${\bm{x}}$ with regression coefficients $\bm{\beta}$ and a measurement error $\xi$ as follows:

The Cox proportional hazards model (Cox 1972) expresses the hazard function $h(t)$ of the survival time of each individual of the target population as the product of a common baseline hazard function $h_{0}(t)$, which determines the shape of $h(t)$, and an exponential term which includes the relevant covariates ${\bm{x}}$ with regression coefficients $\bm{\beta}$ as follows:

The estimation of the regression coefficients in the Cox model under the frequentist approach can be obtained without specifying a model for the baseline hazard function by using partial likelihood methodology (Cox 1972). This is not the case of Bayesian analysis which in general needs to specify some model for the baseline hazard (Christensen et al. 2011). Depending on the context of the study, baseline hazard misspecification can imply a loss of valuable model information that makes it impossible to fully report the estimation of the outcomes of interest, such as probabilities or survival curves for relevant covariate patterns (Royston 2011). This is specially important in survival studies where $h_{0}(t)$ represents the natural course of a disease or an infection, or even the control group when comparing several treatments.

Baseline hazard functions can be defined through parametric or semi-parametric approaches. Parametric models give restricted shapes which do not allow the presence of irregular behaviour. Semi-parametric choices result in more flexible baseline hazard shapes that allow for multimodal patterns such as piecewise constant functions or spline functions. Theoretical and methodological aspects of more flexible approaches to define the baseline hazard function can be found in several specific references such as Ibrahim et al. 2001, Lin et al. 2015, Mitra and Müller 2015, Bogaerts et al. 2017, and Lázaro et al. 2020.

Cure models deal with target populations in which a part of the individuals cannot experience the event of interest. This type of models has widely been matured as a consequence of the discovery and development of new treatments against cancer. The rationale of considering a cure subpopulation comes from the idea that a successful treatment removes totally the original tumor and the individual cannot experience any recurrence of the disease. These models allow to estimate the probability of cure, a key and valuable outcome in cancer research.

The effect of a baseline covariate vector ${\bm{x}}$ on the cure fraction $\eta$ for each individual is typically modelled by means of a logistic link function, logit$(\eta)$, but the probit link or the complementary log-log link could also be used. Covariates for modelling $T$ in the uncured subpopulation are usually considered via Cox models. Cure fraction model is usually known as the incidence model and the survival model (i.e., time-to-event $T$ in the uncured subpopulation) as the latency model Klein et al. 2013.

We consider a bone marrow transplant dataset, referred to as bmt. It is available from the smcure package (Cai et al. 2012):

R> library("smcure")
R> data("bmt")
R> str(bmt)
’data.frame’:Ψ91 obs. of  3 variables:
 $ Time  : num  11 14 23 31 32 35 51 59 62 78 ...
 $ Status: num  1 1 1 1 1 1 1 1 1 1 ...
 $ TRT   : num  0 0 0 0 0 0 0 0 0 0 ...

This dataset refers to a bone marrow transplant study for the refractory acute lymphoblastic leukemia patients, in which 91 patients were divided into two treatment groups (Kersey et al. 1987). The following variables were observed for each patient:

• •

  Time: time to death (in days).

• •

  Status: censoring indicator (1: if patient is uncensored; 0: otherwise).

• •

  TRT: treatment group indicator (0: allogeneic; 1: autologous).

The (cure subpopulation) incidence model for each individual is expressed in terms of a logistic regression:

Survival time for each individual in the uncured subpopulation is modelled from a proportional hazards specification:

We assume prior independence and specify prior marginal distributions based on non-informative distributions commonly employed in the literature. The $\beta$’s follow a N($0,0.001$), while $\lambda$ and $\alpha$ follow a Gamma($0.01,0.01$) and a Un($0,10$), respectively.

The likelihood function for mixture cure models is not implemented in JAGS, so the “zeros trick” approach using a Poisson distribution has also been used to specify it indirectly (Ntzoufras 2009; Lunn et al. 2012).

We have created two design matrices, one for model (10) and another for model (11), with the TRT covariate:

R> XC <- model.matrix(~ TRT, data = bmt) # Reference = allogeneic
R> XU <- model.matrix(~ TRT, data = bmt)
R> XU <- matrix(XU[,-1], ncol = 1) # Remove intercept

Listing 3 shows a generic implementation of a mixture cure model in BUGS syntax using bmt data.

Once the variables have been defined, a list with all the elements required in the model is created:

R> d.jags <- list(n = nrow(bmt), t = bmt$Time, XC = XC, XU = XU,
+    delta = bmt$Status, zeros = rep(0, nrow(bmt)), NbetasC = ncol(XC), NbetasU = ncol(XU))

The initial values for each mixture cure model parameter are passed to JAGS using a function that returns a list of random values:

R> i.jags <- function(){
+    list(betaC = rnorm(ncol(XC)), betaU = rnorm(ncol(XU)), lambda = runif(1), alpha = runif(1))
+  }

The vector of monitored/saved parameters is:

R> p.jags <- c("betaC", "betaU", "alpha", "lambda")

Next, the JAGS model is compiled:

R> library("rjags")
R> m3 <- jags.model(data = d.jags, file = "Cure.txt", inits = i.jags, n.chains = 3)

We now run the model for 10000 burn-in simulations:

R> update(m3, 10000)

Finally, the model is run for 100000 additional simulations to keep one in 100 so that a proper thinning is done:

R> res <- coda.samples(m3, variable.names = p.jags, n.iter = 100000, n.thin = 100)

Similarly to the first example (Section 2.1), numerical and graphical summaries of the model parameters can be obtained using the summary and densplot functions, respectively. Gelman and Rubin’s convergence diagnostic can be calculated with the gelman.diag function, and the traceplot function provides a visual way to inspect sampling behaviour and assesses mixing across chains and convergence.

Next, simulations from the three Markov chains are merged together for inference:

R> result <- as.mcmc(do.call(rbind, res))

The posterior samples of each parameter are obtained by:

R> alpha <- result[,1]; betaC1 <- result[,2]; betaC2 <- result[,3]
R> betaU <- result[,4]; lambda <- result[,5]

Table 2 shows posterior summaries for the mixture cure model parameters using bmt data.

As we discussed earlier, the cure fraction ($\eta$) is a relevant quantity for mixture cure models. For allogeneic (TRT=0) or autologous (TRT=1) treated patients, it is modelled as:

We can easily summarise the posterior distribution of the cure fraction for individuals in both groups:

R> CP.allo <- exp(betaC1) / (1 + exp(betaC1))
R> CP.auto <- exp(betaC1 + betaC2) / (1 + exp(betaC1 + betaC2))
R> summary(cbind(CP.allo, CP.auto))
    CP.allo           CP.auto
 Min.   :0.01604   Min.   :0.02277
 1st Qu.:0.22426   1st Qu.:0.15680
 Median :0.26850   Median :0.19477
 Mean   :0.27146   Mean   :0.19925
 3rd Qu.:0.31515   3rd Qu.:0.23699
 Max.   :0.65573   Max.   :0.55898

The uncured survival curve based on posterior samples is another relevant information in this type of studies. So, from the posterior samples obtained above, we can summarise the posterior distribution of the mean value of the uncured survival curve for allogeneic and autologous treated patients in a grid of points as follows:

R> grid <- 100
R> time <- seq(0, bmt$Time, len = grid)
R> surv.allo <- surv.auto <- vector()
R> for(l in 1:grid){
+     surv.allo[l] <- mean(exp(-lambda * time[l]^alpha))
+     surv.auto[l] <- mean(exp(-lambda * exp(betaU) * time[l]^alpha))
+  }

Figure 3 shows the difference between both curves using the code below:

R> library("ggplot2")
R> treat.col <- rep(0:1, each = grid)
R> treat.col[treat.col == 0] <- "allogeneic"
R> treat.col[treat.col == 1] <- "autologous"
R> df <- data.frame(time = rep(time, 2), survival = c(surv.allo, surv.auto),
+    treatment = treat.col)
R> ggplot(data = df, aes(x = time, y = survival, group = treatment, colour = treatment)) +
+    geom_line() + theme_bw() + theme(legend.position = "top")

We consider a stem-cell transplanted patients dataset, referred to as okiss. It is available from the compeir package (Grambauer and Neudecker 2011):

R> library("compeir")
R> data("okiss")
R> str(okiss)
’data.frame’:Ψ1000 obs. of  4 variables:
 $ time  : ’times’ num  21 6 14 12 11 18 10 8 17 10 ...
  ..- attr(*, "format")= chr "h:m:s"
 $ status: num  2 1 2 2 2 2 7 2 2 2 ...
 $ allo  : num  1 1 1 1 1 1 1 0 1 1 ...
 $ sex   : Factor w/ 2 levels "f","m": 1 1 1 2 2 2 2 1 2 2 ...

This dataset provides information about a sub-sample of 1000 patients enrolled in the ONKO-KISS programme, which is part of the German National Reference Centre for Surveillance of Hospital-Acquired Infections (Dettenkofer et al. 2005). These patients have been treated by peripheral blood stem-cell transplantation, which has become a successful therapy for severe hematologic diseases. After transplantation, patients are neutropenic, i.e., they have a low count of white blood cells, which are the cells that primarily avert infections (Beyersmann et al. 2007). Occurrence of bloodstream infection during neutropenia is a severe complication. The following variables were observed for each patient:

• •

  time: time (in days) of neutropenia until first event.

• •

  status: event status indicator (1: infection; 2: end of neutropenia; 7: death; 11: censored observation).

• •

  allo: transplant type indicator (0: autologous; 1: allogeneic).

• •

  sex: sex of each patient (m: if patient is male; f: if patient is female).

We have redefined the status variable using an auxiliary one (delta) in a matrix format:

R> delta <- matrix(c(as.integer(okiss$status == 1), as.integer(okiss$status == 2),
+    as.integer(okiss$status == 7)), ncol = 3)
R> head(delta)
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    0
[3,]    0    1    0
[4,]    0    1    0
[5,]    0    1    0
[6,]    0    1    0

Cause-specific hazard functions for infection ($k=1$), end of neutropenia ($k=2$), and death ($k=3$) are modelled from a proportional hazard specification:

We have created a design matrix X with the covariates allo and sex:

R> X <- model.matrix(~ allo + sex, data = okiss) # Reference = female
R> X <- X[,-1] # Remove intercept

Listing 4 shows a generic implementation of a competing risks model in BUGS syntax using okiss data.

Once the variables have been defined, a list with all the elements required in the model is created:

R> d.jags <- list(n = nrow(okiss), t = as.vector(okiss$time), X = X,
+    delta = delta, zeros = rep(0, nrow(okiss)), Nbetas = ncol(X), Nrisks = ncol(delta))

The initial values for each competing risks model parameter are passed to JAGS using a function that returns a list of random values:

R> i.jags <- function(){
+    list(beta = matrix(rnorm(ncol(X) * ncol(delta)), ncol = ncol(delta)),
+      lambda = runif(ncol(delta)), alpha = runif(ncol(delta)))
+  }

The vector of monitored/saved parameters is:

R> p.jags <- c("beta", "alpha", "lambda")

Next, the JAGS model is compiled:

R> library("rjags")
R> m4 <- jags.model(data = d.jags, file = "CR.txt", inits = i.jags, n.chains = 3)

We now run the model for 1000 burn-in simulations:

R> update(m4, 1000)

Finally, the model is run for 10000 additional simulations to keep one in 10 so that a proper thinning is done:

R> res <- coda.samples(m4, variable.names = p.jags, n.iter = 10000, n.thin = 10)

Similarly to the first example (Section 2.1), numerical and graphical summaries of the model parameters can be obtained using the summary and densplot functions, respectively. Gelman and Rubin’s convergence diagnostic can be calculated with the gelman.diag function, and the traceplot function provides a visual way to inspect sampling behaviour and assesses mixing across chains and convergence.

Next, simulations from the three Markov chains are merged together for inference:

R> result <- as.mcmc(do.call(rbind, res))

The posterior samples of each parameter are obtained by:

R> alpha1 <- result[,1]; alpha2 <- result[,2]; alpha3 <- result[,3]; beta11 <- result[,4]
R> beta21 <- result[,5]; beta12 <- result[,6]; beta22 <- result[,7]; beta13 <- result[,8]
R> beta23 <- result[,9]; lambda1 <- result[,10]; lambda2 <- result[,11]; lambda3 <- result[,12]

Table 3 shows posterior summaries for the competing risks model parameters using okiss data.

As discussed in Section 4, the cumulative incidence function $F_{k}(t)$ in ( 13) is the most appropriate way to analyse the evolution of each cause $k$ over time. For our proportional hazard specification (14), it is given by:

The integral in (15) has no closed form, so some approximate method of integration is required. To do this, we first have created a function fk which describes the integrand of (15):

R> fk <- function(u.vect, lambda, alpha, beta, x, k){
+  res <- sapply(u.vect, function(u){
+    # Cause-specific hazard
+    hk <- lambda[k] * alpha[k] * (u^(alpha[k] - 1)) * exp(sum(unlist(beta[,k]) * x))
+    # Cumulative cause-specific hazard
+    Hk <- lambda * (rep(u, length(lambda))^alpha) * exp((t(beta) %*% matrix(x, ncol = 1))[,1])
+    # Cause-specific hazard x Overall survival
+    aux <- hk * exp(-sum(Hk))
+    return(aux)
+  })
+  return(res) }

Next, we have created a function cif that computes $F_{k}(t)$ in (15) by integrating out the fk function using the integral function available from the pracma package (Borchers 2019).

R> library("pracma")
R> cif <- function(tt, lambda, alpha, beta, x, k){
+    return(integral(fk, xmin = 0, xmax = tt, method = "Simpson", lambda = lambda, alpha = alpha,
+      beta = beta, x = x, k = k)) }

Finally, we have constructed a function mcmc$\_$cif that takes the output from JAGS (variable obj) and computes $F_{k}(t)$ in (15) for a vector of times (variable t.pred) using covariates x. Note that mcmc$\_$cif is based on the mclapply function, available from the parallel package to speed computations up (R Core Team 2020).

R> library("parallel")
R> options(mc.cores = detectCores())
R> mcmc_cif <- function(obj, t.pred, x){
+    var.names <- names(obj)
+    # Indices of beta’s, alpha’s, and lambda’s
+    b.idx <- which(substr(var.names, 1, 4) == "beta")
+    a.idx <- which(substr(var.names, 1, 5) == "alpha")
+    l.idx <- which(substr(var.names, 1, 6) == "lambda")
+    # Number of causes and number of covariates
+    K <- length(a.idx)
+    n.b <- length(b.idx) / K
+    # Sub-sample to speed up computations
+    samples.idx <- sample(1:nrow(obj), 200)
+
+    res <- lapply(1:K, function(k){
+      sapply(t.pred, function(tt){
+        aux <- mclapply(samples.idx, function(i){
+          cif(tt, alpha = unlist(obj[i, a.idx]), lambda = unlist(obj[i, l.idx]),
+             beta = matrix(unlist(c(res[i, b.idx])), nrow = n.b), x = x, k = k)
+        })
+        return(mean(unlist(aux)))
+      })
+    })
+    return(res)
+  }

Hence, we redefine the MCMC output as a data.frame and set a vector of times to evaluate the cumulative incidence function. In this example, we are interested in this function when both covariates are 1 (i.e., allogeneic transplant and male).

R> res <- as.data.frame(result)
R> t.pred <- seq(0, 100, by = 2.5)
R> cum_inc <- mcmc_cif(res, t.pred, c(1, 1))

Figure 4, generated with the code below, shows the posterior mean of the cumulative incidence function for a man with an allogeneic transplant for the three types of events considered in the okiss data.

R> library("ggplot2")
R> df <- data.frame(cif = unlist(cum_inc), time = t.pred,
+    cause = rep(c("infection", "end of neutropenia", "death"), each = length(t.pred)))
R> ggplot(data = df, aes(x = time, y = cif, group = cause)) + geom_line(aes(color = cause)) +
+    ylab("cumulative incidence") + ylim(c(0,1)) + theme_bw() + theme(legend.position = "top")

We consider a heart transplant dataset, referred to as heart2. It is available from the p3state.msm package (Meira-Machado and Roca-Pardinas 2012):

R> library("p3state.msm")
R> data("heart2")
R> str(heart2)
’data.frame’:Ψ103 obs. of  8 variables:
 $ times1 : num  50 6 1 36 18 3 51 40 85 12 ...
 $ delta  : int  0 0 1 1 0 0 1 0 0 1 ...
 $ times2 : num  0 0 15 3 0 0 624 0 0 46 ...
 $ time   : int  50 6 16 39 18 3 675 40 85 58 ...
 $ status : int  1 1 1 1 1 1 1 1 1 1 ...
 $ age    : num  -17.16 3.84 6.3 -7.74 -27.21 ...
 $ year   : num  0.123 0.255 0.266 0.49 0.608 ...
 $ surgery: int  0 0 0 0 0 0 0 0 0 0 ...

This dataset provides information about a sample of 103 patients of the Stanford Heart Transplant Program (Crowley and Hu 1977). The patients are initially on the waiting list (state 1) and can either be transplanted (state 2, non-terminal event) and then die (state 3, terminal event), or just one or none of them because they continue to be on the waiting list. The following variables were observed for each patient:

• •

  times1: time of transplant/censoring time (state 2).

• •

  delta: transplant indicator (1: yes; 0: no).

• •

  times2: time to death since the transplant/censoring time (state 3).

• •

  time: times1 + times2.

• •

  status: censoring indicator (1: dead; 0: alive).

• •

  age: age - 48 years.

• •

  year: year of acceptance (in years after 1 Nov 1967).

• •

  surgery: prior bypass surgery (1: yes; 0: no).

The patients had the following characteristics: 4 were censored for both events (transplant and death), 24 moved from state 1 (waiting list) to state 2 (transplant) and survived, 30 moved from state 1 to state 3 (death) without going through state 2, and 45 moved from state 1 to state 2 and then to state 3.

We have redefined delta and status variables using an auxiliary one (event) in a matrix format:

R> event <- matrix(c(heart2$delta, heart2$status * (1 - heart2$delta),
+    heart2$delta * heart2$status), ncol = 3)
R> head(event)
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    0    1    0
[3,]    1    0    1
[4,]    1    0    1
[5,]    0    1    0
[6,]    0    1    0

Hazard functions for survival times $T_{12}$, $T_{13}$ and $T_{23}$ are modelled using a proportional hazard specification:

We have created a design matrix X with the covariates age, year and surgery, and defined a time3 $=t-t_{12}$ variable according to semi-Markovian specification:

R> X <- model.matrix(~ age + year + surgery, data = heart2)
R> X <- X[,-1] # Remove intercept
R> time3 <- heart2$times2
R> time3[time3 == 0] <- 0.0001

The values of time3 equal to zero have been replaced by 0.0001 to avoid computational problems when calculating $(t-t_{12})^{\alpha_{3}-1}$.

Listing 5 shows a generic implementation of an illness-death model in BUGS syntax using heart2 data.

Once the variables have been defined, a list with all the elements required in the model is created:

R> d.jags <- list(n = nrow(heart2), t1 = heart2$times1, t2 = heart2$time, t3 = time3,
+    X = X, event = event, zeros = rep(0, nrow(heart2)), Nbetas = ncol(X))

The initial values for each illness-death model parameter are passed to JAGS using a function that returns a list of random values:

R> i.jags <- function(){
+    list(beta = matrix(rnorm(3 * ncol(X)), ncol = 3), lambda = runif(3), alpha = runif(3))
+  }

The vector of monitored/saved parameters is:

R> p.jags <- c("beta", "alpha", "lambda")

Next, the JAGS model is compiled:

R> library("rjags")
R> m5 <- jags.model(data = d.jags, file = "IllDeath.txt", inits = i.jags, n.chains = 3)

We now run the model for 1000 burn-in simulations:

R> update(m5, 1000)

Finally, the model is run for 10000 additional simulations to keep one in 10 so that a proper thinning is done:

R> res <- coda.samples(m5, variable.names = p.jags, n.iter = 10000, n.thin = 10)

Similarly to the first example (Section 2.1), numerical and graphical summaries of the model parameters can be obtained using the summary and densplot functions, respectively. Gelman and Rubin’s convergence diagnostic can be calculated with the gelman.diag function, and the traceplot function provides a visual way to inspect sampling behaviour and assesses mixing across chains and convergence.

Next, simulations from the three Markov chains are merged together for inference:

R> result <- as.mcmc(do.call(rbind, res))

The posterior samples of each parameter are obtained by:

R> alpha1 <- result[,1]; alpha2 <- result[,2]; alpha3 <- result[,3]
R> beta11 <- result[,4]; beta21 <- result[,5]; beta31 <- result[,6]
R> beta12 <- result[,7]; beta22 <- result[,8]; beta32 <- result[,9]
R> beta13 <- result[,10]; beta23 <- result[,11]; beta33 <- result[,12]
R> lambda1 <- result[,13]; lambda2 <- result[,14]; lambda3 <- result[,15]

Table 4 shows posterior summaries for the illness-death model parameters using heart2 data.

As discussed in Section 5, the transition probabilities (16)–(21) are the most appropriate way to analyse the evolution of each state over time. The implementation of the posterior distribution for these transition probabilities requires auxiliary functions, similar to the calculation of the cumulative incidence function in Section 4. To avoid unnecessary repetitions on how to calculate these quantities of interest, we will omit their implementation here. However, the code to reproduce Figure 5 is available in Appendix A.

We consider a kidney infection dataset, referred to as kidney. It is available from the frailtyHL package (Ha et al. 2018):

R> library("frailtyHL")
R> data("kidney")
R> str(kidney)
’data.frame’:Ψ76 obs. of  10 variables:
 $ id     : num  1 1 2 2 3 3 4 4 5 5 ...
 $ time   : num  8 16 23 13 22 28 447 318 30 12 ...
 $ status : num  1 1 1 0 1 1 1 1 1 1 ...
 $ age    : num  28 28 48 48 32 32 31 32 10 10 ...
 $ sex    : num  1 1 2 2 1 1 2 2 1 1 ...
 $ disease: Factor w/ 4 levels "Other","GN","AN",..: 1 1 2 2 1 1 1 1 1 1 ...
 $ frail  : num  2.3 2.3 1.9 1.9 1.2 1.2 0.5 0.5 1.5 1.5 ...
 $ GN     : num  0 0 1 1 0 0 0 0 0 0 ...
 $ AN     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ PKD    : num  0 0 0 0 0 0 0 0 0 0 ...

This dataset consists of times to the first and second recurrences of infection in 38 kidney patients using a portable dialysis machine. Infections can occur at the location of insertion of the catheter. The catheter is later removed if infection occurs and can be removed for other reasons, in which case the observation is censored (McGilchrist and Aisbett 1991). The following variables were repeated twice for each patient:

• •

  id: patient number.

• •

  time: time (in days) from insertion of the catheter to infection in kidney patients using portable dialysis machine.

• •

  status: censoring indicator (1: if patient is uncensored; 0: otherwise).

• •

  age: age (in years) of each patient.

• •

  sex: sex of each patient (1: if patient is male; 2: if patient is female).

• •

  disease: disease type (GN; AN; PKD; Other).

• •

  frail: frailty estimate from original paper.

• •

  GN: indicator for disease type GN.

• •

  AN: indicator for disease type AN.

• •

  PKD: indicator for disease type PKD.

Note that the status variable equal to zero in any of the recurrences means that until that moment of observation there was no infection. We would also like to mention that we will use a modelling based on the clock-reset approach (Kleinbaum and Klein 2012), in which time starts at zero again after each recurrence.

Time, in days, from insertion of a catheter in patient $i$th to infection is modelled trough a proportional hazard specification with multiplicative frailties:

Our variable of interest is time, which represents the time from insertion of the catheter to infection. The status covariate plays an important role in the codification of the survival and censoring times:

R> # Number of patients and catheters
R> n <- length(unique(kidney$id))
R> J <- 2
R> # Survival and censoring times
R> time <- kidney$time
R> cens <- time
R> time[kidney$status == 0] <- NA # Censored
R> is.censored <- as.numeric(is.na(time))
R> # Matrix format
R> time <- matrix(time, n, J, byrow = TRUE)
R> cens <- matrix(cens, n, J, byrow = TRUE)
R> is.censored <- matrix(is.censored, n, J, byrow = TRUE)

Without loss of generality, we have created a design matrix X with the sex covariate:

R> sex <- kidney$sex[seq(1, 2 * n, 2)] - 1 # Reference = male
R> X <- model.matrix(~ sex)

Listing 6 shows a generic implementation of a frailty model in BUGS syntax using kidney data.