Tutorial on anoint package

FDA PACES Workshop, June 25, 2013

Stephanie Kovalchik
Research Fellow, National Cancer Institute

anoint, for ANalysis Of INTeractions, is a package written in the R language, which provides a suite of tools for investigating heterogeneity of treatment effect in a clinical trial.

Description of anoint Package

  • Open-source software

  • Written in the R language

  • Provides methods for assessing heterogeneity of treatment effect, including:

    • Proportional interactions modeling
    • Unrestricted multiple interaction regression
    • Conventional subgroup analyses
    • Forest plots

First...A Very Brief Introduction to R

  • R is a statistical programming environment

  • It is maintained by the R Development Core Team

  • It can be run on Windows, Linux, and Mac OS platiforms

  • It can be downloaded from http://cran.r-project.org/

You Can Use R as a Calculator

(1:3)^2

## [1] 1 4 9

(1:3) * (4:6)

## [1]  4 10 18

You Can Use R to Store & Manipulate Data

object <- data.frame(x = 1:6, y = rep(1, 6))

object

##   x y
## 1 1 1
## 2 2 1
## 3 3 1
## 4 4 1
## 5 5 1
## 6 6 1


object[1, ]

##   x y
## 1 1 1

You Can Use R to Import Data

thedata <- read.csv("filename")

You Can Use R to Perform Statistical Analyses of All Kinds

quantile(object$x)

##   0%  25%  50%  75% 100% 
## 1.00 2.25 3.50 4.75 6.00

You Can Use R to Perform Statistical Analyses of All Kinds

qqnorm(rnorm(100))

plot of chunk unnamed-chunk-5

You Can Extend R by Writing Functions

is.missing <- function(x) NA %in% x

is.missing(object$x)

## [1] FALSE

sapply(object, is.missing)

##     x     y 
## FALSE FALSE

You Can Also Extend R by Installing Contributed Packages

  • The rapid development of R is owing to the contributions of a growing community of programmers (mostly academics), who are willing to share their work with others.

  • After passing some quality checks, these packages are posted to the Comprehensive R Archive Network (CRAN), which can be accessed directly from the R environment.

install.packages("package-name")

