Journal of Clinical Epidemiology 65 (2012) 855e862

Visualizing interaction effects: a proposal for presentation

and interpretation
Claudia Lamina*, Gisela Sturm, Barbara Kollerits, Florian Kronenberg
Division of Genetic Epidemiology, Department of Medical Genetics, Molecular and Clinical Pharmacology, Innsbruck Medical University,
Schopfstr. 41, A-6020 Innsbruck, Austria
Accepted 15 February 2012; Published online 30 May 2012

Abstract
Objective: Interaction terms are often included in regression models to test whether the impact of one variable on the outcome is modified by another variable. However, the interpretation of these models is often not clear. We propose several graphical presentations and
corresponding statistical tests alleviating the interpretation of interaction effects.
Study Design and Setting: We implemented functions in the statistical program R that can be used on interaction terms in linear,
logistic, and Cox Proportional Hazards models. Survival data were simulated to show the functionalities of our proposed graphical visualization methods.
Results: The mutual modifying effect of the interaction term is grasped by our presented figures and methods: the combined effect of
both continuous variables is shown by a two-dimensional surface mimicking a 3D-Plot. Furthermore, significance regions were calculated
for the two variables involved in the interaction term, answering the question for which values of one variable the effect of the other variable
significantly differs from zero and vice versa.
Conclusion: We propose several graphical visualization methods to ease the interpretation of interaction effects making arbitrary
categorizations unnecessary. With these approaches, researchers and clinicians are equipped with the necessary information to assess
the clinical relevance and implications of interaction effects. 2012 Elsevier Inc. All rights reserved.
Keywords: Interaction; Visualization; Categorization; Cox models; Linear models; Logistic models

1. Background
In epidemiology, interaction terms are often incorporated into multivariable regression models to test whether
the impact of one variable on the outcome is modified by
another variable. However, the interpretation of interaction
effects is often not clear. It depends on the scale of the variables included in the interaction term (continuous and/or
categorical) and the scale of the model (linear regression,
logistic regression, Cox Proportional Hazards model). The
interpretation is rather straightforward, if interaction effects
between the two categorical variables or between one continuous and one categorical variable are considered. One
such example would be an effect that can only be seen in
men, but not in women, in Europeans but not in Asians
etc. For such situations, there are also graphical solutions implemented in standard statistical programs [1]. If

* Corresponding author. Tel.: 43-512-9003-70565; fax: 43-5129003-73561.

E-mail address: [email protected] (C. Lamina).
0895-4356/$ - see front matter 2012 Elsevier Inc. All rights reserved.
doi: 10.1016/j.jclinepi.2012.02.013

researchers are interested in interaction between two continuous variables, the interpretation is less clear. It seems
that interaction terms between continuous variables are
avoided or continuous variables are categorized beforehand. If significant P-values of continuous interaction terms
are reported, most authors switch to categorizations of these
former continuous variables for interpretation purposes. In
general, categorizations of continuous variables are often
done rather arbitrarily by choosing percentiles of the respective interacting variables, for example, by dichotomizing using the median [2]. Categorization of continuous
variables leads to significant loss of information and power,
though [3]. Sometimes, clinically recognized and predefined cutpoints are chosen, such as a body mass index of
O30 kg/m2 for the definition of obesity. Even in this case,
information is lost because the definition of these cutpoints
has been designed for a totally different purpose. This approach can lead to considerable bias. The problem even
worsens, if both continuous variables that constitute to an
interaction term are categorized [4].
In other cases, regression coefficients and P-values of interaction terms are shown, but the reader is left alone with the

856

C. Lamina et al. / Journal of Clinical Epidemiology 65 (2012) 855e862

What is new?
1. Epidemiological publications on interaction effects in survival analyses mostly use categorizations of variables for the ease of interpretation.
This approach might lead to considerable bias,
though.
2. To ease the interpretation of interaction effects between two continuous variables, we propose
a graphical visualization approach which we implemented in functions in the statistical program
R. They can be used on Generalized Linear
Models and Cox Proportional Hazards models.
3. The mutual modifying effect of both variables
constituting the interaction effect can be shown
by a two-dimensional twisted surface and interaction plots showing the effect estimator of one variable for varying values of the other variable.

interpretation. The effects and P-values of the main effects

