Non-linear vs linear measurements
Elias Eythorsson
January 2nd, 2024
1 Including total IgG, IgA and IgM, and M protein or FLC ratio as linear or non-linear predictors
This script assumes that you have already ran mgus_bm_prediction_20231125. This analysis has been completed in response to the revision request by the Annals of Internal Medicine, and will test two independent hypotheses: 1) 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 modeling these as linear variables. 2) the hypothesis that an ordinal regression containing all predictors as non-linear variables using restricted cubic splines with four knots adds significant predictive value, compared to modeling all predictors (except for FLC ratio) as linear variables. Evidence for this will be supplied by a comparison of the pre-test variance (linear) and post-test variance (non-linear), a likelihood ratio test and various measures of increased explained variance, after confirming that both models are well calibrated. # Total IgG, IgA and IgM as non-linear compared to linear variables
Our rationale for testing these two hypotheses separately is based on the prior literature and clinical experience. Total IgG, IgA and IgM have a much smaller evidence base for being predictive of smoldering multiple myeloma and multiple myeloma risk. In the prior literature, studies has included these variables in a very unusual and likely inefficient (and likely unjustifiable) statistical manor, by first dichotomizing whether total IgG, IgA and IgM were below the lower limit of normal, and then counting the number outside of this limit, other than the MGUS isotype, and evaluating this number as a risk factor. We therefore have the lowest a priori reason to believe that total IgG, IgA and IgM add significant predictive information, and accordingly, the lowest reason to believe that modeling them non-linearly should add predictive information compared to modeling them linearly. We 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, B = 200))
##
## n=1043 Mean absolute error=0.013 Mean squared error=0.00019
## 0.9 Quantile of absolute error=0.02Next we fit the linear model. The linear model is a nested model of the non-linear model, in which the coefficients of the restricted cubic spline terms for total IgG, IgA, and IgM are set to zero.
fit_pred_linear <- lrm(
formula = 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
## Again, we must visually be satisfied that the calibration is equally as good as the non-linear model. We believe this is satisfied.
plot(calibrate(fit_pred_linear, B = 200))
##
## n=1043 Mean absolute error=0.016 Mean squared error=0.00032
## 0.9 Quantile of absolute error=0.0272 Predictor effect sizes
Before we evaluate the consequences of modeling total IgG, IgA and IgM as linear or non-linear predictors with regards metrics important for prediction models (precision, explained variability, calibration, discrimination and decision curve analysis), we compare the effect sizes of the predictors for both models. It is important to reiterate that effect sizes are not important when developing a prediction model and should not weigh heavily in deciding whether or not to use the linear or non-linear model.
First we examine total IgG on the absolute predicted risk scale
Predict(fit_pred_linear,
kint = 2,
igg,
fun = plogis) %>%
select(igg, yhat, lower, upper) %>%
mutate(model = "Linear") %>%
union_all(
Predict(fit_pred,
kint = 2,
igg,
fun = plogis) %>%
select(igg, yhat, lower, upper) %>%
mutate(model = "Non-linear")
) %>%
as_data_frame() %>%
ggplot(aes(x = igg, y = yhat, color = model, fill = model)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .5) +
coord_cartesian(ylim = c(0, NA)) +
geom_vline(xintercept = 6.1, lty = 2) +
geom_vline(xintercept = 15.7, lty = 2) +
scale_y_continuous(labels = scales::percent, limits = c(0, 0.14)) +
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")
Then we examine total IgG on the odds ratio scale
Predict(fit_pred_linear,
kint = 2,
igg,
ref.zero = TRUE, fun = exp) %>%
select(igg, yhat, lower, upper) %>%
mutate(model = "Linear") %>%
union_all(
Predict(fit_pred,
kint = 2,
igg,
ref.zero = TRUE, fun = exp) %>%
select(igg, yhat, lower, upper) %>%
mutate(model = "Non-linear")
) %>%
as_data_frame() %>%
ggplot(aes(x = igg, y = yhat, color = model, fill = model)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .5) +
geom_hline(yintercept = 1, lty = 2) +
geom_vline(xintercept = 6.1, lty = 2) +
geom_vline(xintercept = 15.7, lty = 2) +
scale_y_continuous(
breaks = c(0.1, 0.25, 0.5, 0.75, 1, 1.5, 2, 4, 8, 16, 32),
limits = c(NA, NA)) +
scale_x_continuous(limits = c(5, 22)) +
coord_trans(y = "log10") +
labs(x = "IgG (g/L)", y = "OR of 10% BMPC or more", title = "Total IgG") +
theme_bw() +
theme(legend.position = "bottom")
Now we examine IgA on the absolute predicted risk scale
Predict(fit_pred_linear,
kint = 2,
iga,
fun = plogis) %>%
select(iga, yhat, lower, upper) %>%
mutate(model = "Linear") %>%
union_all(
Predict(fit_pred,
kint = 2,
iga,
fun = plogis) %>%
select(iga, yhat, lower, upper) %>%
mutate(model = "Non-linear")
) %>%
as_data_frame() %>%
ggplot(aes(x = iga, y = yhat, color = model, fill = model)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .5) +
coord_cartesian(ylim = c(0, NA)) +
geom_vline(xintercept = 0.4, lty = 2) +
geom_vline(xintercept = 2.3, lty = 2) +
scale_y_continuous(labels = scales::percent, limits = c(0, 0.15)) +
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")
Then we examine total IgA on the odds ratio scale
Predict(fit_pred_linear,
kint = 2,
iga,
ref.zero = TRUE, fun = exp) %>%
select(iga, yhat, lower, upper) %>%
mutate(model = "Linear") %>%
union_all(
Predict(fit_pred,
kint = 2,
iga,
ref.zero = TRUE, fun = exp) %>%
select(iga, yhat, lower, upper) %>%
mutate(model = "Non-linear")
) %>%
as_data_frame() %>%
ggplot(aes(x = iga, y = yhat, color = model, fill = model)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .5) +
geom_hline(yintercept = 1, lty = 2) +
geom_vline(xintercept = 0.7, lty = 2) +
geom_vline(xintercept = 3.7, lty = 2) +
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") +
labs(x = "IgA (g/L)", y = "OR of 10% BMPC or more", title = "Total IgA") +
theme_bw() +
theme(legend.position = "bottom")
Finally we examine IgM on the absolute predicted risk scale
Predict(fit_pred_linear,
kint = 2,
igm,
fun = plogis) %>%
select(igm, yhat, lower, upper) %>%
mutate(model = "Linear") %>%
union_all(
Predict(fit_pred,
kint = 2,
igm,
fun = plogis) %>%
select(igm, yhat, lower, upper) %>%
mutate(model = "Non-linear")
) %>%
as_data_frame() %>%
ggplot(aes(x = igm, y = yhat, color = model, fill = model)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .5) +
coord_cartesian(ylim = c(0, NA)) +
geom_vline(xintercept = 0.4, lty = 2) +
geom_vline(xintercept = 2.3, lty = 2) +
scale_y_continuous(labels = scales::percent, limits = c(0, NA)) +
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")
and on the odds ratio scale
Predict(fit_pred_linear,
kint = 2,
igm,
ref.zero = TRUE, fun = exp) %>%
select(igm, yhat, lower, upper) %>%
mutate(model = "Linear") %>%
union_all(
Predict(fit_pred,
kint = 2,
igm,
ref.zero = TRUE, fun = exp) %>%
select(igm, yhat, lower, upper) %>%
mutate(model = "Non-linear")
) %>%
as_data_frame() %>%
ggplot(aes(x = igm, y = yhat, color = model, fill = model)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .5) +
geom_hline(yintercept = 1, lty = 2) +
geom_vline(xintercept = 0.4, lty = 2) +
geom_vline(xintercept = 2.3, lty = 2) +
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, 8)) +
coord_trans(y = "log10") +
labs(x = "IgM (g/L)", y = "OR of 10% BMPC or more", title = "Total IgM") +
theme_bw() +
theme(legend.position = "bottom")
3 Likelihood ratio test
Because the linear model is nested within the non-linear model, we can use the likelihood ratio test to test the null hypothesis that linear and non-linear models fit the data equally well, and therefore that the non-linear terms for 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_linear, fit_pred)##
## Model 1: outcome ~ igg + iga + igm + 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
## 2.292153e+01 6.000000e+00 8.231908e-04The degree to which the non-linear model fits the data better can be quantified by calculating the adequacy of the linear 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."] # linear 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 linear model compared to the non-linear model, 0.942 out of a maximum of 1
add_base <- pre_test_lr/post_test_lr
add_base## Model L.R.
## 0.942855In other words, the fraction of new information provided by non-linear model by this metric is 0.0571, which is substantial.
frac_new <- 1 - add_base
frac_new## Model L.R.
## 0.05714496The contribution of the non-linear components of total IgG, IgA and IgM can can be shown for each immunoglobulin individually, with the overall non-linear likelihood ratio Chi-squared statistic of 22.85 on 6 degrees of freedom equaling the overall likelihood ratio test above.
anova(fit_pred, igg, iga, igm)## 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
## TOTAL NONLINEAR 22.85 6 0.0008
## TOTAL 81.32 9 <.00014 Akaike’s information criterion
The question remains whether the increased fit of the non-linear model to the data compared to the linear 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)## [1] 1969.557AIC(fit_pred_linear)## [1] 1980.479Although not directly answering the question of whether we should prefer modelling immunoglobulins linearly compared to using restricted cubic slines, it should be noted that we did allow for penalized maximum likelihood estimation in the main manuscript, inculding an added penalty on non-linear terms. The full model had the lowest Akaike’s information criterion, supporting the above calculations suggesting that the non-linear model fits the data well.
penalty <- pentrace(
fit_pred,
list(
simple = seq(0, 1, by = .1),
nonlinear = seq(0, 1, by = .1)
)
)
penalty##
## Best penalty:
##
## simple nonlinear df
## 0 0 18
##
## simple nonlinear df aic bic aic.c
## 0.0 0.0 18.00000 365.11201 276.014592 364.44404
## 0.0 0.1 14.32334 75.69983 4.801354 75.27269
## 0.1 0.1 14.25522 75.87344 5.312146 75.45025
## 0.0 0.2 13.94986 73.95926 4.909454 73.55354
## 0.1 0.2 13.88750 74.12680 5.385684 73.72461
## 0.2 0.2 13.82737 74.28682 5.843318 73.88801
## 0.0 0.3 13.73193 72.89448 4.923418 72.50100
## 0.1 0.3 13.67250 73.05411 5.377168 72.66394
## 0.2 0.3 13.61513 73.20676 5.813809 72.81977
## 0.3 0.3 13.55966 73.35300 6.234615 72.96907
## 0.0 0.4 13.58051 72.14550 4.923898 71.76042
## 0.1 0.4 13.52307 72.29755 5.360291 71.91563
## 0.2 0.4 13.46755 72.44305 5.780624 72.06418
## 0.3 0.4 13.41382 72.58253 6.186068 72.20659
## 0.4 0.4 13.36176 72.71644 6.577668 72.34333
## 0.0 0.5 13.46558 71.58355 4.930884 71.20478
## 0.1 0.5 13.40965 71.72867 5.352852 71.35296
## 0.2 0.5 13.35554 71.86761 5.759601 71.49484
## 0.3 0.5 13.30314 72.00085 6.152216 71.63091
## 0.4 0.5 13.25233 72.12881 6.531670 71.76162
## 0.5 0.5 13.20302 72.25189 6.898837 71.88734
## 0.0 0.6 13.37331 71.14496 4.948989 70.77122
## 0.1 0.6 13.31864 71.28379 5.358456 70.91302
## 0.2 0.6 13.26571 71.41676 5.753419 71.04884
## 0.3 0.6 13.21441 71.54430 6.134892 71.17914
## 0.4 0.6 13.16463 71.66683 6.503786 71.30433
## 0.5 0.6 13.11629 71.78468 6.860923 71.42477
## 0.6 0.6 13.06930 71.89820 7.207042 71.54079
## 0.0 0.7 13.29631 70.79314 4.978331 70.42357
## 0.1 0.7 13.24272 70.92624 5.376673 70.55956
## 0.2 0.7 13.19081 71.05375 5.761140 70.68986
## 0.3 0.7 13.14046 71.17610 6.132686 70.81490
## 0.4 0.7 13.09159 71.29365 6.492169 70.93505
## 0.5 0.7 13.04409 71.40674 6.840360 71.05067
## 0.6 0.7 12.99789 71.51567 7.177958 71.16204
## 0.7 0.7 12.95292 71.62071 7.505600 71.26945
## 0.0 0.8 13.23016 70.50507 5.017681 70.13906
## 0.1 0.8 13.17755 70.63292 5.405941 70.26973
## 0.2 0.8 13.12655 70.75543 5.780896 70.39498
## 0.3 0.8 13.07706 70.87301 6.143443 70.51519
## 0.4 0.8 13.02899 70.98600 6.494389 70.63072
## 0.5 0.8 12.98224 71.09471 6.834466 70.74191
## 0.6 0.8 12.93676 71.19943 7.164335 70.84902
## 0.7 0.8 12.89246 71.30041 7.484598 70.95233
## 0.8 0.8 12.84928 71.39789 7.795803 71.05207
## 0.0 0.9 13.17207 70.26534 5.065509 69.90245
## 0.1 0.9 13.12035 70.38835 5.444518 70.02822
## 0.2 0.9 13.07018 70.50626 5.810741 70.14881
## 0.3 0.9 13.02147 70.61945 6.165025 70.26457
## 0.4 0.9 12.97413 70.72823 6.508133 70.37585
## 0.5 0.9 12.92808 70.83291 6.840758 70.48295
## 0.6 0.9 12.88325 70.93375 7.163527 70.58615
## 0.7 0.9 12.83956 71.03099 7.477012 70.68567
## 0.8 0.9 12.79696 71.12485 7.781735 70.78176
## 0.9 0.9 12.75539 71.21554 8.078175 70.87460
## 0.0 1.0 13.12015 70.06318 5.120345 69.70307
## 0.1 1.0 13.06925 70.18172 5.490784 69.82431
## 0.2 1.0 13.01986 70.29537 5.848917 69.94058
## 0.3 1.0 12.97188 70.40449 6.195543 70.05223
## 0.4 1.0 12.92522 70.50938 6.531387 70.15957
## 0.5 1.0 12.87981 70.61032 6.857104 70.26290
## 0.6 1.0 12.83558 70.70756 7.173291 70.36245
## 0.7 1.0 12.79246 70.80134 7.480493 70.45848
## 0.8 1.0 12.75040 70.89186 7.779206 70.55118
## 0.9 1.0 12.70934 70.97931 8.069886 70.64075
## 1.0 1.0 12.66924 71.06386 8.352954 70.727375 Comparing variances
Next we compare the variances of the predicted risks from the linear and non-linear 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, type='fitted')
post <- post[,2]
par(mgp=c(4,1,0), mar=c(6,5,2,2))
xlab <- c(
paste0("Linear\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 linear model compared tothe non-linear model is 0.974 (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 non-linear model by this metric is only 0.0256
frac_new2 <- 1 - relativ_expl
frac_new2 ## [1] 0.02564103We can show that the predicted risks between the linear and non-linear models differ most at moderate predicted risks.
plot(pre, post-pre, xlab='Predicted risk (linear immunoglobulins)',
ylab='Difference in predicted risks (non-linear - linear)', 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_linear_nonlinear <- 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: 0x7fea1100d698>df_linear_nonlinear$nonlinear_pred <- with(
df_linear_nonlinear,
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_linear) # next we extract the linear predictor from the linear model## 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: 0x7fe9f0263738>df_linear_nonlinear$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") )
)We create the decision curve analysis models using the dca pacakge. We examine the net-benefit of the linear 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, there is no perceivable difference in net benefit between the two models in this setting.
nonlinear_model <- decision_curve(perc_10 ~ nonlinear_pred,
data = df_linear_nonlinear,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 1000)
linear_model <- decision_curve(perc_10 ~ linear_pred,
data = df_linear_nonlinear,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 1000)
nonlinear_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
union_all(
linear_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
filter(model == "perc_10 ~ linear_pred")
) %>%
mutate(model = factor(model, levels = c("perc_10 ~ nonlinear_pred", "perc_10 ~ linear_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 linear 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 of Total IgG, IgA and IgM as non-linear compared to linear variables
We have shown that the model that includes immunoglobulins as non-linear variables i) fits the data significantly better on the likelihood ratio test, providing roughly a 5% increase in the fraction of new information on this metric compared to the linear model. ii) has a lower Akaike’s information criterion than does the linear model. iii) 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 linear model.
Arguments for the removal of the non-linear components of the total immunoglobulins include the principle of parsimony and the general observation that simpler models tend to perform better on external validation (bias-variance trade off). This is also supported by the linear and non-linear models performing equally well on decision curve analysis.
Arguments against the removal of the non-linear components include the fact that they were pre-specified and only considered for removal after testing (which may introduce testimation bias), plausibility (the fact that natural relationships are almost certainly never linear), and finally the rationale provided by the statistical exploration shown in this document. Added to this is the fact that modelling total immunoglobulins as non-linear functions does not add any cost to the clinician or patient (does not increased the number of measured variables), and the added complexity is not relevant given that a clinical calculator will be provided.
After taking into account the comments by the editorial team and the statistical reviewer Dr. Rita Popat, including meeting with Dr. Wee and Dr. Popat, we are willing to remove the non-linear components of the total IgG, IgA and IgM. # All variables modeled as non-linear compared to linear variables (except FLC ratio)
M protein is a strong predictor (if not the strongest predictor) of progression from MGUS to SMM or MM. We therefore consider this separately as any misspecification of M protein will likely have drastic consequences for the predictive ability of the model. From the above chapter we have made the decision, after strong recommendations from Drs. Wee and Popat, to model total IgG, IgA and IgM as linear variables. We will therefore use this as the baseline. We will explore modeling M protein as a five- four- and three-knot restricted cubic spline and compare the AIC of these as recommended by Dr. Harrell’s Regression Modeling Strategies. We will only deviate from the four-knot restricted cubic spline if there is substantial evidence of increased predictive power, as the four-knot restricted cubic spline was pre-specified. It is self evident from the prior literature and clinical experience that the relationship between FLC ratio and the outcome will be highly non-linear. We therefore do not evaluate modeling this as a linear relationship.
fit_pred_mprotein_5 <- lrm(
formula = outcome ~ igg + iga + igm + rcs(m_protein, 5) +
rcs(flc_ratio, 4) + isotype, data = df_pred, x = TRUE, y = TRUE
)
fit_pred_mprotein_5## Logistic Regression Model
##
## lrm(formula = outcome ~ igg + iga + igm + rcs(m_protein, 5) +
## 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.46 R2 0.341 C 0.760
## max |deriv| 1e-07 d.f. 13 R2(13,1043)0.296 Dxy 0.520
## Pr(> chi2) <0.0001 R2(13,882.4)0.339 gamma 0.520
## Brier 0.086 tau-a 0.327
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=5-9% 1.4741 0.4269 3.45 0.0006
## y>=10-14% -1.1455 0.4271 -2.68 0.0073
## y>=>15% -2.5262 0.4422 -5.71 <0.0001
## igg 0.0121 0.0231 0.53 0.5995
## iga 0.2840 0.0452 6.28 <0.0001
## igm -0.1751 0.0465 -3.77 0.0002
## m_protein -0.7752 0.2701 -2.87 0.0041
## m_protein' 70.1534 35.8432 1.96 0.0503
## m_protein'' -97.2704 55.3084 -1.76 0.0786
## m_protein''' 26.5642 23.7929 1.12 0.2642
## flc_ratio -2.1670 0.3534 -6.13 <0.0001
## flc_ratio' 39.7564 6.4897 6.13 <0.0001
## flc_ratio'' -93.6762 15.3845 -6.09 <0.0001
## isotype=IgA -0.1086 0.2526 -0.43 0.6671
## isotype=IgG -0.4988 0.2009 -2.48 0.0130
## isotype=Light-chain -0.4693 0.3223 -1.46 0.1454
## The four-knot model is the fit_pred_linear model shown in the first chapter
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
## fit_pred_mprotein_3 <- lrm(
formula = outcome ~ igg + iga + igm + rcs(m_protein, 3) +
rcs(flc_ratio, 4) + isotype, data = df_pred, x = TRUE, y = TRUE
)
fit_pred_mprotein_3## Logistic Regression Model
##
## lrm(formula = outcome ~ igg + iga + igm + rcs(m_protein, 3) +
## 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 360.35 R2 0.327 C 0.752
## max |deriv| 2e-09 d.f. 11 R2(11,1043)0.285 Dxy 0.505
## Pr(> chi2) <0.0001 R2(11,882.4)0.327 gamma 0.505
## Brier 0.087 tau-a 0.317
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=5-9% 0.9978 0.4089 2.44 0.0147
## y>=10-14% -1.5869 0.4116 -3.86 0.0001
## y>=>15% -2.9723 0.4273 -6.96 <0.0001
## igg 0.0229 0.0228 1.01 0.3147
## iga 0.3060 0.0448 6.83 <0.0001
## igm -0.1649 0.0463 -3.56 0.0004
## m_protein 0.0485 0.0663 0.73 0.4648
## m_protein' 0.3274 0.1112 2.95 0.0032
## flc_ratio -2.3610 0.3510 -6.73 <0.0001
## flc_ratio' 44.0973 6.4129 6.88 <0.0001
## flc_ratio'' -104.1289 15.1938 -6.85 <0.0001
## isotype=IgA -0.1397 0.2516 -0.56 0.5787
## isotype=IgG -0.4458 0.1989 -2.24 0.0250
## isotype=Light-chain -0.1425 0.3121 -0.46 0.6478
## fit_pred_mprotein_linear <- lrm(
formula = outcome ~ igg + iga + igm + m_protein +
rcs(flc_ratio, 4) + isotype, data = df_pred, x = TRUE, y = TRUE
)
fit_pred_mprotein_linear## Logistic Regression Model
##
## lrm(formula = outcome ~ igg + iga + igm + m_protein + 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 351.52 R2 0.320 C 0.747
## max |deriv| 1e-10 d.f. 10 R2(10,1043)0.279 Dxy 0.494
## Pr(> chi2) <0.0001 R2(10,882.4)0.321 gamma 0.494
## Brier 0.088 tau-a 0.310
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=5-9% 0.5858 0.3841 1.52 0.1273
## y>=10-14% -1.9852 0.3895 -5.10 <0.0001
## y>=>15% -3.3283 0.4093 -8.13 <0.0001
## igg 0.0364 0.0223 1.63 0.1024
## iga 0.2775 0.0434 6.39 <0.0001
## igm -0.1624 0.0463 -3.51 0.0004
## m_protein 0.2307 0.0257 8.99 <0.0001
## flc_ratio -2.2961 0.3485 -6.59 <0.0001
## flc_ratio' 42.9422 6.3675 6.74 <0.0001
## flc_ratio'' -101.4123 15.0865 -6.72 <0.0001
## isotype=IgA -0.0408 0.2495 -0.16 0.8702
## isotype=IgG -0.4615 0.1993 -2.32 0.0206
## isotype=Light-chain 0.1435 0.2969 0.48 0.6289
## As in the first chapter, we argue that calibration is essential, and if one of the models is obviously poorly calibrated it should not be considered in the following steps. First we show the current four-knot model with linear IgG, IgA and IgM.
plot(calibrate(fit_pred_linear, B = 200))
##
## n=1043 Mean absolute error=0.016 Mean squared error=0.00029
## 0.9 Quantile of absolute error=0.025Next we show the five-knot model, which is qualitatively identical.
plot(calibrate(fit_pred_mprotein_5, B = 200))##
## Divergence or singularity in 1 samples
##
## n=1043 Mean absolute error=0.016 Mean squared error=0.00031
## 0.9 Quantile of absolute error=0.026Next we show the three-knot model, which is qualitatively identical.
plot(calibrate(fit_pred_mprotein_3, B = 200))
##
## n=1043 Mean absolute error=0.017 Mean squared error=0.00034
## 0.9 Quantile of absolute error=0.024And lastly, the linear M protein model, which is qualitatively identical.
plot(calibrate(fit_pred_mprotein_linear, B = 200))
##
## n=1043 Mean absolute error=0.016 Mean squared error=3e-04
## 0.9 Quantile of absolute error=0.0218 Predictor effect sizes
All prior disclaimors about effect sizes and prediction models apply. First we examine M protein on the absolute predicted risk scale and show that the effect sizes are essentially similar.
Predict(fit_pred_mprotein_5,
kint = 2,
m_protein,
fun = plogis) %>%
select(m_protein, yhat, lower, upper) %>%
mutate(model = "Five-knot") %>%
union_all(
Predict(fit_pred_linear,
kint = 2,
m_protein,
fun = plogis) %>%
select(m_protein, yhat, lower, upper) %>%
mutate(model = "Four-knot")
) %>%
union_all(
Predict(fit_pred_mprotein_3,
kint = 2,
m_protein,
fun = plogis) %>%
select(m_protein, yhat, lower, upper) %>%
mutate(model = "Three-knot")
) %>%
union_all(
Predict(fit_pred_mprotein_linear,
kint = 2,
m_protein,
fun = plogis) %>%
select(m_protein, yhat, lower, upper) %>%
mutate(model = "Linear")
) %>%
as_data_frame() %>%
ggplot(aes(x = m_protein, y = yhat, color = model, fill = model)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper, color = NULL), alpha = .2) +
coord_cartesian(ylim = c(0, NA)) +
scale_y_continuous(labels = scales::percent) +
scale_x_continuous(limits = c(5, 22)) +
labs(x = "M protein (g/L)", y = "Proportion (%) having \n≥10% BMPC", title = "M protein") +
theme_bw() +
theme(legend.position = "bottom")
Then we examine total IgG on the odds ratio scale
Predict(fit_pred_mprotein_5,
kint = 2,
m_protein,
ref.zero = TRUE, fun = exp) %>%
select(m_protein, yhat, lower, upper) %>%
mutate(model = "Five-knot") %>%
union_all(
Predict(fit_pred_linear,
kint = 2,
m_protein,
ref.zero = TRUE, fun = exp) %>%
select(m_protein, yhat, lower, upper) %>%
mutate(model = "Four-knot")
) %>%
union_all(
Predict(fit_pred_mprotein_3,
kint = 2,
m_protein,
ref.zero = TRUE, fun = exp) %>%
select(m_protein, yhat, lower, upper) %>%
mutate(model = "Three-knot")
) %>%
union_all(
Predict(fit_pred_mprotein_linear,
kint = 2,
m_protein,
ref.zero = TRUE, fun = exp) %>%
select(m_protein, yhat, lower, upper) %>%
mutate(model = "Linear")
) %>%
as_data_frame() %>%
ggplot(aes(x = m_protein, y = yhat, color = model, fill = model)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper, color = NULL), alpha = .3) +
geom_hline(yintercept = 1, lty = 2) +
scale_y_continuous(
breaks = c(0.1, 0.25, 0.5, 0.75, 1, 1.5, 2, 4, 8, 16, 32, 64, 128),
limits = c(NA, NA)) +
coord_trans(y = "log10") +
labs(x = "M protein (g/L)", y = "OR of 10% BMPC or more", title = "M protein") +
theme_bw() +
theme(legend.position = "bottom")
9 Likelihood ratio test
We will do a series of comparisons. Our main aim is to decide whether the M protein can justifiably be modeled as a linear variable compared to the current prespecified four-knot restricted cubic spline. From the following, we can say that we have strong evidence that a four-knot restriced cubic spline transformation fits the data better than including M protein as a linear variable.
lrtest(fit_pred_mprotein_linear, fit_pred_linear)##
## Model 1: outcome ~ igg + iga + igm + m_protein + rcs(flc_ratio, 4) + isotype
## Model 2: outcome ~ igg + iga + igm + rcs(m_protein, 4) + rcs(flc_ratio,
## 4) + isotype
##
## L.R. Chisq d.f. P
## 2.666601e+01 2.000000e+00 1.620133e-06We further evaluate whether there is strong evidence that a five-knot restriced cubic spline fits the data better than four. There is no evidence of this.
lrtest(fit_pred_mprotein_5, fit_pred_linear)##
## Model 1: outcome ~ igg + iga + igm + rcs(m_protein, 5) + rcs(flc_ratio,
## 4) + isotype
## Model 2: outcome ~ igg + iga + igm + rcs(m_protein, 4) + rcs(flc_ratio,
## 4) + isotype
##
## L.R. Chisq d.f. P
## 0.2733163 1.0000000 0.6011158and lastly, whether the four-knot restricted cubic spline fits the data better than a three-knot restriced cubic spline. There is strong evidence for this.
lrtest(fit_pred_mprotein_3, fit_pred_linear)##
## Model 1: outcome ~ igg + iga + igm + rcs(m_protein, 3) + rcs(flc_ratio,
## 4) + isotype
## Model 2: outcome ~ igg + iga + igm + rcs(m_protein, 4) + rcs(flc_ratio,
## 4) + isotype
##
## L.R. Chisq d.f. P
## 1.784498e+01 1.000000e+00 2.396509e-05The degree to which the four-knot restricted cubic spline model fits the data better compared to the others, can be quantifiable by calculating the adequacy models using the models’ likelihood ratio Chi-squared statistics. First we compare the linear M protein model to the four-knot restricted cubic spline model.
pre_test_lr <- fit_pred_mprotein_linear$stats["Model L.R."] # linear likelihood chi-squared statistic
post_test_lr <- fit_pred_linear$stats["Model L.R."] # 4-knot likelihood chi-squared statisticBelow is a measure of the adequacy of the linear M protein model compared to the four-knot restricted cubic spline model, 0.93 out of a maximum of 1
add_base <- pre_test_lr/post_test_lr
add_base## Model L.R.
## 0.9294905In other words, the fraction of new information provided by the four-knot model by this metric is 0.070, which is substantial.
frac_new <- 1 - add_base
frac_new## Model L.R.
## 0.07050946Now we show the four-knot compared to the three-knot model.
pre_test_lr <- fit_pred_mprotein_3$stats["Model L.R."] # 3-knot likelihood chi-squared statistic
post_test_lr <- fit_pred_linear$stats["Model L.R."] # 4-knot likelihood chi-squared statisticBelow is a measure of the adequacy of the three-knot M protein model compared to the four-knot restricted cubic spline model, 0.929 out of a maximum of 1
add_base <- pre_test_lr/post_test_lr
add_base## Model L.R.
## 0.9528148In other words, the fraction of new information provided by the 4-knot model by this metric is 0.070, which is substantial.
frac_new <- 1 - add_base
frac_new## Model L.R.
## 0.0471851610 Akaike’s information criterion
The question remains whether the increased fit of the non-linear model to the data compared to the linear 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. The AIC of the four-knot restricted cubic spline is by far the lowest, after which the five-knot restricted cubic spline is the next lowest.
AIC(fit_pred_mprotein_5)## [1] 1982.206AIC(fit_pred_linear)## [1] 1980.479AIC(fit_pred_mprotein_3)## [1] 1996.324AIC(fit_pred_mprotein_linear)## [1] 2003.14511 Comparing variances
Next we compare the variances of the predicted risks from the four-knot model and the linear M protein model. At this point, we do not feel that comparing the five-knot and three-knot is needed. In the prediction setting, increased variance in risk estimates is generally desirable. As shown below, the variance is slightly higher with the four-knot model.
pre <- predict(fit_pred_mprotein_linear, type='fitted')
pre <- pre[,2]
post <- predict(fit_pred_linear, type='fitted')
post <- post[,2]
par(mgp=c(4,1,0), mar=c(6,5,2,2))
xlab <- c(
paste0("Linear\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 linear M protein model compared to the four-knot model is 0.92 (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.9210526In other words, the fraction of new information contained in the 4-knot M protein model by this metric is 0.079
frac_new2 <- 1 - relativ_expl
frac_new2 ## [1] 0.0789473712 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 the models are similar. If not, we should heavily favor the model that leads to increased net-benefit over plausible risk thresholds. Again, we feel that we no longer need to compare the 5-knot and 3-knot models and focus on the 4-knot model compared to the linear M protein model. First we extract the linear predictor from the 4-knot model. This has already been done in the chapter above and is saved in ‘df_linear_nonlinear$nonlinear_pred’. Next we extract the linear predictor from the linear M protein model.
Function(fit_pred_mprotein_linear) # next we extract the linear predictor from the linear model## function(igg = 10.21,iga = 2.167,igm = 0.718,m_protein = 1.91,flc_ratio = 1.237794,isotype = "IgG") {+0.036413502*igg+0.27748825*iga-0.16241595*igm+0.23069192*m_protein-2.2960801* flc_ratio+0.82010851*pmax(flc_ratio-0.2041965,0)^3-1.9367662*pmax(flc_ratio-1.0186109,0)^3+1.1214264*pmax(flc_ratio-1.6415064,0)^3-0.0047686997*pmax(flc_ratio-7.440332,0)^3-0.040766989*(isotype=="IgA")-0.46147411*(isotype=="IgG")+0.14347237*(isotype=="Light-chain") }
## <environment: 0x7fe9eea3bd30>df_linear_nonlinear$linear_mprotein_pred <- with(
df_linear_nonlinear,
expit(
-1.9852 +0.036413502 * igg + 0.27748825 * iga - 0.16241595 * igm + 0.23069192 * m_protein -
2.2960801 * flc_ratio + 0.82010851 * pmax(flc_ratio - 0.2041965, 0) ^ 3 -
1.9367662 * pmax(flc_ratio - 1.0186109, 0) ^ 3 + 1.1214264 *
pmax(flc_ratio - 1.6415064, 0) ^ 3 - 0.0047686997 * pmax(flc_ratio - 7.440332, 0) ^ 3
- 0.040766989 * (isotype == "IgA") - 0.46147411 * (isotype == "IgG")
+ 0.14347237 * (isotype == "Light-chain") )
)We create the decision curve analysis models using the dca pacakge. We examine the net-benefit of the 4-knot and linear M protein 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 figure, the 4-knot model is superior to the linear M protein model over the majority of plausible low-risk thresholds.
linear_mprotein_model <- decision_curve(perc_10 ~ linear_mprotein_pred,
data = df_linear_nonlinear,
policy = 'opt-out',
fitted = TRUE,
thresholds = seq(0, 1, by = .001),
bootstraps = 1000)
linear_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
union_all(
linear_mprotein_model$derived.data %>%
select(thresholds, NB, NB_lower, NB_upper, model) %>%
filter(model == "perc_10 ~ linear_mprotein_pred")
) %>%
mutate(model = factor(model, levels = c("perc_10 ~ linear_pred", "perc_10 ~ linear_mprotein_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 4-knot model",
"Defer based on linear M protein 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)
)
13 Summary of all variables modeled as non-linear compared to linear variables (except FLC ratio)
We have shown that the model that includes M protein as a non-linear variable compared to a linear variable i) fits the data significantly better on the likelihood ratio test, providing roughly a 7% increase in the fraction of new information on this metric compared to the linear M protein model. ii) has a lower Akaike’s information criterion than does the linear M protein model iii) significantly increases the variance of predicted risk, providing a roughly 7.8% increase in the fraction of new information on this metric compared to the linear M protein model. iv) provides superior net benefit to a model that includes M protein as a linear variable over the vast majority of plausible low-risk thresholds.
It is clear that the model containing M protein as a four knot restricted cubic spline transformation is superior to one that contains M protein as a linear variable. Furthermore, we found no evidence that models with three or five knot outperformed the four knot model, which should then be chosen by default, as the four knot model was pre-specified.