Getting Started

  1. Install R (http://cran.r-project.org/)

  2. Install anoint package from CRAN (Version 1.3)

  • Run R
  • If you have an active web connection, use the following code:
install.packages("anoint")
  • If needed, install any dependent packages

Getting Started with anoint

All sessions with the anoint begin by loading the package.

library(anoint)
  • anoint's dependent packages will need to be installed

Getting Help

An index of all of the functions of anoint can be obtained by using the help function.

Click the link for any function to bring up documentation describing its use.

help(package = "anoint")

Getting Help

You can also use `?' as a shortcut to a function's documentation.

`?`(anoint)

Overview of anoint package

Function Description
anoint Creates an analysis of interactions object
obo Extract one-by-one subgroup analyses from anoint object
uim Extract unrestricted multiple interaction regression from anoint object
pim.subsets Perform all subsets proportional interactions modeling
pim.fit Fit a specific proportional interactions model
forest.subsets Forest plot of all subsets procedure

Examples with bootstrap SOLVD-T

We will demonstrate the major tools of the anoint package with a bootstrap SOLVD-T dataset.

If this is located in the current directory, it can be loaded into the R session with the function load.

setwd("~/master/project/anoint/tutorial")  # On My System

load("solvd.RData")

Description of SOLVD-T Dataset

str(solvd)  # Show Structure

## 'data.frame':    2569 obs. of  10 variables:
##  $ enal      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ deathdays : num  723 835 1547 1051 1449 ...
##  $ death     : num  0 0 0 0 0 1 1 0 0 0 ...
##  $ vaso      : int  0 0 1 0 1 0 0 1 0 1 ...
##  $ ischemia  : int  1 1 1 1 1 0 1 1 0 1 ...
##  $ sodium    : int  144 140 142 140 135 143 139 136 142 136 ...
##  $ lvef      : int  35 28 24 26 24 15 17 35 22 16 ...
##  $ nyha      : int  2 2 3 3 3 2 2 2 3 2 ...
##  $ sodium.cat: Ord.factor w/ 3 levels "<140"<"140-141"<..: 3 2 3 2 1 3 1 1 3 1 ...
##  $ lvef.cat  : Ord.factor w/ 3 levels "<23"<"23-29"<..: 3 2 2 2 2 1 1 3 1 1 ...

Description of SOLVD-T Variables

  • enal: Indicator of enalapril treatment group (1=enalapril, 0=placebo)
  • deathdays: Days to death or last date of follow-up if still alive
  • death: Indicator of death
  • vaso: Indicator of prior vasodilator use (1=use, 0=no use)
  • ischemia: Indicator of ischemic cause of congestive heart failure
  • sodium: Sodium level (mmol/liter)
  • lvef: Left ventricular ejection fraction (as percent)
  • nyha: NYHA prognostic class (1 - 4, 1=best and 4=worst)

Syntax of anoint Object

anoint(formula, data, family, ...)
Argument Description
formula Model formula
data Data frame containing formula variables
family coxph or any of glm families, i.e. gaussian, binomial, etc.

anoint Formula Specification

The formula used to create an anoint object requires

  • a response (appropriate for the model)

  • covariates

  • a treatment variable

Example: anoint Formula Specification

Suppose:

  • response = y

  • covariates = a, b, c

  • treatment variable = trt

formula = y ~ (a + b + c) * trt

Example: anoint

fit <- anoint(
    formula = Surv(deathdays, death) ~ (vaso + ischemia + lvef + sodium + nyha) * enal, 
    data = solvd, 
    family = "coxph")

Example: anoint

class(fit)

## [1] "anoint"
## attr(,"package")
## [1] "anoint"

Methods for anoint Object

Method Description
print Shows formula
summary Shows formula, covariates
obo Returns results of one-by-one subgroup analyses
uim Returns unrestricted interaction model

Example: print Method for anoint

fit

## Surv(deathdays, death) ~ (vaso + ischemia + lvef + sodium + nyha) * 
##     enal

The obo Method for anoint

Returns a list with the following named components:

  • fit: list of model objects; length equal to the number of covariates

  • LRT: vector of likelihood ratio test (LRT) statistics for interaction

  • pvalue: vector of unadjusted p-values for LRTs

Example: obo Method

subgroups <- obo(fit)

names(subgroups)

## [1] "fit"    "LRT"    "pvalue"

Example: obo Method

subgroups$LRT

##     vaso ischemia     lvef   sodium     nyha 
##   2.4068   1.4669   4.1930   0.2081   3.7393

subgroups$pvalue

##     vaso ischemia     lvef   sodium     nyha 
##  0.12081  0.22584  0.04059  0.64828  0.05315

Example: obo method

subgroups$fit[[1]]  ## Cox model for vasodilator

## Call:
## coxph(formula = f, data = anoint@data)
## 
## 
##               coef exp(coef) se(coef)       z     p
## vaso      -0.00504     0.995   0.0902 -0.0559 0.960
## enal      -0.17119     0.843   0.0968 -1.7678 0.077
## vaso:enal  0.20119     1.223   0.1297  1.5516 0.120
## 
## Likelihood ratio test=5.45  on 3 df, p=0.142  n= 2569, number of events= 966

Example: obo method

betas <- sapply(subgroups$fit, function(x) x$coef[1])  ## Covariate effects in placebo group
betas

##         Surv(deathdays, death) ~ vaso *  enal.vaso 
##                                          -0.005044 
## Surv(deathdays, death) ~ ischemia *  enal.ischemia 
##                                          -0.057539 
##         Surv(deathdays, death) ~ lvef *  enal.lvef 
##                                          -0.047876 
##     Surv(deathdays, death) ~ sodium *  enal.sodium 
##                                          -0.040647 
##         Surv(deathdays, death) ~ nyha *  enal.nyha 
##                                           0.549418

Example: obo Method

interactions <- sapply(subgroups$fit, function(x) x$coef[3])  # Interaction
interactions

##         Surv(deathdays, death) ~ vaso *  enal.vaso:enal 
##                                                0.201191 
## Surv(deathdays, death) ~ ischemia *  enal.ischemia:enal 
##                                                0.174479 
##         Surv(deathdays, death) ~ lvef *  enal.lvef:enal 
##                                                0.019561 
##     Surv(deathdays, death) ~ sodium *  enal.sodium:enal 
##                                                0.009408 
##         Surv(deathdays, death) ~ nyha *  enal.nyha:enal 
##                                               -0.189664

Example: obo Method

interactions/betas + 1  # implied proportional effects

##         Surv(deathdays, death) ~ vaso *  enal.vaso:enal 
##                                                -38.8847 
## Surv(deathdays, death) ~ ischemia *  enal.ischemia:enal 
##                                                 -2.0324 
##         Surv(deathdays, death) ~ lvef *  enal.lvef:enal 
##                                                  0.5914 
##     Surv(deathdays, death) ~ sodium *  enal.sodium:enal 
##                                                  0.7685 
##         Surv(deathdays, death) ~ nyha *  enal.nyha:enal 
##                                                  0.6548

The uim Method for anoint

Returns a list with the same named components as obo:

  • fit: fitted model object

  • LRT: value of global likelihood ratio test (LRT) for any interaction

  • pvalue: value of corresponding LRT

Example: uim Method

unrestricted <- uim(fit)  # Cox multiple interaction model

names(unrestricted)

## [1] "fit"    "LRT"    "pvalue"

Example: uim Method

unrestricted$LRT

## [1] 7.693

unrestricted$pvalue

## [1] 0.174

Example: uim Method

unrestricted$fit

## Call:
## coxph(formula = object@formula@formula, data = object@data)
## 
## 
##                    coef exp(coef) se(coef)         z       p
## vaso          -3.15e-02    0.9690  0.09271 -0.339303 7.3e-01
## ischemia       1.16e-05    1.0000  0.10174  0.000114 1.0e+00
## lvef          -3.82e-02    0.9626  0.00677 -5.636172 1.7e-08
## sodium        -4.07e-02    0.9601  0.01434 -2.837351 4.5e-03
## nyha           4.64e-01    1.5897  0.06782  6.835390 8.2e-12
## enal          -2.45e+00    0.0863  2.89633 -0.845774 4.0e-01
## vaso:enal      1.73e-01    1.1889  0.13286  1.302541 1.9e-01
## ischemia:enal  1.44e-01    1.1547  0.14823  0.970554 3.3e-01
## lvef:enal      9.85e-03    1.0099  0.00974  1.010600 3.1e-01
## sodium:enal    1.65e-02    1.0167  0.02062  0.801963 4.2e-01
## nyha:enal     -1.44e-01    0.8655  0.09909 -1.457243 1.5e-01
## 
## Likelihood ratio test=154  on 11 df, p=0  n= 2569, number of events= 966

Example: uim Method

betas <- unrestricted$fit$coef[1:5]  # covariate placebo effects
betas

##       vaso   ischemia       lvef     sodium       nyha 
## -3.146e-02  1.157e-05 -3.815e-02 -4.069e-02  4.636e-01

Example: uim Method

interactions <- unrestricted$fit$coef[7:11]  # interaction effects
interactions

##     vaso:enal ischemia:enal     lvef:enal   sodium:enal     nyha:enal 
##      0.173049      0.143865      0.009846      0.016537     -0.144392

Example: uim Method

interactions/betas + 1  # implied proportional effect

##     vaso:enal ischemia:enal     lvef:enal   sodium:enal     nyha:enal 
##       -4.5014    12438.0717        0.7419        0.5936        0.6885

Proportional Interactions All Subsets

  • We see from the subgroup and unrestricted interaction effects that there is some suggestion of proportionality among the SOLVD-T candidate effect modifiers.

  • To investigate this formally, we can use the function pim.subsets

Syntax of pim.subsets Function

The syntax of the pim.subsets function is very similar to the function anoint.

pim.subsets(formula, data, family, trt, ...)

Syntax of pim.subsets Function

However, rather than specifying the treatment variable as part of formula, the name of the treatment variable is indicated in the argument trt.

pim.subsets(formula, data, family, trt, ...)

Example: pim.subsets Function

fit.subsets <- pim.subsets(
    formula = Surv(deathdays, death) ~ vaso + ischemia + lvef + sodium + nyha, 
    trt = "enal", data = solvd, family = "coxph")

##    subset interaction num      LRT    lower   upper    pvalue reject
## 9   01101      0.7434   3  8.90024   0.5952  0.9016 0.0028513   TRUE
## 24  11101      0.7343   4  9.42195   0.5937  0.8863 0.0021440   TRUE
## 2   00101      0.7331   2  9.65005   0.6021  0.8891 0.0018935   TRUE
## 16  10101      0.7245   3 10.15835   0.5874  0.8781 0.0014365   TRUE
## 11  01111      0.7261   4 10.78335   0.5929  0.8809 0.0010242   TRUE
## 26  11111      0.7163   5 11.50786   0.5897  0.8653 0.0006930   TRUE
## 4   00111      0.7150   3 11.65697   0.5799  0.8680 0.0006396   TRUE
## 18  10111      0.7061   4 12.29694   0.5736  0.8647 0.0004537   TRUE
## 13  10010      1.0887   2  0.01920   0.2238  5.9074 0.8897986  FALSE
## 21  11010      1.1774   3  0.06677   0.2932  5.1624 0.7961048  FALSE
## 6   01010      0.8401   2  0.08374   0.2047  3.3398 0.7722924  FALSE
## 23  11100      0.7154   3  1.79286   0.3812  1.1946 0.1805778  FALSE
## 25  11110      0.7049   4  2.35769   0.4067  1.1278 0.1246666  FALSE
## 15  10100      0.6615   2  2.67703   0.3518  1.0393 0.1018056  FALSE
## 8   01100      0.6609   2  2.68945   0.3408  1.0461 0.1010148  FALSE
## 19  11000     -8.2577   2  3.13799 -25.3649 -0.3323 0.0764882  FALSE
## 10  01110      0.6636   3  3.23664   0.3783  1.0168 0.0720084  FALSE
## 17  10110      0.6545   3  3.35297   0.3678  0.9975 0.0670837  FALSE
## 20  11001      0.6633   3  3.37300   0.3876  1.0071 0.0662730  FALSE
## 12  10001      0.6630   2  3.41665   0.3870  1.0331 0.0645417  FALSE
## 5   01001      0.6599   2  3.54650   0.3933  1.0205 0.0596713  FALSE
## 7   01011      0.6540   3  4.62729   0.4251  0.9492 0.0314672  FALSE
## 22  11011      0.6476   4  4.73537   0.4124  0.9786 0.0295486  FALSE
## 3   00110      0.6019   2  4.76715   0.3432  0.9735 0.0290078  FALSE
## 1   00011      0.6524   2  4.77758   0.4271  0.9901 0.0288327  FALSE
## 14  10011      0.6451   3  4.87995   0.4230  0.9415 0.0271704  FALSE

pim.subsets Returned Object

A list with the named components:

  • subset: indicators of included covariates
  • interaction: proportionality parameter estimates
  • LRT: LRTs for proportionality
  • num: number of included covariates
  • lower: lower endpoints of 95% CI for proportionality parameter
  • upper: upper endpoints of 95% CI for proportionality parameter
  • pvalue: p-values of corresponding LRT
  • covariates: covariate names
  • reject: vector of logical values indicating rejection of LRT at multiplicity corrected criterion (set with fwer argument)

Example: pim.subsets Function

fit.subsets$covariates

## [1] "vaso"     "ischemia" "lvef"     "sodium"   "nyha"

Example: pim.subsets Function

cbind(fit.subsets$subset, fit.subsets$LRT, fit.subsets$reject)

##       [,1]    [,2]                 [,3]   
##  [1,] "01101" "8.90023811334097"   "TRUE" 
##  [2,] "11101" "9.42194817028956"   "TRUE" 
##  [3,] "00101" "9.65005022459761"   "TRUE" 
##  [4,] "10101" "10.1583522349928"   "TRUE" 
##  [5,] "01111" "10.7833529645664"   "TRUE" 
##  [6,] "11111" "11.5078553943072"   "TRUE" 
##  [7,] "00111" "11.6569696787207"   "TRUE" 
##  [8,] "10111" "12.2969384312842"   "TRUE" 
##  [9,] "10010" "0.0191986181248966" "FALSE"
## [10,] "11010" "0.0667661675164961" "FALSE"
## [11,] "01010" "0.0837394419787868" "FALSE"
## [12,] "11100" "1.79286210618331"   "FALSE"
## [13,] "11110" "2.35769090031863"   "FALSE"
## [14,] "10100" "2.67702689475455"   "FALSE"
## [15,] "01100" "2.68944818922313"   "FALSE"
## [16,] "11000" "3.13798679990886"   "FALSE"
## [17,] "01110" "3.23663577543446"   "FALSE"
## [18,] "10110" "3.35297358414815"   "FALSE"
## [19,] "11001" "3.37299889754784"   "FALSE"
## [20,] "10001" "3.41665296927792"   "FALSE"
## [21,] "01001" "3.54649768064868"   "FALSE"
## [22,] "01011" "4.62728953837953"   "FALSE"
## [23,] "11011" "4.73536618393853"   "FALSE"
## [24,] "00110" "4.76714889609609"   "FALSE"
## [25,] "00011" "4.77757528714653"   "FALSE"
## [26,] "10011" "4.87995034896806"   "FALSE"

forest.subsets Plotting Method

The results of the all subsets procedure can be graphically displayed with a forest plot using the function forest.subsets.

The syntax is:

forest.subsets(object, index, labels, ...)

forest.subsets Plotting Method

forest.subsets(object, index, labels, ...)
Argument Description
object Result of pim.subsets or anoint.subgroups
index Vector specifying the PIM subsets to include
... Additional arguments for plot rendering

Example: forest.subsets Plotting Method

forest.subsets(fit.subsets, xlimits = c(0.5, 1.25))  # With limits for x-axis

plot of chunk unnamed-chunk-42

anoint.subgroups Function

Bonferroni-corrected subgroup analyses can be obtained with anoint.subgroups, which takes the same arguments as pim.subsets.

anoint.subgroups(formula, data, family, trt, ...)

Example: anoint.subgroups

fit.subgroups <- anoint.subgroups(
    formula = Surv(deathdays, death) ~ vaso + ischemia + lvef + sodium + nyha, 
    trt = "enal", data = solvd, family = "coxph")

##          subset interaction    LRT      lower    upper  pvalue covariates
## lvef      00100    0.019561 4.1930  0.0008367 0.038285 0.04059       lvef
## nyha      00001   -0.189664 3.7393 -0.3817445 0.002416 0.05315       nyha
## vaso      10000    0.201191 2.4068 -0.0529541 0.455337 0.12081       vaso
## ischemia  01000    0.174479 1.4669 -0.1080974 0.457055 0.22584   ischemia
## sodium    00010    0.009408 0.2081 -0.0310088 0.049825 0.64828     sodium
##          reject
## lvef      FALSE
## nyha      FALSE
## vaso      FALSE
## ischemia  FALSE
## sodium    FALSE

Example: anoint.subgroups

cbind(fit.subgroups$cov, fit.subgroups$LRT, fit.subgroups$reject)

##      [,1]   [,2] [,3]
## [1,]    2 4.1930    0
## [2,]    3 3.7393    0
## [3,]    5 2.4068    0
## [4,]    1 1.4669    0
## [5,]    4 0.2081    0

Example: Forest Plot of Subgroup Treatment-Covariate Interactions

forest.subsets(fit.subgroups, subgroup = TRUE)

plot of chunk unnamed-chunk-46

Fitting Specific PIMs with pim.fit

If we wanted to consider a single proportional interactions model, we can do it with the pim.fit function.

object <- pim.fit(
	formula = Surv(deathdays, death) ~ vaso + lvef + sodium + nyha, 
	trt = "enal", data = solvd, family = "coxph")

pim.fit Returned Object

A list with the named components:

  • interaction: proportional effect
  • LRT: LRT statistic
  • lower: lower endpoint of 95% CI for proportionality parameter
  • upper: upper endpoint of 95% CI for proportionality parameter
  • pvalue: p-value of LRT
  • model0: covariate model object in placebo group
  • model1: covariate model object in treatment group

Example: pim.fit Returned Object

object$interaction

## [1] 0.7061

object$LRT

## [1] 12.3

object$pvalue

## [1] 0.0004537

Example: pim.fit Returned Object

object$model1

## Call:
## coxph(formula = formula, data = treatment)
## 
## 
##           coef exp(coef) se(coef)     z       p
## vaso    0.1607     1.174  0.09393  1.71 8.7e-02
## lvef   -0.0268     0.974  0.00691 -3.88 1.0e-04
## sodium -0.0240     0.976  0.01481 -1.62 1.0e-01
## nyha    0.3263     1.386  0.07260  4.49 7.0e-06
## 
## Likelihood ratio test=46  on 4 df, p=2.51e-09  n= 1285, number of events= 467

References

  • R Core Team (2013). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/

  • Varadhan R, Kovalchik S (2013). anoint: Analysis of interactions. R package Version 1.3. http://CRAN.R-project.org/package=anoint

  • Kovalchik S, Weiss CO, Varadhan R (2013). Assessing heterogeneity of treatment effect in a clinical trial with the proportional interactions model. Statistics in Medicine, In press.