are also not interpretable without taking into account the
concurrent levels of the other constituting variable [5].
By including multiplicative interaction terms, conditional hypotheses are tested: An increase in one variable
is associated with an increase in the outcome variable,
when the second interacting variable is equal to a specific
value. Inevitably, all constitutive terms of an interaction
term should be included in the regression model also individually. With the inclusion of an interaction term, these
former average effects are conditional effects then and
change just by definition. Therefore, changes in effects
and/or P-values after including interaction terms cannot
be interpreted usefully. The conditional nature also applies
for the P-values of the conditional effects, which are given
in the output tables of statistical programs: they test,
whether the slope of one variable is significantly different
from zero, if the other constituting variable is zero. In most
cases, this test is not very useful, depending on the scale of
the variable.
The typical reader will not make the necessary linear
combinations to calculate the effect of the first variable
for specific values of the second variable. It would be possible, though, given the typical output table of a regression
model, including effect estimates, standard errors, and
P-values. However, it is not possible to conclude on the statistical significance of such a linear combination.
Therefore, it is hardly possible to interpret interaction
effects simply by looking at the output table from regression models. The results of interaction terms have to be presented in a different, more reader-friendly way.
Such examples can only rarely be found. We conducted
a literature search on survival analyses publications

including interaction terms between two continuous variables (see Appendix on the journals Web site at www.
jclinepi.com). It revealed that only one study group used
the complete range of both continuous variables to interpret
the interaction effects: The authors presented interaction
plots, showing the modifying effect of waist circumference
on the relationship between several parameters (cholesterol,
leptin, and adiponectin) and mortality [6,7].
The lack of such forms of presentations is certainly because of unfamiliarity with the methodological background
and lack of easy to use methods. Therefore, such methods
are required in the clinical epidemiological setting to avoid
common misinterpretations in interaction models.
To ease the interpretation of interaction effects between
two continuous variables, we propose several possibilities
for graphical presentation, which we implemented in functions in the statistical program R: These functions have
been implemented for linear regression, logistic regression,
and Cox Proportional Hazards models. To demonstrate the
functionalities, one data set including survival data has
been simulated and analyzed in Cox models including linear interaction effects between two continuous predictors.

2. Methods
2.1. Formulating and interpreting multiplicative
interaction terms in linear, logistic, and Cox models
For explanatory purposes, a linear regression model is
assumed first. Interaction terms between two continuous
variables x1 and x2 on the continuous outcome variable Y
can be modeled in the following way:
EYjx1 ; x2 5 Xb 5 b0 b1 x1 b2 x2 b3 x1  x2
1
with b0 being the intercept, b1 and b2 the conditional effects of x1 and x2 and b3 being the interaction effect on
the outcome Y. The linear combination on the right side
of the equation is the linear predictor Xb, to which further
covariates can be added. The regression coefficients for x1
and x2 are conditional effects because they depend on the
values of the other constituting variable:
A one unit change of x1 corresponds to a change of
b1 b3x2 on Y.
A one unit change of x2 corresponds to a change of
b2 b3x1 on Y.
This means that b1 itself is the effect of x1 on Y, if x2 is
equal to zero and vice versa. Thus, the interpretation of regression coefficients cannot be done without taking the
levels of the interacting variable into account. According
to Rothman [8], the statistical term interaction is used
to refer to departure from the underlying form of a statistical
model. For a linear model, a significant interaction term
thus implies deviation from additivity, meaning the additive
effect of the individual effects b1 x1 b2 x2 on the outcome.

C. Lamina et al. / Journal of Clinical Epidemiology 65 (2012) 855e862

The interpretation based on the direction of effect estimates

is only unambiguous, if both the conditional effect, for example of x1, and the interaction effect have the same effect
direction. In this case, the absolute value of the overall x1
effect increases for increasing values of x2. If all effect estimates are positive, this implies that the interaction leads to
more than an additive effect. If the interaction effect and
conditional effect have opposing algebraic signs, the conditional effect is attenuated by increasing values of the other
interacting variable and can even be switched into the other
direction, depending on the size of the corresponding estimates. Obviously, interpretation of interaction effects based
on the table of regression coefficients alone is hard to accomplish and requires at least some linear combinations
to be made.
The corresponding P-values are interpreted accordingly:
They also reflect the x1 effect on Y for x2 5 0 and vice versa. This may or may not be useful depending on the scale
and ranges of the variables. One may want to test, however,
whether the effect of x1 significantly differs from zero given
a specific value of x2. The equation for this simple slopes
test requires the regression coefficients to be known and
the corresponding variance-covariance matrix [5], which
can be given by the estimated regression model.
If you want to test the effect of x1 for a specific x2 value,
the slope to be tested would be given by:
slopex1 5 b1 b3 x2
and the corresponding standard error by
p
seslopex1 5 s11 2  x2 s13 x2 s33

