Including interaction terms with total immunoglobulins
Elias Eythorsson
January 2nd, 2024
1 Interaction between isotype and total IgG, IgA and IgM
This script assumes that you have already run mgus_bm_prediction_20231125. This analysis has been completed in response to the revision request by the Annals of Internal Medicine. Members of the editorial team Drs. Wee and Popat had previously requested that we simplify the original model such that total IgG, IgA, and IgM were modeled as linear variables instead of restricted cubic splines with four-knots and we have agreed to this. Therefore, this script will explore possible interaction effects between MGUS isotype and total IgG, IgA and IgM modeled as linear variables and compare this to a model without these interactions. We have completed a similar script that evaluates interaction effects in the original non-linear model and have shown that this lead to an unacceptably high number of predictor parameters that was not justifiable.
First we fit the newly hypothesized ‘full’ model that includes interactions between total IgG, IgA and IgM and MGUS isotype.
fit_pred_inter <- lrm(outcome ~ isotype + igg + iga + igm +
rcs(m_protein, 4) + rcs(flc_ratio, 4) +
igg * isotype + iga * isotype +
igm * isotype, data = df_pred,
x = TRUE, y = TRUE)As seen below the model likelihood ratio chi-squared statistic is 406.00 on 21 degrees of freedom, with an R squared of 0.361 and a concordance index of 0.768.
fit_pred_inter## Logistic Regression Model
##
## lrm(formula = outcome ~ isotype + igg + iga + igm + rcs(m_protein,
## 4) + rcs(flc_ratio, 4) + igg * isotype + iga * isotype +
## igm * isotype, data = df_pred, x = TRUE, y = TRUE)
##
##
## Frequencies of Responses
##
## 0-4% 5-9% 10-14% >15%
## 455 430 91 67
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 1043 LR chi2 406.00 R2 0.361 C 0.768
## max |deriv| 4e-08 d.f. 21 R2(21,1043)0.309 Dxy 0.537
## Pr(> chi2) <0.0001 R2(21,882.4)0.354 gamma 0.537
## Brier 0.084 tau-a 0.337
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=5-9% 0.3553 0.7992 0.44 0.6566
## y>=10-14% -2.3223 0.8043 -2.89 0.0039
## y>=>15% -3.7255 0.8132 -4.58 <0.0001
## isotype=IgA 2.2009 1.0047 2.19 0.0285
## isotype=IgG -0.4903 0.8193 -0.60 0.5496
## isotype=Light-chain 2.3445 0.9892 2.37 0.0178
## igg 0.0807 0.0568 1.42 0.1550
## iga 0.3766 0.1104 3.41 0.0006
## igm -0.1041 0.0533 -1.95 0.0507
## m_protein -0.7044 0.1816 -3.88 0.0001
## m_protein' 47.9214 10.9594 4.37 <0.0001
## m_protein'' -60.3934 13.9100 -4.34 <0.0001
## flc_ratio -1.9265 0.3621 -5.32 <0.0001
## flc_ratio' 35.4470 6.6334 5.34 <0.0001
## flc_ratio'' -83.5600 15.7193 -5.32 <0.0001
## isotype=IgA * igg -0.2013 0.0850 -2.37 0.0178
## isotype=IgG * igg 0.0381 0.0619 0.62 0.5382
## isotype=Light-chain * igg -0.1124 0.0777 -1.45 0.1482
## isotype=IgA * iga -0.0674 0.1310 -0.51 0.6070
## isotype=IgG * iga -0.1084 0.1389 -0.78 0.4354
## isotype=Light-chain * iga -0.4420 0.1945 -2.27 0.0231
## isotype=IgA * igm -0.3324 0.3076 -1.08 0.2799
## isotype=IgG * igm -0.1632 0.1616 -1.01 0.3126
## isotype=Light-chain * igm -0.5762 0.4223 -1.36 0.1724
## Importantly for all further comparisons, both models must be shown to be equally well calibrated, as we would heavily favor a calibrated over a poorly calibrated model. What constitutes ‘equally well calibrated’ is subjective, and is best judged by examining calibration plots for both the interaction model and the non-interaction models. In the case of the interaction model with linear IgG, IgA and IgM (this script), calibration is excellent.
plot(calibrate(fit_pred_inter, B = 200))
##
## n=1043 Mean absolute error=0.011 Mean squared error=0.00016
## 0.9 Quantile of absolute error=0.017Next we fit the non-interaction model that includes total IgG, IgA and IgM as linear variables.
fit_pred_linear <- lrm(outcome ~ igg + iga + igm + rcs(m_protein, 4) +
rcs(flc_ratio, 4) + isotype, data = df_pred, x = TRUE, y = TRUE)As seen below the model likelihood ratio chi-square statistic is 378.19 on 12 degrees of freedom, with an R squared of 0.341 and a concordance index of 0.760.
fit_pred_linear## Logistic Regression Model
##
## lrm(formula = outcome ~ igg + iga + igm + rcs(m_protein, 4) +
## rcs(flc_ratio, 4) + isotype, data = df_pred, x = TRUE, y = TRUE)
##
##
## Frequencies of Responses
##
## 0-4% 5-9% 10-14% >15%
## 455 430 91 67
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 1043 LR chi2 378.19 R2 0.341 C 0.760
## max |deriv| 4e-08 d.f. 12 R2(12,1043)0.296 Dxy 0.519
## Pr(> chi2) <0.0001 R2(12,882.4)0.340 gamma 0.519
## Brier 0.086 tau-a 0.326
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=5-9% 1.4481 0.4232 3.42 0.0006
## y>=10-14% -1.1724 0.4231 -2.77 0.0056
## y>=>15% -2.5540 0.4380 -5.83 <0.0001
## igg 0.0130 0.0230 0.56 0.5722
## iga 0.2854 0.0451 6.32 <0.0001
## igm -0.1745 0.0464 -3.76 0.0002
## m_protein -0.6690 0.1797 -3.72 0.0002
## m_protein' 48.8852 10.8537 4.50 <0.0001
## m_protein'' -61.5660 13.7741 -4.47 <0.0001
## flc_ratio -2.1699 0.3533 -6.14 <0.0001
## flc_ratio' 39.8168 6.4875 6.14 <0.0001
## flc_ratio'' -93.8202 15.3791 -6.10 <0.0001
## isotype=IgA -0.1051 0.2525 -0.42 0.6771
## isotype=IgG -0.4976 0.2008 -2.48 0.0132
## isotype=Light-chain -0.4562 0.3211 -1.42 0.1554
## The calibration of the non-interaction model.
plot(calibrate(fit_pred, B = 200))
##
## n=1043 Mean absolute error=0.015 Mean squared error=0.00026
## 0.9 Quantile of absolute error=0.0242 Effect differences for the interaction model
ggplot(
Predict(
fit_pred_inter,
kint = 2,
igg, isotype,
ref.zero = TRUE, fun = exp
), adj.subtitle = FALSE) +
geom_hline(yintercept = 1, lty = 2) +
geom_vline(xintercept = 6.1, lty = 2) +
geom_vline(xintercept = 15.7, lty = 2) +
geom_rug(data = df_pred, aes(x = igg, y = NULL)) +
scale_y_continuous(breaks = c(0.1, 0.25, 0.5, 0.75, 1, 1.5, 2, 4, 8, 16, 32)) +
scale_x_continuous(limits = c(5, 22)) +
coord_trans(y = "log10", ylim = c(0.1, 32)) +
labs(x = "IgG (g/L)", y = "OR of 10% BMPC or more", title = "Total IgG") +
theme_bw() +
theme(legend.position = "bottom")
ggplot(
Predict(
fit_pred_inter,
kint = 2,
igg, isotype,
fun = plogis
), adj.subtitle = FALSE) +
coord_cartesian(ylim = c(0, 0.25)) +
geom_vline(xintercept = 6.1, lty = 2) +
geom_vline(xintercept = 15.7, lty = 2) +
geom_rug(data = df_pred, aes(x = igg, y = NULL)) +
scale_y_continuous(labels = scales::percent) +
scale_x_continuous(limits = c(5, 22)) +
labs(x = "IgG (g/L)", y = "Proportion (%) having \n≥10% BMPC", title = "Total IgG") +
theme_bw() +
theme(legend.position = "bottom")
ggplot(
Predict(
fit_pred_inter,
kint = 2,
igm, isotype,
ref.zero = TRUE, fun = exp
)) +
geom_hline(yintercept = 1, lty = 2) +
geom_vline(xintercept = 0.4, lty = 2) +
geom_vline(xintercept = 2.3, lty = 2) +
geom_rug(data = df_pred, aes(x = igm, y = NULL)) +
scale_y_continuous(breaks = c(0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 1, 1.5, 2, 4, 8, 16, 32)) +
scale_x_continuous(limits = c(NA, 8)) +
coord_trans(y = "log10", ylim = c(0.01, 32)) +
labs(x = "IgM (g/L)", y = "OR of 10% BMPC or more", title = "Total IgM") +
theme_bw() +
theme(legend.position = "bottom")
ggplot(
Predict(
fit_pred_inter,
kint = 2,
igm, isotype,
fun = plogis
)) +
coord_cartesian(ylim = c(0, 0.1)) +
geom_vline(xintercept = 0.4, lty = 2) +
geom_vline(xintercept = 2.3, lty = 2) +
geom_rug(data = df_pred, aes(x = igm, y = NULL)) +
scale_y_continuous(labels = scales::percent) +
scale_x_continuous(limits = c(NA, 8)) +
labs(x = "IgM (g/L)", y = "Proportion (%) having \n≥10% BMPC", title = "Total IgM") +
theme_bw() +
theme(legend.position = "bottom")
ggplot(
Predict(
fit_pred_inter,
kint = 2,
iga, isotype,
ref.zero = TRUE, fun = exp
), adj.subtitle = FALSE) +
geom_hline(yintercept = 1, lty = 2) +
geom_vline(xintercept = 0.7, lty = 2) +
geom_vline(xintercept = 3.7, lty = 2) +
geom_rug(data = df_pred, aes(x = iga, y = NULL)) +
scale_y_continuous(breaks = c(0.1, 0.25, 0.5, 0.75, 1, 1.5, 2, 4, 8, 16, 32)) +
scale_x_continuous(limits = c(NA, 10)) +
coord_trans(y = "log10", ylim = c(0.1, 32)) +
labs(x = "IgA (g/L)", y = "OR of 10% BMPC or more", title = "Total IgA") +
theme_bw() +
theme(legend.position = "bottom")
ggplot(
Predict(
fit_pred_inter,
kint = 2,
iga, isotype,
fun = plogis
), adj.subtitle = FALSE) +
coord_cartesian(ylim = c(0, 0.2)) +
geom_vline(xintercept = 0.4, lty = 2) +
geom_vline(xintercept = 2.3, lty = 2) +
geom_rug(data = df_pred, aes(x = iga, y = NULL)) +
scale_y_continuous(labels = scales::percent) +
scale_x_continuous(limits = c(NA, 6)) +
labs(x = "IgA (g/L)", y = "Proportion (%) having \n≥10% BMPC", title = "Total IgA") +
theme_bw() +
theme(legend.position = "bottom")
3 Likelihood ratio test
Because the non-interaction model is nested within the interaction model, we can use the likelihood ratio test to test the null hypothesis that both models fit the data equally well, and therefore that the interaction terms are redundant. As seen below, we find some evidence against this null hypothesis.
lrtest(fit_pred_linear, fit_pred_inter)##
## Model 1: outcome ~ igg + iga + igm + rcs(m_protein, 4) + rcs(flc_ratio,
## 4) + isotype
## Model 2: outcome ~ isotype + igg + iga + igm + rcs(m_protein, 4) + rcs(flc_ratio,
## 4) + igg * isotype + iga * isotype + igm * isotype
##
## L.R. Chisq d.f. P
## 27.811232694 9.000000000 0.001025651We can quantify the degree to which the interaction model fits the data better, by calculating the adequacy of the non-interaction model compared to the non-linear model using the models’ likelihood ratio Chi-squared statistics.
pre_test_lr <- fit_pred_linear$stats["Model L.R."] # non-interaction likelihood chi-squared statistic
post_test_lr <- fit_pred_inter$stats["Model L.R."] # interactoni likelihood chi-suqared statisticBelow is a measure of the adequacy of the linear model compared to the non-linear model, 0.931 out of a maximum of 1
add_base <- pre_test_lr/post_test_lr
add_base## Model L.R.
## 0.9314997In other words, the fraction of new information provided by interaction model by this metric is 0.069, which is substantial.
frac_new <- 1 - add_base
frac_new## Model L.R.
## 0.06850028The contribution of the interaction components of total IgG, IgA and IgM (or alternatively MGUS isotype) can be shown for each immunoglobulin individually, with the overall interaction likelihood ratio Chi-squared statistic of 27.1 on 9 degrees of freedom approximating the overall likelihood ratio test above.
anova(fit_pred_inter)## Wald Statistics Response: outcome
##
## Factor Chi-Square d.f. P
## isotype (Factor+Higher Order Factors) 34.91 12 0.0005
## All Interactions 27.10 9 0.0013
## igg (Factor+Higher Order Factors) 15.16 4 0.0044
## All Interactions 12.87 3 0.0049
## iga (Factor+Higher Order Factors) 40.85 4 <.0001
## All Interactions 5.55 3 0.1360
## igm (Factor+Higher Order Factors) 11.91 4 0.0181
## All Interactions 3.70 3 0.2960
## m_protein 47.72 3 <.0001
## Nonlinear 24.20 2 <.0001
## flc_ratio 50.26 3 <.0001
## Nonlinear 29.79 2 <.0001
## isotype * igg (Factor+Higher Order Factors) 12.87 3 0.0049
## isotype * iga (Factor+Higher Order Factors) 5.55 3 0.1360
## isotype * igm (Factor+Higher Order Factors) 3.70 3 0.2960
## TOTAL NONLINEAR 58.26 4 <.0001
## TOTAL INTERACTION 27.10 9 0.0013
## TOTAL NONLINEAR + INTERACTION 96.75 13 <.0001
## TOTAL 317.69 21 <.00014 Akaike’s information criterion
The question remains whether the potentially increased fit of the interaction model to the data compared to the non-interaction model is worthwhile, given its added complexity. Akaike’s information criterion provides a method for penalizing the log likelihood achieved by a given model for its complexity to obtain a more unbiased assessment of the model’s worth. The penalty is to add twice the number ofparameters to the −2 log likelihood. On this metric, the model with the lower AIC should generally be favored. In this case, the interaction model has a considerably lower AIC than the model that contains total IgG, IgA and IgM as linear variables.
AIC(fit_pred_linear)## [1] 1980.479AIC(fit_pred_inter)## [1] 1970.6685 Comparing variances
Next we compare the variances of the predicted risks from the non-interaction and interaction models. In the prediction setting, increased variance in risk estimates is generally desirable. As shown below, the variance is almost identical between the linear and non-linear models and the high-resolution histogram of predicted risks is visually identical.
pre <- predict(fit_pred_linear, type='fitted')
pre <- pre[,2]
post <- predict(fit_pred_inter, type='fitted')
post <- post[,2]
par(mgp=c(4,1,0), mar=c(6,5,2,2))
xlab <- c(
paste0("Non-interaction\nvariance =", round(var(pre), 3)),
paste0("Interaction\nvariance=",round(var(post), 3))
)
histbackback(pre, post, brks=seq(0, 1, by=0.005), xlab=xlab, ylab= "Proportion (%) having \n≥10% BMPC")
We show that the relative explained variation of predicted risk of the non-interaction model compared to the interaction model is 0.984 (theoretically out of a maximum of 1)
var_pre <- round(var(pre), 3) #linear
var_post <- round(var(post), 3) #non-linear
relativ_expl <- var_pre/var_post
relativ_expl## [1] 0.974359In other words, the fraction of new information contained in the interaction model by this metric is only 0.0255
frac_new2 <- 1 - relativ_expl
frac_new2 ## [1] 0.02564103We can show that the predicted risks between the non-interaction and interaction models differ most at moderate predicted risks.
plot(pre, post-pre, xlab='Predicted risk (non-interaction)',
ylab='Difference in predicted risks (interaction - non-interaction)', pch=46)
lo <- Rq(post-pre ~ rcs(pre, 7), tau=0.1) # 0.1 quantile
hi <- Rq(post-pre ~ rcs(pre, 7), tau=0.9) # 0.9 quantile
at <- seq(0, 1, length=200)
lines(at, Predict(lo, pre=at)$yhat, col='red', lwd=1.5)
lines(at, Predict(hi, pre=at)$yhat, col='red', lwd=1.5)
abline(h=0, col='red', lty = 2)
6 Decision curve analysis
Finally, as this is a prediction model that we intend for clinical use, we should be satisfied that the decision curves associated with both models are similar. If not, we should heavily favor the model that leads to increased net-benefit over plausible risk thresholds.
df_interaction <- df_pred # temporary dataframe only used for this script
Function(fit_pred_linear) # we extract the linear predictor from the non-interaction model first## function(igg = 10.21,iga = 2.167,igm = 0.718,m_protein = 1.91,flc_ratio = 1.237794,isotype = "IgG") {+0.012987234*igg+0.28536845*iga-0.17450545*igm-0.66897437* m_protein+0.47127836*pmax(m_protein-0.481261,0)^3-0.59352797*pmax(m_protein-0.881164,0)^3+0.12313286*pmax(m_protein-2.48194,0)^3-0.00088325512*pmax(m_protein-10.665996,0)^3-2.1699146* flc_ratio+0.76041857*pmax(flc_ratio-0.2041965,0)^3-1.7917724*pmax(flc_ratio-1.0186109,0)^3+1.0353426*pmax(flc_ratio-1.6415064,0)^3-0.0039887145*pmax(flc_ratio-7.440332,0)^3-0.10513226*(isotype=="IgA")-0.49763544*(isotype=="IgG")-0.45621863*(isotype=="Light-chain") }
## <environment: 0x7fe9ec189ac0>df_interaction$linear_pred <- with(
df_linear_nonlinear,
expit(
-1.17235498 +0.012987234 * igg + 0.28536845 * iga - 0.17450545 * igm - 0.66897437 * m_protein +
0.47127836 * pmax(m_protein - 0.481261, 0) ^ 3 - 0.59352797 * pmax(m_protein - 0.881164, 0) ^ 3
+ 0.12313286 * pmax(m_protein - 2.48194, 0) ^ 3 - 0.00088325512 *
pmax(m_protein - 10.665996, 0) ^ 3 - 2.1699146 * flc_ratio + 0.76041857 *
pmax(flc_ratio - 0.2041965, 0) ^ 3 - 1.7917724 * pmax(flc_ratio - 1.0186109, 0) ^
3 + 1.0353426 * pmax(flc_ratio - 1.6415064, 0) ^ 3 - 0.0039887145 * pmax(flc_ratio - 7.440332, 0) ^ 3
- 0.10513226 * (isotype == "IgA") - 0.49763544 * (isotype == "IgG")
- 0.45621863 * (isotype == "Light-chain") )
)
Function(fit_pred_inter) # next we extract the linear predictor from the linear model## function(isotype = "IgG",igg = 10.21,iga = 2.167,igm = 0.718,m_protein = 1.91,flc_ratio = 1.237794) {+2.2008812*(isotype=="IgA")-0.49027132*(isotype=="IgG")+2.3445491*(isotype=="Light-chain")+0.080719829*igg+0.37659232*iga-0.10413849*igm-0.70441891* m_protein+0.46198715*pmax(m_protein-0.481261,0)^3-0.58222351*pmax(m_protein-0.881164,0)^3+0.12117983*pmax(m_protein-2.48194,0)^3-0.00094347138*pmax(m_protein-10.665996,0)^3-1.9264969* flc_ratio+0.67696428*pmax(flc_ratio-0.2041965,0)^3-1.5958236*pmax(flc_ratio-1.0186109,0)^3+0.92248483*pmax(flc_ratio-1.6415064,0)^3-0.0036255376*pmax(flc_ratio-7.440332,0)^3+igg*(-0.20125818*(isotype=="IgA")+0.03807887*(isotype=="IgG")-0.11239517*(isotype=="Light-chain"))+iga*(-0.067380239*(isotype=="IgA")-0.10837871*(isotype=="IgG")-0.44202957*(isotype=="Light-chain"))+igm*(-0.33237601*(isotype=="IgA")-0.16323079*(isotype=="IgG")-0.57616915*(isotype=="Light-chain")) }
## <environment: 0x7fea029c2ed8>df_interaction$inter_pred <- with(
df_interaction,
expit(
-2.3223 +2.2008812 * (isotype == "IgA") - 0.49027132 * (isotype == "IgG") + 2.3445491 *
(isotype == "Light-chain") + 0.080719829 * igg + 0.37659232 * iga - 0.10413849 *
igm - 0.70441891 * m_protein + 0.46198715 * pmax(m_protein - 0.481261, 0) ^
3 - 0.58222351 * pmax(m_protein - 0.881164, 0) ^ 3 + 0.12117983 * pmax(m_protein -2.48194, 0) ^ 3 - 0.00094347138 * pmax(m_protein - 10.665996, 0) ^ 3 - 1.9264969 * flc_ratio +
0.67696428 * pmax(flc_ratio - 0.2041965, 0) ^ 3 - 1.5958236 * pmax(flc_ratio -1.0186109, 0) ^ 3 + 0.92248483 * pmax(flc_ratio - 1.6415064, 0) ^ 3 - 0.0036255376 *
pmax(flc_ratio - 7.440332, 0) ^ 3 + igg * (-0.20125818 * (isotype == "IgA") + 0.03807887 * (isotype == "IgG") - 0.11239517 * (isotype == "Light-chain")
) + iga * (-0.067380239 * (isotype == "IgA") - 0.10837871 * (isotype == "IgG") - 0.44202957 *(isotype == "Light-chain")
) + igm * (-0.33237601 * (isotype == "IgA") - 0.16323079 * (isotype == "IgG") - 0.57616915 *(isotype == "Light-chain")
)
)
)We create the decision curve analysis models using the dca pacakge. We examine the net-benefit of the non-interaction and interaction models in a opt-out strategy, in which net benefit is expressed in the relative proportion of true negatives. As demonstrated by the resulting decision curve analysis figure, there is no perceivable difference in net benefit between the two models in this setting.
noninteraction_model <- decision_curve(perc_10 ~ linear_pred,
data = df_interaction,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 500)
interaction_model <- decision_curve(perc_10 ~ inter_pred,
data = df_interaction,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 500)
noninteraction_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
union_all(
interaction_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
filter(model == "perc_10 ~ inter_pred")
) %>%
mutate(model = factor(model, levels = c("perc_10 ~ linear_pred", "perc_10 ~ inter_pred", "All", "None"))) %>%
ggplot(aes(x = thresholds, y = NB, color = model, fill = model)) +
geom_line(size = 1.1, position = position_dodge(width = 0.005)) +
geom_vline(aes(xintercept = 0.05), lty = 2) +
geom_vline(aes(xintercept = 0.10), lty = 2) +
scale_color_manual(
name = "Follow-up strategy:",
labels = c(
"Defer based on non-interaction model",
"Defer based on interaction model",
"Obtain sample from all",
"Obtain sample from none"
),
values = c(nejm_palette[c(2, 3)], "black", nejm_palette[c(1)])
) +
scale_x_continuous(breaks = seq(0, 1, by = 0.02), labels = scales::percent_format(accuracy = 1)) +
coord_cartesian(xlim = c(0, 0.25), ylim = c(-0.2, 0.7)) +
labs(x = NULL, y = "Net Benefit") +
theme_bw() +
theme(
legend.title = element_blank(),
legend.background = element_rect(fill = NA),
panel.grid.minor.x = element_blank(),
legend.position = c(0.2, 0.82)
)
7 Summary
We have shown that the model that includes an interaction term between MGUS isotype and total IgG, IgA and IgM modeled as linear variables: i) is similarly calibrated ii) fits the data significantly better on the likelihood ratio test at an alpha of 0.05, providing a roughly 6.9% increase in the fraction of new information on this metric compared to the non-interaction model model. iii) has a considerably lower Akaike’s information criterion than does the non-interaction model iv) slightly increases the variance of predicted risk, providing a roughly 2.5% increase in the fraction of new information on this metric compared to the non-interaction model v) has similar net-benefit on decision curve analysis.
This may at first glance suggest that an interaction model is superior. However, we do not feel this to be the case. Our argument is based on the linear interaction model being computationally, statistically and practically much less intuitive than the original pre-specified, fully non-linear model, without providing any measurable benefits compared to that model. We show this below. We have previously conceded to simplify our pre-specified model to include total IgG, IgA and IgM as linear terms. Any of the above arguments for the linear interaction model should apply to the pre-specified fully non-linear model, which has the added benefit of being fully pre-specified. Comparison of the linear interaction model and the fully pre-specified model ## Likelihood ratio test No evidence for the interaction model fitting the data better than the original.
lrtest(fit_pred_inter, fit_pred)##
## Model 1: outcome ~ isotype + igg + iga + igm + rcs(m_protein, 4) + rcs(flc_ratio,
## 4) + igg * isotype + iga * isotype + igm * isotype
## Model 2: outcome ~ igg + iga + igm + rcs(m_protein, parms = c(0.5, 1,
## 2.5, 10.5)) + rcs(flc_ratio, parms = c(0.2, 1, 1.6, 7.5)) +
## isotype
##
## L.R. Chisq d.f. P
## 27.425714799 9.000000000 0.001188976pre_test_lr <- fit_pred$stats["Model L.R."] # non-interaction likelihood chi-squared statistic
post_test_lr <- fit_pred_inter$stats["Model L.R."] # interactoni likelihood chi-suqared statisticBelow is a measure of the adequacy of the original model compared to the linear interaction model, 0.987 out of a maximum of 1
add_base <- pre_test_lr/post_test_lr
add_base## Model L.R.
## 0.9324493In other words, the fraction of new information provided by interaction model by this metric is 0.012.
frac_new <- 1 - add_base
frac_new## Model L.R.
## 0.067550747.1 Akaike’s information criterion
The AIC of the original model is lower.
AIC(fit_pred)## [1] 1980.093AIC(fit_pred_inter)## [1] 1970.6687.2 Comparing variances
pre <- predict(fit_pred, type='fitted')
pre <- pre[,2]
post <- predict(fit_pred_inter, type='fitted')
post <- post[,2]
par(mgp=c(4,1,0), mar=c(6,5,2,2))
xlab <- c(
paste0("Non-interaction\nvariance =", round(var(pre), 3)),
paste0("Interaction\nvariance=",round(var(post), 3))
)
histbackback(pre, post, brks=seq(0, 1, by=0.005), xlab=xlab, ylab= "Proportion (%) having \n≥10% BMPC")
We show that the relative explained variation of predicted risk of the non-interaction model compared to the interaction model is 1 (theoretically out of a maximum of 1)
var_pre <- round(var(pre), 3) #linear
var_post <- round(var(post), 3) #non-linear
relativ_expl <- var_pre/var_post
relativ_expl## [1] 0.974359In other words, the fraction of new information contained in the interaction model by this metric is 0.
frac_new2 <- 1 - relativ_expl
frac_new2 ## [1] 0.025641038 Decision curve analysis
Finally, as this is a prediction model that we intend for clinical use, we should be satisfied that the decision curves associated with both models are similar. If not, we should heavily favor the model that leads to increased net-benefit over plausible risk thresholds.
Function(fit_pred) # we extract the linear predictor from the non-interaction model first## function(igg = 10.21,iga = 2.167,igm = 0.718,m_protein = 1.91,flc_ratio = 1.237794,isotype = "IgG") {+0.013109259*igg+0.28583096*iga-0.17434261*igm-0.62628967* m_protein+0.36358352*pmax(m_protein-0.5,0)^3-0.48979428*pmax(m_protein-1,0)^3+0.12715131*pmax(m_protein-2.5,0)^3-0.00094054869*pmax(m_protein-10.5,0)^3-2.2468523* flc_ratio+0.81545513*pmax(flc_ratio-0.2,0)^3-1.9401717*pmax(flc_ratio-1,0)^3+1.1285244*pmax(flc_ratio-1.6,0)^3-0.003807771*pmax(flc_ratio-7.5,0)^3-0.10647038*(isotype=="IgA")-0.49770728*(isotype=="IgG")-0.4587627*(isotype=="Light-chain") }
## <environment: 0x7fea13c98d50>df_interaction$pred <- with(
df_linear_nonlinear,
expit(
1.5889 -0.35437745 * igg + 0.012363328 * pmax(igg - 6.111, 0) ^ 3 - 0.034311919 *
pmax(igg - 9.2, 0) ^ 3 + 0.023794703 * pmax(igg - 11.446, 0) ^ 3 - 0.0018461125 *
pmax(igg - 17.462, 0) ^ 3 - 0.13301362 * iga + 0.098986509 * pmax(iga -
0.5914, 0) ^ 3 - 0.22264315 * pmax(iga - 1.641, 0) ^ 3 + 0.13057034 * pmax(iga -
2.6581, 0) ^ 3 - 0.006913695 * pmax(iga - 5.8221, 0) ^ 3 + 1.0889587 * igm -
3.1076573 * pmax(igm - 0.2561, 0) ^ 3 + 6.3686218 * pmax(igm - 0.5751, 0) ^
3 - 3.3346563 * pmax(igm - 0.9166, 0) ^ 3 + 0.073691742 * pmax(igm - 2.5759, 0) ^
3 - 0.66018976 * m_protein + 0.44395199 * pmax(m_protein - 0.481261, 0) ^
3 - 0.55838218 * pmax(m_protein - 0.881164, 0) ^ 3 + 0.11511927 * pmax(m_protein -
2.48194, 0) ^ 3 - 0.00068906911 * pmax(m_protein - 10.665996, 0) ^ 3 - 2.1257973 * flc_ratio +
0.73647737 * pmax(flc_ratio - 0.2041965, 0) ^ 3 - 1.7333446 * pmax(flc_ratio -
1.0186109, 0) ^ 3 + 1.0005139 * pmax(flc_ratio - 1.6415064, 0) ^ 3 - 0.0036466553 *
pmax(flc_ratio - 7.440332, 0) ^ 3 - 0.27663453 * (isotype == "IgA") - 0.53499524 *
(isotype == "IgG") - 0.56401632 * (isotype == "Light-chain") )
)We create the decision curve analysis models using the dca pacakge. We examine the net-benefit of the non-interaction and interaction models in a opt-out strategy, in which net benefit is expressed in the relative proportion of true negatives. As demonstrated by the resulting decision curve analysis figure, there is no perceivable difference in net benefit between the two models in this setting.
noninteraction_model <- decision_curve(perc_10 ~ pred,
data = df_interaction,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 500)
interaction_model <- decision_curve(perc_10 ~ inter_pred,
data = df_interaction,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 500)
noninteraction_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
union_all(
interaction_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
filter(model == "perc_10 ~ inter_pred")
) %>%
mutate(model = factor(model, levels = c("perc_10 ~ pred", "perc_10 ~ inter_pred", "All", "None"))) %>%
ggplot(aes(x = thresholds, y = NB, color = model, fill = model)) +
geom_line(size = 1.1, position = position_dodge(width = 0.005)) +
geom_vline(aes(xintercept = 0.05), lty = 2) +
geom_vline(aes(xintercept = 0.10), lty = 2) +
scale_color_manual(
name = "Follow-up strategy:",
labels = c(
"Defer based on non-interaction model",
"Defer based on interaction model",
"Obtain sample from all",
"Obtain sample from none"
),
values = c(nejm_palette[c(2, 3)], "black", nejm_palette[c(1)])
) +
scale_x_continuous(breaks = seq(0, 1, by = 0.02), labels = scales::percent_format(accuracy = 1)) +
coord_cartesian(xlim = c(0, 0.25), ylim = c(-0.2, 0.7)) +
labs(x = NULL, y = "Net Benefit") +
theme_bw() +
theme(
legend.title = element_blank(),
legend.background = element_rect(fill = NA),
panel.grid.minor.x = element_blank(),
legend.position = c(0.2, 0.82)
)