Should total IgG, IgA and IgM be included in the model?
Elias Eythorsson
November 26th, 2023
1 Including 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, and will test the hypothesis that an ordinal regression containing total IgG, IgA and IgM modeled non-linearly using restricted cubic splines with four knots adds significant predictive value, compared to an ordinal model that does not include total IgG, IgA and IgM. Evidence for this will be supplied by a comparison of the pre-test variance (without total IgG, IgA and IgM) and post-test variance (non-linear immunoglobulins), a likelihood ratio test and various measures of increased explained variance, after confirming that both models are well calibrated and provide similar net-benefit on decision curve analysis.
First we fit the non-linear ‘full’ model
fit_pred <- lrm(outcome ~ rcs(igg, 4) + rcs(iga, 4) + rcs(igm, 4) + 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 401.11 on 18 degrees of freedom, with an R squared of 0.358 and a concordance index of 0.766.
fit_pred## Logistic Regression Model
##
## lrm(formula = outcome ~ rcs(igg, 4) + rcs(iga, 4) + rcs(igm,
## 4) + 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 401.11 R2 0.358 C 0.766
## max |deriv| 4e-08 d.f. 18 R2(18,1043)0.307 Dxy 0.533
## Pr(> chi2) <0.0001 R2(18,882.4)0.352 gamma 0.533
## Brier 0.085 tau-a 0.335
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=5-9% 4.2561 0.8788 4.84 <0.0001
## y>=10-14% 1.5889 0.8686 1.83 0.0674
## y>=>15% 0.1966 0.8704 0.23 0.8213
## igg -0.3544 0.0937 -3.78 0.0002
## igg' 1.5930 0.4062 3.92 <0.0001
## igg'' -4.4209 1.1912 -3.71 0.0002
## iga -0.1330 0.2827 -0.47 0.6380
## iga' 2.7083 1.7328 1.56 0.1181
## iga'' -6.0916 3.9008 -1.56 0.1184
## igm 1.0890 0.7754 1.40 0.1602
## igm' -16.7238 8.5061 -1.97 0.0493
## igm'' 34.2726 17.1540 2.00 0.0457
## m_protein -0.6602 0.1832 -3.60 0.0003
## m_protein' 46.0506 11.3215 4.07 <0.0001
## m_protein'' -57.9203 14.3889 -4.03 <0.0001
## flc_ratio -2.1258 0.3555 -5.98 <0.0001
## flc_ratio' 38.5632 6.5415 5.90 <0.0001
## flc_ratio'' -90.7608 15.5093 -5.85 <0.0001
## isotype=IgA -0.2766 0.2610 -1.06 0.2892
## isotype=IgG -0.5350 0.2038 -2.63 0.0087
## isotype=Light-chain -0.5640 0.3257 -1.73 0.0833
## Importantly for all further comparisons, both models must be shown to be equally well calibrated, as we would heavily favor a calibrated over a non-calibrated model. What constitutes ‘equally well calibrated’ is subjective, and is best judged by examining calibration plots for both models. The non-linear model is well calibrated.
plot(calibrate(fit_pred))
##
## n=1043 Mean absolute error=0.013 Mean squared error=0.00021
## 0.9 Quantile of absolute error=0.021Next we fit the non-immunoglobulin model (non-Ig). The non-Ig is a nested model of the non-linear model, in which the coefficients for total IgG, IgA, and IgM are set to zero.
fit_pred_nonig <- lrm(
formula = outcome ~ 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 314.21 on 9 degrees of freedom, with an R squared of 0.291 and a concordance index of 0.734.
fit_pred_nonig## Logistic Regression Model
##
## lrm(formula = outcome ~ 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 314.21 R2 0.291 C 0.734
## max |deriv| 2e-07 d.f. 9 R2(9,1043)0.254 Dxy 0.467
## Pr(> chi2) <0.0001 R2(9,882.4)0.292 gamma 0.468
## Brier 0.090 tau-a 0.294
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=5-9% 1.5474 0.3410 4.54 <0.0001
## y>=10-14% -0.9605 0.3395 -2.83 0.0047
## y>=>15% -2.2279 0.3551 -6.27 <0.0001
## m_protein -0.7066 0.1777 -3.98 <0.0001
## m_protein' 54.9494 10.6527 5.16 <0.0001
## m_protein'' -69.4951 13.5178 -5.14 <0.0001
## flc_ratio -1.8674 0.3381 -5.52 <0.0001
## flc_ratio' 34.2502 6.2171 5.51 <0.0001
## flc_ratio'' -80.6915 14.7444 -5.47 <0.0001
## isotype=IgA 0.8050 0.2171 3.71 0.0002
## isotype=IgG -0.3047 0.1848 -1.65 0.0991
## isotype=Light-chain 0.0037 0.3057 0.01 0.9904
## Again, we must visually be satisfied that the calibration is equally as good as the non-linear model. This is not the case, which is a strong arguement against removing total IgG, IgA and IgM from the model.
plot(calibrate(fit_pred_nonig))
##
## n=1043 Mean absolute error=0.021 Mean squared error=0.00068
## 0.9 Quantile of absolute error=0.0372 Likelihood ratio test
Because the non-Ig model is nested within the non-linear model, we can use the likelihood ratio test to test the null hypothesis that the non-Ig and non-linear models fit the data equally well, e.g. that IgG, IgA and IgM are redundant. As seen below, there is fairly strong evidence against the null hypothesis, and we may conclude that the non-linear model fits the data better.
lrtest(fit_pred_nonig, fit_pred)##
## Model 1: outcome ~ rcs(m_protein, 4) + rcs(flc_ratio, 4) + isotype
## Model 2: outcome ~ rcs(igg, 4) + rcs(iga, 4) + rcs(igm, 4) + rcs(m_protein,
## 4) + rcs(flc_ratio, 4) + isotype
##
## L.R. Chisq d.f. P
## 8.690152e+01 9.000000e+00 6.772360e-15The degree to which the non-linear model fits the data better can be quantified by calculating the adequacy of the non-Ig model compared to the non-linear model using the models’ likelihood ratio Chi-squared statistics.
pre_test_lr <- fit_pred_nonig$stats["Model L.R."] # non-Ig likelihood chi-squared statistic
post_test_lr <- fit_pred$stats["Model L.R."] # non-linear likelihood chi-suqared statisticBelow is a measure of the adequacy of the non-Ig model compared to the non-linear model, 0.783 out of a maximum of 1
add_base <- pre_test_lr/post_test_lr
add_base## Model L.R.
## 0.7833485In other words, the fraction of new information provided by non-linear model by this metric is 0.217, which is substantial.
frac_new <- 1 - add_base
frac_new## Model L.R.
## 0.2166515The contribution of total IgG, IgA and IgM can can be shown for each immunoglobulin individually, and is fairly significant.
anova(fit_pred)## Wald Statistics Response: outcome
##
## Factor Chi-Square d.f. P
## igg 17.23 3 0.0006
## Nonlinear 16.82 2 0.0002
## iga 40.20 3 <.0001
## Nonlinear 2.45 2 0.2945
## igm 21.20 3 0.0001
## Nonlinear 4.67 2 0.0968
## m_protein 69.79 3 <.0001
## Nonlinear 25.48 2 <.0001
## flc_ratio 62.21 3 <.0001
## Nonlinear 37.76 2 <.0001
## isotype 7.60 3 0.0552
## TOTAL NONLINEAR 92.93 10 <.0001
## TOTAL 320.66 18 <.0001plot(anova(fit_pred))
3 Akaike’s information criterion
The question remains whether the increased fit of the non-linear model to the data compared to the non-Ig 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 of parameters to the −2 log likelihood. On this metric, the model with the lower AIC should generally be favored. In this case, the non-linear model has a considerably lower AIC.
AIC(fit_pred) # AIC of the non-linear model## [1] 1969.557AIC(fit_pred_nonig) # AIC of the non-Ig model## [1] 2038.4594 Comparing variances
Next we compare the variances of the predicted risks from the non-Ig and non-linear models. In the prediction setting, increased variance in risk estimates is generally desirable. As shown below, the variance is visibly slightly higher in the model containing total IgG, IgA and IgM compared to the model without them.
pre <- predict(fit_pred_nonig, type='fitted')
pre <- pre[,2]
post <- predict(fit_pred, type='fitted')
post <- post[,2]
par(mgp=c(4,1,0), mar=c(6,5,2,2))
xlab <- c(
paste0("Non-immunoglobulin\nvariance =", round(var(pre), 3)),
paste0("Non-linear\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-Ig model compared to the non-linear model is 0.795 (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.7948718In other words, the fraction of new information contained in the non-linear model by this metric is 0.205, which is again substantial.
frac_new2 <- 1 - relativ_expl
frac_new2 ## [1] 0.2051282We can show that the predicted risks between the non-Ig and non-linear models differ fairly significantly, and that the way the differ is not symmetrical.
plot(pre, post-pre, xlab='Predicted risk (non-immunoglobulin)',
ylab='Difference in predicted risks (non-linear - non-immunoglobulin)', 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)
5 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_nonig <- df_pred # temporary dataframe only used for this script
Function(fit_pred) # we extract the linear predictor from the non-linear model first## function(igg = 10.21,iga = 2.167,igm = 0.718,m_protein = 1.91,flc_ratio = 1.237794,isotype = "IgG") {-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") }
## <environment: 0x7fe9cc99db28>df_nonig$nonlinear_pred <- with(
df_nonig,
expit(
1.5888896 -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")
)
)
Function(fit_pred_nonig) # next we extract the linear predictor from the non-immunoglobulin model## function(m_protein = 1.91,flc_ratio = 1.237794,isotype = "IgG") {-0.7066287* m_protein+0.52974103*pmax(m_protein-0.481261,0)^3-0.66996903*pmax(m_protein-0.881164,0)^3+0.14177108*pmax(m_protein-2.48194,0)^3-0.0015430737*pmax(m_protein-10.665996,0)^3-1.8673982* flc_ratio+0.6541076*pmax(flc_ratio-0.2041965,0)^3-1.5410426*pmax(flc_ratio-1.0186109,0)^3+0.89034136*pmax(flc_ratio-1.6415064,0)^3-0.0034064036*pmax(flc_ratio-7.440332,0)^3+0.80496408*(isotype=="IgA")-0.30468884*(isotype=="IgG")+0.0036949171*(isotype=="Light-chain") }
## <environment: 0x7fe9e6d2d028>df_nonig$nonig_pred <- with(
df_nonig,
expit(
-0.9605 -0.7066287 * m_protein + 0.52974103 * pmax(m_protein - 0.481261, 0) ^ 3 -
0.66996903 * pmax(m_protein - 0.881164, 0) ^ 3 + 0.14177108 * pmax(m_protein - 2.48194, 0) ^ 3
- 0.0015430737 * pmax(m_protein - 10.665996, 0) ^ 3 - 1.8673982 * flc_ratio +
0.6541076 * pmax(flc_ratio - 0.2041965, 0) ^ 3 - 1.5410426 *
pmax(flc_ratio -1.0186109, 0) ^ 3 + 0.89034136 * pmax(flc_ratio - 1.6415064, 0) ^ 3
- 0.0034064036 * pmax(flc_ratio - 7.440332, 0) ^ 3 + 0.80496408 * (isotype == "IgA")
- 0.30468884 * (isotype == "IgG") + 0.0036949171 * (isotype == "Light-chain") )
)We create the decision curve analysis models using the dca pacakge. First we examine the net-benefit of the non-Ig and non-linear 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, the non-linear model outperforms the non-Ig model over a large proportion of the figure.
nonlinear_model <- decision_curve(perc_10 ~ nonlinear_pred,
data = df_nonig,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 500)
nonig_model <- decision_curve(perc_10 ~ nonig_pred,
data = df_nonig,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 500)
nonlinear_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
union_all(
nonig_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
filter(model == "perc_10 ~ nonig_pred")
) %>%
mutate(model = factor(model, levels = c("perc_10 ~ nonlinear_pred", "perc_10 ~ nonig_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-linear model",
"Defer based on non-Ig 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)
)
For posterity, we repeat this for an opt-in strategy. It is more difficult to see visually, but the same result is obtained.
nonlinear_model <- decision_curve(perc_10 ~ nonlinear_pred,
data = df_nonig,
policy = 'opt-in',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 500)
nonig_model <- decision_curve(perc_10 ~ nonig_pred,
data = df_nonig,
policy = 'opt-in',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 500)
nonlinear_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
union_all(
nonig_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
filter(model == "perc_10 ~ nonig_pred")
) %>%
mutate(model = factor(model, levels = c("perc_10 ~ nonlinear_pred", "perc_10 ~ nonig_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-linear model",
"Defer based on non-Ig 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.2)) +
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.2)
)
6 Summary
We have shown that the model that includes immunoglobulins as non-linear variables i) is well calibrated, whilst the non-immunoglobulin model is visibly less so ii) fits the data significantly better on the likelihood ratio test, providing roughly a 20% increase in the fraction of new information on this metric compared to the non-Ig model. iii) has a lower Akaike’s information criterion than does the non-Ig model. iv) substantially increases the variance of predicted risks, providing a roughly 20% in the fraction of new information on this metric compared to the non-Ig model. v) performs significantly better on decision curve analysis, than does the non-Ig model
Given these results, we have strong evidence for the inclusion of total IgG, IgA and IgM in the prediction model and we deem measuring these additional parameters worthwhile.