with s11 being the variance of b1, s33 the variance of b3, and
s13 the covariance of b1 and b3.
Given the slope and its standard error, a t-test can then
be constructed:

857

effects in a linear model, which shows departure from

additivity.
The Cox Proportional Hazards model [10] models the
time to a specific event such as death accounting for censored observation. The hazard function h(t) is the instantaneous probability of an event at time t given that the
individual has survived until time t and is defined as
follows:
ht 5 h0 texpXb

with h0(t) being the underlying baseline hazard, which cannot be estimated. For each of the covariates in the linear predictor, the factor exp(bi) gives the so-called hazard ratio
(HR), which is a measure of relative risk. It characterizes
the shift in the hazard function as a result of a shift of one
unit in the corresponding variable xi. The interpretation on
interaction effects on the linear predictor Xb in a Cox model
can be done analogical to the linear regression model. For
Cox models, it might be more meaningful to look at the
change in HR, though. As for logistic regression models,
on the scale of the relative risk estimate, the inclusion of interaction terms tests the deviation from multiplicativity.
2.2. Description of the simulated survival data set
including interaction effects
We simulated survival data to explain the interpretation
of interaction effects and show the utilities of our graphical
presentations in a Cox model. One thousand right-censored
observations were generated with three continuous and normally distributed variables with mean 0 and variance 1. In
addition, one N(0,1)-distributed nuisance parameter has
been included, that is ignored in the following Cox model.
Therefore, the linear predictor was constructed as follows:

Linear predictor Xb 5  2  x1 0  x2 1  x1  x2
1  x3 N0; 1

with n being the number of cases and k the number of variables not including the intercept. The simple slopes test of
x2 effects for specific x1 values can be derived likewise.
The interpretations and formulas on effect estimates
given in this paragraph can be generalized to logistic regression or Cox models by substituting the expectation of
the outcome variable by a general function. The linear
predictor Xb5b0 b1 x1 b2 x2 b3 x1  x2 is then modeled by a monotone link function as the logit function
lnp=1  p for logistic regression models, with p being
the probability of the outcome. In this case, the intuitive
effect measure is the odds ratio (OR), which is estimated
by ORi 5 exp(bi) for variable xi. Thus, the regression function can be rewritten to p=1  p5expb0  expb1 x1 
expb2 x2  expb3 x1  x2 . On this scale, a significant
interaction effect would imply a deviation from multiplicativity of individual effects [9], in contrast to interaction

The effect estimates were chosen to represent two possible interaction scenarios: a strong main effect that attenuates for increasing values of the other interacting variable
is reflected by b1, whereas b2 5 0 reflects an effect, that
would not be observed in a model not including an interaction term because the x2 effect even changes direction for
varying values of x1.
The simulation was performed using function Simsurv of the library prodlim in R [11]. For each observation, a Weibull distributed survival time was generated with
baseline risk set to 0.2. The censoring time was simulated
to follow an exponential distribution with a baseline censoring risk set to 3. The baseline risks for censoring and survival were chosen to yield a high censoring rate. The
minimum of the two simulated times was taken as observation time. The event indicator was set to 1, if the survival
time was shorter than the censoring time, otherwise the indicator was set to 0.

slopex1
|tn  k  1
t5
seslopex1

858

C. Lamina et al. / Journal of Clinical Epidemiology 65 (2012) 855e862

2.3. Implementation of the graphical utilities and

