The interactionRCS
package is designed to facilitate
interpretation and presentation of results from a regression model
(linear, logistic, Cox) where an interaction between the main predictor
of interest X (binary or
continuous) and another continuous covariate Z has been specified. Specifically, the
package will provide point estimates of the main effect of X over levels of Z, allowing for settings where Z is flexibly modeled with restricted
cubic splines, and provide a graphical display of this interaction. Two
methods for deriving and plotting confidence intervals are also
implemented, including the delta method and bootstrap.
Functions within the interactionRCS
package require that
a regression model has already been estimated and model results be
provided as an object. Models with linear and log-linear interactions
can be run with standard function (i.e. glm
for linear and
logistic, coxph
for survival). To model interactions with
restricted cubic splines, however, interactionRCS
requires
models to be run with more flexible functions from the rms
package (specifically Glm
, lrm
, and
cph
).
The main function of interactionRCS
is
intEST
, which provides point estimates and confidence
intervals for the effect of X over
levels of Z when an interaction is
included in the model, either as a (log-)linear term or flexibly modeled
with rcs
. The following options must be specified:
model
: the model previously run
(Glm
,lrm
, cph
, or
glm
, coxph
from models without splines)var1
: the name of the main predictor of interest (X)var2
: the name of the continuous predictor interacting
with var1
(Z)var2values
: the values of var2 for which the HR of var1 should be calculatedAdditional options include:
data
: the same dataset used for fitting the model (only
used for bootstrap CIs). If data=NULL, we will search for model$xci
(default TRUE) : whether a confidence interval for
each effect estimate should also be providedci.method
: either "delta"
or
"bootstrap"
. Default "delta"
ci.boot.method
(default= “percentile” - only if
method="bootstrap"
) : see boot.ci
type
parameterR
(default=100 - only if
method="bootstrap"
) : number of bootstrap iterationsparallel
(default “multicore” - only if
method="bootstrap"
): see boot.ci
referenceAn additional function plotINT
is also implemented to
provide a graphical display of the results and only require the
rcs-
, loglin-
, or lin-
results as
object.
The first example is based on a study on drug relapse among 575 patients enrolled in a clinical trial of residential treatment for drug abuse. The main exposure of interest is the binary indicator of assigned treatment (0/1) and a treatment*age interaction is specified.
The following code provides an estimate of the treatment HR at
different ages, when age is modeled with restricted cubic splines. The
model is further adjusted for tumor site, race, and previous use of IV
drug. It is recommended to check the distribution of Z (here, age) to define a realistic range
of var2values
. Finally, the code is replicated by using the
delta method or bootstrap for obtaining confidence intervals. Note that
when ci.method = "bootstrap"
is specified, additional
options can be specified.
myformula <- Surv(time, censor) ~ treat*rcs(age, 3) + site + nonwhite + ivdrug
model_rcs <- cph(myformula , data = umaru , x = TRUE , y=TRUE)
HR_rcs_delta <- intEST( var2values = c(20:50)
, model = model_rcs , data = umaru , var1 ="treat", var2="age" ,ci.method = "delta")
plotINT(HR_rcs_delta , xlab = "Age")
# note that the user could directly specify the function of interest (here rcsHR) and obtain the same results
HR_rcs_delta <- rcsHR( var2values = c(20:50)
, model = model_rcs , data = umaru , var1 ="treat", var2="age" ,ci.method = "delta")
HR_rcs_boot <- intEST( var2values = c(20:50)
, model = model_rcs , data = umaru , var1 ="treat", var2="age" , ci.method = "bootstrap"
, ci.boot.method = "norm" , R = 500 , parallel = "multicore")
plotINT(HR_rcs_boot , xlab = "Age")
The following code replicates the same analysis in the setting where
the interaction model was specified without modeling the continuous
covariate Z with restricted cubic
splines (log-linear interaction model). The loglinHR
function used in this context requires the same options of
rcsHR
myformula <- Surv(time, censor) ~ treat*age + site + nonwhite + ivdrug
model_loglin <- cph(myformula , data = umaru , x = TRUE , y=TRUE)
HR_loglin_delta <- intEST( var2values = c(20:50)
, model = model_loglin , data = umaru , var1 ="treat", var2="age", ci.method = "delta")
plotINT(HR_loglin_delta , xlab = "Age")
HR_loglin_boot <- intEST( var2values = c(20:50)
, model = model_loglin , data = umaru , var1 ="treat", var2="age", ci.method = "bootstrap"
, ci.boot.method = "norm" , R = 500 , parallel = "multicore")
plotINT(HR_loglin_boot , xlab = "Age")
In this example the log-linear and spline interaction models provide
similar results, with the treatment HR decreasing over levels of age
(note that, in regular practice, a thorough model comparison to decide
whether the log-linear or spline interaction model better fits the data
should be conducted before using interactionRCS
).
This second example is based on a larger dataset of 17549 individuals from a population study of non-alcoholic fatty liver disease (NAFLD). Here the main exposure X of interest is age, and an interaction with BMI is specified. The code will be the same in the presence of a continuous covariate, but the interpretation will be slightly different as the functions will estimate the HR for a unit increase in X (age) over levels of BMI. The following code will provide estimates of the HR for a unit increase of age over levels of BMI, modeled with restricted cubic spline.
myformula <- Surv(futime, status) ~ age*rcs(bmi, 3) + male
model2_rcs <- cph(myformula , data = nafld1 , x = TRUE , y=TRUE)
HR2_rcs_delta <- intEST( var2values = c(15:50)
, model = model2_rcs , data = nafld1 , var1 ="age", var2="bmi" , ci=TRUE , conf = 0.95 , ci.method = "delta")
plotINT(HR2_rcs_delta , xlab = "BMI")
If the user is interested in a different HR (e.g. the HR for a 10-unit increase, or the HR for a 1-sd unit increase), the predictor should be standardized accordingly before fitting the model. For example, here we estimate the HR for a 5-years increase in age
nafld1$age5<-nafld1$age/5
myformula <- Surv(futime, status) ~ age5*rcs(bmi, 3) + male
model3_rcs <- cph(myformula , data = nafld1 , x = TRUE , y=TRUE)
HR3_rcs_delta <- intEST( var2values = c(15:50)
, model = model3_rcs , data = nafld1 , var1 ="age5", var2="bmi" , ci=TRUE , conf = 0.95 , ci.method = "delta")
plotINT(HR3_rcs_delta , xlab = "BMI")
The effect of age changes over levels of BMI in a quadratic form, with a lower effect for both those with increasingly higher and lower BMI. Including BMI without any spline transformation (i.e. the log-linear interaction model) would fail to capture this functional form, as shown from the following code. A comparison of the models through the use of AIC would also confirm that the spline interaction model better fits the data
myformula <- Surv(futime, status) ~ age*bmi + male
model2_loglin <- cph(myformula , data = nafld1 , x = TRUE , y=TRUE)
HR2_loglin_delta <- intEST( var2values = c(15:50)
, model = model2_loglin , data = nafld1 , var1 ="age", var2="bmi" , ci=TRUE , conf = 0.95 , ci.method = "delta")
plotINT(HR2_loglin_delta , xlab = "BMI")
AIC(model2_loglin)
## [1] 15979.93
AIC(model2_rcs)
## [1] 15924.4
interactionRCS
requires results from a regression model
where an interaction between a main predictor (binary or continuous)
X and a continuous predictor Z has been specified. This interaction
can be included as a simple product term between the 2 predictors, or by
flexibly modeling Z with restricted
cubic splines. For both interaction settings, the main exposure of
interest X has to either be binary
or continuous.
A basic Cox model with 2 predictors and their interaction takes the form:
h(t|x,z)=h0⋅exp(β1x+β2z+β3x⋅z)
After estimation of the model, we want to predict the HR for X over levels of Z. This is given by
HR10=h(t|x=1,z)h(t|x=0,z)=h0⋅exp(β1x+β2z+β3x⋅z)h0⋅exp(β1x+β2z+β3x⋅z)=exp(β1+β2z+β3z)exp(β2z) ,
which will be plotted against Z
To estimate 95% confidence intervals SE(HR10) is required. This is obtained by focusing on log(HR10) and calculating lower and upper bounds for this quantity, which are then exponentiated.
SE(log(HR10))=SE(log(exp(β1+β2z+β3z)exp(β2z)))=SE(log(exp(β1+β2z+β3z)−log(exp(β2z)))=SE(β1+β3z)
This is calculated by using the delta method. Upper and lower bounds are then derived with exp(log(HR10)±1.96⋅SE(log(HR10)))
A logistic regression model with 2 predictors and their interaction takes the form:
logit(p|x,z)=β0+β1x+β2z+β3x⋅z
After estimation of the model, we want to predict the OR for X over levels of Z. This is given by
OR10=odds(p|x=1,z)odds(p|x=0,z)=exp(β0+β1x+β2z+β3x⋅z)exp(β0+β1x+β2z+β3x⋅z)=exp(β0+β1+β2z+β3z)exp(β0+β2z) ,
which will be plotted against Z
To estimate 95% confidence intervals SE(OR10) is required. This is obtained by focusing on log(OR10) and calculating lower and upper bounds for this quantity, which are then exponentiated.
SE(log(HR10))=SE(log(exp(β0+β1+β2z+β3z)exp(β0+β2z)))=SE(log(exp(β0+β1+β2z+β3z)−log(exp(β0+β2z)))=SE(β1+β3z)
This is calculated by using the delta method. Upper and lower bounds are then derived with exp(log(OR10)±1.96⋅SE(log(OR10)))
A linear regression model with 2 predictors and their interaction takes the form:
E[Y|x,z]=β0+β1x+β2z+β3x⋅z
After estimation of the model, we want to predict the effect for X over levels of Z. This is given by
E[Y|x=1,z]−E[Y=0|x,z]=(β0+β1+β2z+β3⋅z)−(β0+β2z)=β1+β3⋅z ,
which will be plotted against Z
To estimate 95% confidence intervals we need the basic linear combination of standard errors calculated by SE(β1+β3z)=√SE(β1)2+SE(β3z)2+2cov(β1,β3)z,
or by using the Delta method
interactionRCS
allows the continuous covariate Z to be flexibly modeled with restricted
cubic splines, with any given of number (p) of knots (t1,t2,t3,…,tp).
The interaction model, in this setting, takes the form:
h(t|x,z,c)=h0⋅exp(β1x+sp(z)+sp2(z⋅x))
If 3 knots are specified:
sp(z)=α1z+α2⋅[(z−t1)3+−(z−t2)3+⋅t3−t1t3−t2+(z−t3)3+⋅t2−t1t3−t2(t3−t1)2]
and
sp2(z⋅x)=γ1x+γ2x⋅[(z−t1)3+−(z−t2)3+⋅t3−t1t3−t2+(z−t3)3+⋅t2−t1t3−t2(t3−t1)2]
After estimation of the model, we want to predict the HR for X over levels of Z. This is given by
HR10=h(t|X=x1,Z)h(t|X=x2,Z)=h0⋅exp(β1x1+sp(z)+sp2(z)⋅x1)h0⋅exp(β1x2+sp(z)+sp2(z)⋅x2)=h(t|X=x1,Z)h(t|X=x2,Z)=exp(β1+sp(z)+sp2(z))exp(sp(z)) ,
which will be plotted against Z
Similarly to the previous situation, the SE can be derived by using the delta method to calculate SE(β1+sp2(z))
The interaction model, in this setting, takes the form:
logit(p|x,z)=β0+β1x+sp(z)+sp2(z⋅x)
After estimation of the model, we want to predict the HR for X over levels of Z. This is given by
OR10=exp(β0+β1+sp(z)+sp2(z))exp(β0+sp(z))=exp(β1+sp2(z)) ,
which will be plotted against Z
Similarly to the previous situation, the SE can be derived by using the delta method to calculate SE(β1+sp2(z))
The interaction model, in this setting, takes the form:
E[Y|x,z]=β0+β1x+sp(z)+sp2(z⋅x)
After estimation of the model, we want to predict
E[Y|x=1,z]−E[Y=0|x,z]=(β0+β1+sp(z)+sp2(z))−(β0+sp(z))=β1+sp2(z) ,
which will be plotted against Z
The SE can be derived by using the delta method to calculate SE(β1+sp2(z))
Based on Harrell’s book (chapter 2-23), we derive equations for more than 3 knots
The general splines function for k knots is:
f(X)=β0+β1⋅X1+β2⋅X2+⋯+βk−1⋅Xk−1
where X1=X and for j=1,…,k−2,
Xj+1=(z−tj)3+−(z−tk−1)3+⋅tk−tjtk−tk−1+(z−tk)3+⋅tk−1−tjtk−tk−1(tk−t1)2
The Cox model with k=4 will therefore be:
h(t|x,z,c)=h0⋅exp(β1x+sp(z)+sp2(z⋅x))
where
sp(z)=α1z+α2⋅[(z−t1)3+−(z−t3)3+⋅t4−t1t4−t3+(z−t3)3+⋅t3−t1t4−t3(t4−t1)2]+α3⋅[(z−t2)3+−(z−t3)3+⋅t4−t2t4−t3+(z−t4)3+⋅t3−t2t4−t3(t4−t1)2]
and
sp2(z⋅x)=γ1x+γ2x⋅[(z−t1)3+−(z−t3)3+⋅t4−t1t4−t3+(z−t4)3+⋅t3−t2t4−t3(t3−t1)2]+γ3x⋅[(z−t2)3+−(z−t3)3+⋅t4−t2t4−t3+(z−t4)3+⋅t3−t2t4−t3(t4−t1)2]
The generic Cox model will be:
h(t|x,z,c)=h0⋅exp(β1x+sp(z)+sp2(z⋅x))
where
sp(z)=α1z+α2⋅[(z−tj)3+−(z−tk−1)3+⋅tk−tjtk−tk−1+(z−tk)3+⋅tk−1−tjtk−tk−1(tk−t1)2]+α3⋅[(z−tj)3+−(z−tk−1)3+⋅tk−tjtk−tk−1+(z−tk)3+⋅tk−1−tjtk−tk−1(tk−t1)2]
and
sp2(z⋅x)=γ1x+γ2x⋅[(z−tj)3+−(z−tk−1)3+⋅tk−tjtk−tk−1+(z−tk)3+⋅tk−1−tjtk−tk−1(tk−t1)2]+γ3x⋅[(z−tj)3+−(z−tk−1)3+⋅tk−tjtk−tk−1+(z−tk)3+⋅tk−1−tjtk−tk−1(tk−t1)2]