calculations in R
Several functions have been implemented in R to alleviate the interpretation of multiplicative interaction effects.
They can be downloaded from http://www3.i-med.ac.at/
genepi/ under Analysis tools. First, a generalized linear
model (GLM) or Cox model has to be fitted using the
function glm or cph (from the Design library). Then, the
function PlotInteraction3D (or PlotInteraction3DCox) calculates predicted values of the linear predictor for either
the GLM or Cox model for given ranges of both variables
constituting the interaction term. Additional covariates and
also factor variables can be included in the models. By default, mean values are chosen for continuous covariates. For
the factor variables, one factor level has to be chosen, for
which predicted values should be calculated. The resulting
matrix of predicted values is spanned as a two-dimensional
surface into a three-dimensional grid. Alternatively, a colored contour plot can be given out. Consecutively, twodimensional slices can be cut out of this plot for specific
values of one of the constituting variables (function PlotInteraction2D). Each slice is the linear predictor of the
regression model for varying values of one variable keeping
the other variable fixed. Assuming a linear model, for example, this corresponds to the predicted outcome. Setting
x2 to the fixed value c, it can be derived from formula (1)
by y5b0 b1 x1 b2 c b3 cx1. The function PlotInteraction2D uses the R-function predict, which gives out
these fitted values and the corresponding standard errors
(se.fit). In addition, point-wise 95% confidence intervals
are calculated for each value of, for example, x1 by
y6z0:975  se:fit, which are plotted for the complete range
of x1 values, with z0:975 being the 97.5% quantile of the normal distribution (|1.96). It is recommended to look at those
plots for a set of values spanning the range of the variable,
for example, for the 5th, 25th, 50th, 75th, and 95th quantile.
By means of the functions getRegion or getRegionCox, significance regions can be calculated. For a given
range of values for x1, for example, the slope of the x2 variable with the corresponding P-value is printed. An interaction effect plot [12] is also given out, which shows how the
effect of one variable changes with the changing values of
the other variable. This corresponds to plotting the
effect estimate of one variable (slopex1 ) for a given value
of the other variable as in formula (1). The respective
point-wise 95% confidence intervals for the x1 effects
for given values of x2 are derived from (1) and (2) by
slopex1 6z0:975  seslopex1 . Both, effect estimates for x1

and corresponding confidence intervals are plotted for varying values of x2 over their complete range. For logistic and
Cox models, OR or HR can be given alternatively by exponentiating the effect estimates and the corresponding confidence intervals. Thus, two figures, one for each of the two
interacting variables, are sufficient to summarize the interaction effect over the complete range of both continuous
variables.

3. Results
Eighty-two percent of the simulated observations are
censored. The results of the simulated Cox model are given
in Table 1. The risk experiencing an event is significantly
decreasing for increasing values of x1 (HR 5 0.1521,
P ! 0.001), if x2 5 0. If x1 5 0, x2 does not have an effect
on the hazard. Other implications that can directly be drawn
from the numbers in the table, are that there is a significant
interaction effect, meaning, that the effects of x1 and x2 vary
for varying values of the other variable. Other slopes can be
calculated by linear combinations of the effect estimates:
The effect estimate of x2 for specific x1 values can be calculated by 0.0227 1.0378*(x1 value). For x1 5 0.6456
(525% quantile), this would be 0.6473, for x1 5 1.772
(595% quantile), it is 1.862. For the calculation of
P-values, however, specific covariances are needed from
the output of the model. Thus, they cannot be derived
directly from Table 1, but can be given out by applying
the test given in formula (4) implemented in the function
getRegionCox: P 5 4.4e  14 (for x1 5 0.6456) and
P ! 10e  20 (for x1 5 1.772).
A rough overview of the impact of the interaction can be
given by the twisted surface in Fig. 1: The highest risk can
be found at consecutive low values of both variables. The
variable x1 does have a negative impact on the risk for
low values of x2. For increasing values of x2, the x1 effect
is attenuated. For low values of x1, the variable x2 is positively associated with the outcome variable, which is time
to event in this example (increasing slope), whereas x2 is inversely associated for high values of x1 (decreasing slope).
The surface is twisted, indicating that the direction of the x2
effect even changes direction for varying values of x1. Alternatively, a colored contour plot (Fig. 2) can depict the
same information: The highest risk can be found in areas,
which are displayed in pink (both x1 and x2 low), whereas
the risk decreases for darkening blue color.
Fig. 3 shows slices through the three-dimensional plot
for specific values of x1 or x2. The left panel shows the x1

Table 1. True and estimated effect estimates of the conditional and interaction effects in the simulated Cox model
^
^
Variable
True effect, b
Estimated effect, b
se b
HR
95% CI of HR
x1
x2
x 1 * x2
x3

2
0
1
1

1.8834
0.0227
1.0378
0.8826

Abbreviations: HR, hazard ration; CI, confidence interval.

0.1092
0.0976
0.0883
0.0777

0.1521
1.023
2.823
2.417

[0.1228;
[0.8449;
[2.3744;
[2.0757;

0.1884]
1.2386]
3.3564]
2.8148]

P-value
!0.001
0.816
!0.001
!0.001

C. Lamina et al. / Journal of Clinical Epidemiology 65 (2012) 855e862

859

Fig. 1. Interaction surface plot on the linear predictor of the hazard

function. The z-axis shows the linear predictor of the hazard function
for varying values of x1 and x2, illustrating their interacting effect on
the simulated event variable.

Fig. 2. Contour plot displaying the linear predictor of the hazard function. The linear predictor of the hazard function for varying values of
x1 and x2 is displayed in different colors, ranging from blue (low risk)
to pink (high risk).

effect, given x2 is equal to the 5th, 25th, 75th, and 95th

quantile with its corresponding confidence intervals. The
right panel shows the x2 effects for the same given values
of x1. These plots give the impression of flipping through
a flip book, showing the twisting of effects in mutual dependence of the other constituting variable. Significance regions can be given out by the function getRegionCox.
The x1 effect is significantly different from zero for x2
values 1.55 (594.1% quantile). The x2 effect is significantly different from zero with a negative slope for x1
values !0.19 (542.16% quantile) and significantly different from zero with a positive slope for x1 values
O0.17 (555.9% quantile). The interaction plots in Fig. 4
show the HRs for x1 for varying x2 values with pointwise confidence intervals for each plotted value. It can be
seen that the HR is significantly lower than 1 for the vast
majority of x2 values. The significance region can also be
read from this plot. The x2 effect changes from an HR that
is smaller than 1 to higher than 1 almost exactly at the
median of x1 values. This plot can be printed for the linear
predictor, which is the default option, but also for the HR,
which has been chosen here in this example.

interpretation of an interaction effect: we observed a significant interaction between time-varying values of albumin
and phosphorus on overall mortality in our cohort [13]:
for high albumin values, increasing phosphorus concentrations were significantly associated with an increased mortality risk, whereas the beneficial effect of increasing
albumin values was only observed for low phosphorus
levels. We concluded that decisions about phosphoruslowering therapy should also take albumin values into account and epidemiological studies and guidelines should
consider this interplay. The proposed graphical presentation
and the calculation of the significance regions helped in the
discussion with clinicians and interpretation of results.

4. Discussion
We have implemented functions in R providing several
different graphical presentations and corresponding statistical tests to ease the interpretation of interaction effects between continuous variables. We have successfully applied
one of these approaches in a prospective study on incident dialysis patients, which was very useful for the

4.1. Categorizing continuous variables in interaction

terms: pros and cons
Most of the interaction models between continuous variables that are reported are broken down to subgroup analysis
and/or categorizations of at least one of the two continuous
variables. A qualitative assessment of epidemiological publications [2] found that 86% of papers evaluating continuous
risk factors on health outcomes categorize these risk factors.
Most of them use quantiles, others equally spaced categories
(e.g., 10-year-age categories) or external criteria, such as
clinically respected cutpoints. The reasons for the categorization approaches were not mentioned in most of the publications. There seems to be a gap between statistical
capabilities and the ability to interpret interaction effects
in a meaningful way. Therefore, simplicity of interpretation
might be the most crucial point. However, this simplicity is

860

C. Lamina et al. / Journal of Clinical Epidemiology 65 (2012) 855e862

Fig. 3. Slices through the interaction surface plot for specific values. The y-axis shows the linear predictor of the hazard function for varying values
of x1, holding x2 fixed at specific values (5%, 25%, 50%, 75%, and 95% quantiles; left panel) and the linear predictor of the hazard function for
varying values of x2, holding x1 fixed at specific values (5%, 25%, 50%, 75%, and 95% quantiles; right panel).

gained at the cost of power loss. Dichotomization of data is

even equal to discarding one third of the data [3]. Using cutpoint models in epidemiology is in best cases just unrealistic

and conservative. Considerable variability in the categories

is lost and artificial steps in risk are created. In regression
models where both originally continuous variables are

C. Lamina et al. / Journal of Clinical Epidemiology 65 (2012) 855e862

861

Fig. 4. Interaction plot. The y-axis gives the hazard ratio (HR) with point-wise 95% confidence intervals for x1 on mortality for varying values of x2
(on the left) and for x2 for varying values of x1 (on the right).

dichotomized, however, there is not only the chance for

a conservative error (i.e., power loss), but the probability
of false positive findings might also be increased dramatically [4]. In situations, where both explanatory variables
are correlated with each other, but the partial correlation between one variable and the outcome conditioning on the
other continuous variable is close to zero, dichotomization
leads to bias: the effect of this variable, which would not
be associated in a regression model including both variables
continuously, is then overestimated. Spurious significant interaction can also occur in such a situation. An additional
nonlinear effect in one of the explanatory variables creates
the illusion of interaction after dichotomization.
The extent of bias worsens the higher the continuous
variables are correlated. In epidemiological studies, there
is hardly a situation, where two explanatory variables in
one model are not correlated with each other. These findings are also true for categorizations with more than two
groups [14] or if only one of the continuous variables is categorized. This was shown for building categorized age
groups [15], an often-observed analysis in epidemiological
studies. The more groups are chosen for categorization,
however, the less bias is likely to occur.
As there are graphical methods available for interpretation purposes (see Figs. 1e4) categorization of continuous
data is not necessary, anyway [16]. Another reason for categorizations might be to detect nonlinear effects. To check
for nonlinearity, however, many groups are needed. Nevertheless, the real underlying structure can be missed. Therefore, flexible smoothing techniques should be preferred for
this purpose. Linear and logistic regression models can be
generalized to generalized additive models incorporating
smooth functions in the R-package mgcv [17]. For interaction effects, tensor products of smooth functions can be
constructed [18] leading to interaction surfaces reflecting
the nonlinear relationship. For survival models, this technique is, to our knowledge, not implemented in common
statistics programs. However, this can be done using Bayesian adaptive regression in the program BayesX [19], approximating the surface by a tensor product of Bayesian

P-Splines [20]. These methods, however, are statistically

more complex and this may limit their widespread use in
clinical epidemiology. There is also the problem of overinterpretation and uncertainties introduced by model instabilities. In conclusion, if the linearity assumption holds, the
highest statistical efficacy can be reached by modeling interaction terms continuously.
4.2. Understanding and interpreting interaction effects
As already stated, regression models including multiplicative interaction terms cannot be interpreted usefully just
by looking at the typical output table as in Table 1. Mean
centering of the variables has been proposed to make this
interpretation easier. This does not alleviate the problem
that the effect estimate and P-value only refers to one specific value of the other interacting variable, which is the
mean value in that case. Mean centering has also been said
to reduce the multicollinearity problem in interaction
models. However, centering is just an algebraic transformation and alters nothing important, not even multicollinearity
issues [12,21]. It does result in different coefficients and
standard errors, but the only difference is the following:
Not centering means: b1 is the effect of x1, when x2 is zero.
Centering means: b1 is the effect of x1, when x2 is at its
mean. The latter does also not help with the interpretation
of interaction effects. In our example, the conditional effect
for x2 was simulated to be 0 and estimated as |0.02. This is
the effect for x1 equal to 0. As this is also the mean of x1,
this is equal to the conditional effect after mean centering.
However, Fig. 4 and the output of our proposed functions
reveal that the effect of x2 is significantly different from
0 for 98.1% of the observed data. In conclusion, presenting
such a table as Table 1 at all is not very useful for interpretation purposes, no matter if centering was performed or
not. One needs methods spanning the complete range of
values for proper interpretation of effects on variables,
which are involved in interaction effects.
For this purpose, we propose different graphical visualization methods, which can be used for GLMs or Cox

862

C. Lamina et al. / Journal of Clinical Epidemiology 65 (2012) 855e862

models. The twisted plane in Fig. 1 shows how the effects

of both interacting variables change with changing values
of the other respective variable. However, there is a danger
of overinterpretation at the tails of the distributions. Effects
that can only be seen for very rare high or low values might
be overrated. By cutting out slices out of this twisted plane
for specific values (Fig. 3) together with its confidence
bands, one can see, if the interaction is only triggered by
a small percentage of the data at the tails of the distribution
or if it affects most of the observations. Conditional effects
with their corresponding P-values or confidence intervals
can also be given out for a specified range of values (function getRegion or getRegionCox, see Fig. 4). The interaction plots exemplified in Fig. 4 are sufficient to show
the effects for the complete range of both constituting variables together with their point-wise confidence intervals.
Reporting the significance region and the percentage of
the observations falling into this region can help to better
judge the possible implications: Are the ranges of variables
where the interaction applies realistic and relevant? An interaction effect might not be interesting especially for clinical interpretation, if it is only driven by a small percentage
of observations.
5. Conclusion
As the interpretation of interaction effects between the
two continuous variables is not straightforward, continuous
variables are categorized in most clinical epidemiological
publications leading to a loss of power and information.
We give an overview over common misconceptions in the
interpretation of interaction effects and propose several
graphical visualization methods to ease this interpretation.
Together with the calculation of significance regions, these
plots can give researchers and clinicians the information
needed to better assess the clinical relevance and implications of interaction effects.
Appendix
Supplementary material
Supplementary data related to this article can be found
online at doi:10.1016/j.jclinepi.2012.02.013.

References
[1] Barthel FM-S, Royston P. Graphical representation of interactions.
Stata J 2006;6:348e63.
[2] Turner EL, Dobson JE, Pocock SJ. Categorisation of continuous risk
factors in epidemiological publications: a survey of current practice.
Epidemiol Perspect Innov 2010;7:9.
[3] Royston P, Altman DG, Sauerbrei W. Dichotomizing continuous predictors in multiple regression: a bad idea. Stat Med 2006;25:127e41.
[4] Maxwell SE, Delaney HD. Bivariate median splits and spurious statistical significance. Psychol Bull 1993;113:181e90.
[5] Aiken LS, West SG. Multiple regression: testing and interpreting interactions. Thousand Oaks, CA: SAGE Publications; 1991.
[6] Postorino M, Marino C, Tripepi G, Zoccali C. Abdominal obesity
modifies the risk of hypertriglyceridemia for all-cause and cardiovascular mortality in hemodialysis patients. Kidney Int 2011;79:765e72.
[7] Zoccali C, Postorino M, Marino C, Pizzini P, Cutrupi S, Tripepi G.
Waist circumference modifies the relationship between the adipose
tissue cytokines leptin and adiponectin and all-cause and cardiovascular mortality in haemodialysis patients. J Intern Med 2011;
269:172e81.
[8] Rothman KJ. Epidemiology; an introduction. New York, NY: Oxford
University Press; 2002.
[9] de Mutsert R, Jager KJ, Zoccali C, Dekker FW. The effect of joint
exposures: examining the presence of interaction. Kidney Int
2009;75:677e81.
[10] Therneau TM, Grambsch PM. Modeling survival datadextending
the Cox model. New York, NY: Springer; 2000.
[11] Bender R, Augustin T, Blettner M. Generating survival times to simulate Cox proportional hazards models. Stat Med 2005;24:1713e23.
[12] Brambor T, Clark WR, Golder M. Understanding interaction models:
improving empirical analyses. Polit Anal 2006;14:63e82.
[13] Zitt E, Lamina C, Sturm G, Knoll F, Lins F, Freistatter O, et al. Interaction of time-varying albumin and phosphorus on mortality in incident dialysis patients. Clin J Am Soc Nephrol 2011;6:2650e6.
[14] Taylor JMG, Yu MG. Bias and efficiency loss due to categorizing an
explanatory variable. J Multivar Anal 2002;83:248e63.
[15] Chen H, Cohen P, Chen S. Biased odds ratios from dichotomization
of age. Stat Med 2007;26:3487e97.
[16] Altman DG. Categorizing continuous variables. In: Armitage P,
Colton T, editors. Encyclopedia of biostatistics. Chichester, UK: John
Wiley and Sons; 1998:563e7.
[17] Wood SN. Generalized additive models: an introduction with R. Boca
Raton, FL: Chapmann & Hall; 2006.
[18] Wood SN. Thin plate regression splines. J R Stat Soc B 2003;
65:95e114.
[19] Brezger A, Lang S. Generalized structured additive regression based
on Bayesian P-splines. Comput Stat Data Anal 2006;50:967e91.
[20] Fahrmeir L, Kneib T, Lang S. Regression: Modelle, Methoden und
Anwendungen. Berlin, Germany: Springer; 2007.
[21] Echembadi R, Hess JD. Mean-centering does not alleviate collinearity problems in moderated multiple regression models. Market Sci
2007;26:438e45.