Resampling Validation of a Logistic or Ordinal Regression Model
validate.lrm.RdThe validate function when used on an object created by
lrm or orm does resampling validation of a logistic
regression model,
with or without backward step-down variable deletion. It provides
bias-corrected Somers' \(D_{xy}\) rank correlation, R-squared index,
the intercept and slope of an overall logistic calibration equation, the
maximum absolute difference in predicted and calibrated probabilities
\(E_{max}\), the discrimination index \(D\) (model L.R. \((\chi^2
- 1)/n\)), the unreliability index \(U\) =
difference in -2 log likelihood between un-calibrated \(X\beta\) and \(X\beta\) with overall intercept and slope
calibrated to test sample / n, the overall quality index (logarithmic
probability score) \(Q = D - U\), and the Brier or quadratic
probability score, \(B\) (the last 3 are not computed for ordinal
models), the \(g\)-index, and gp, the \(g\)-index on the
probability scale. The corrected slope can be thought of as shrinkage
factor that takes into account overfitting. For orm fits, a
subset of the above indexes is provided, Spearman's \(\rho\) is
substituted for \(D_{xy}\), and a new index is reported: pdm, the mean
absolute difference between 0.5 and the predicted probability that
\(Y\geq\) the marginal median of \(Y\).
See predab.resample for information about confidence limits.
Usage
# fit <- lrm(formula=response ~ terms, x=TRUE, y=TRUE) or orm
# S3 method for class 'lrm'
validate(fit, method="boot", B=40,
bw=FALSE, rule="aic", type="residual", sls=0.05, aics=0,
force=NULL, estimates=TRUE,
pr=FALSE, kint, Dxy.method,
emax.lim=c(0,1), ...)
# S3 method for class 'orm'
validate(fit, method="boot", B=40, bw=FALSE, rule="aic",
type="residual", sls=.05, aics=0, force=NULL, estimates=TRUE,
pr=FALSE, ...)Arguments
- fit
a fit derived by
lrmororm. The optionsx=TRUEandy=TRUEmust have been specified.- method,B,bw,rule,type,sls,aics,force,estimates,pr
see
validateandpredab.resample- kint
In the case of an ordinal model, specify which intercept to validate. Default is the middle intercept. For
validate.orm, intercept-specific quantities are not validated so this does not matter.- Dxy.method
deprecated and ignored.
lrmthroughlrm.fitcomputes exact rank correlation coefficients as of version 6.9-0.- emax.lim
range of predicted probabilities over which to compute the maximum error. Default is entire range.
- ...
other arguments to pass to
lrm.fitand topredab.resample(note especially thegroup,cluster, andsubsetparameters)
Value
a matrix with rows corresponding to \(D_{xy}\),
\(R^2\), Intercept, Slope, \(E_{max}\), \(D\),
\(U\), \(Q\), \(B\), \(g\), \(gp\), and
columns for the original index, resample estimates, indexes applied to
the whole or omitted sample using the model derived from the resample,
average optimism, corrected index, and number of successful re-samples.
For validate.orm not all columns are provided, Spearman's rho
is returned instead of \(D_{xy}\), and pdm is reported.
Details
If the original fit was created using penalized maximum likelihood estimation,
the same penalty.matrix used with the original
fit are used during validation.
See https://www.fharrell.com/post/bootcal/ for simulations of the accuracy of approximate bootstrap confidence intervals for overfitting-corrected Brier scores.
References
Miller ME, Hui SL, Tierney WM (1991): Validation techniques for logistic regression models. Stat in Med 10:1213–1226.
Harrell FE, Lee KL (1985): A comparison of the discrimination of discriminant analysis and logistic regression under multivariate normality. In Biostatistics: Statistics in Biomedical, Public Health, and Environmental Sciences. The Bernard G. Greenberg Volume, ed. PK Sen. New York: North-Holland, p. 333–343.
Examples
n <- 1000 # define sample size
age <- rnorm(n, 50, 10)
blood.pressure <- rnorm(n, 120, 15)
cholesterol <- rnorm(n, 200, 25)
sex <- factor(sample(c('female','male'), n,TRUE))
# Specify population model for log odds that Y=1
L <- .4*(sex=='male') + .045*(age-50) +
(log(cholesterol - 10)-5.2)*(-2*(sex=='female') + 2*(sex=='male'))
# Simulate binary y to have Prob(y=1) = 1/[1+exp(-L)]
y <- ifelse(runif(n) < plogis(L), 1, 0)
f <- lrm(y ~ sex*rcs(cholesterol)+pol(age,2)+blood.pressure, x=TRUE, y=TRUE)
#> number of knots in rcs defaulting to 5
#Validate full model fit
validate(f, B=10) # normally B=300
#> index.orig training test optimism index.corrected Lower Upper
#> Dxy 0.3549 0.371 0.3391 0.0320 0.3230 0.2849 0.3635
#> R2 0.1230 0.135 0.1131 0.0215 0.1014 0.0783 0.1288
#> Intercept 0.0000 0.000 0.0165 -0.0165 0.0165 -0.0583 0.1555
#> Slope 1.0000 1.000 0.8952 0.1048 0.8952 0.8069 1.0359
#> Emax 0.0000 0.000 0.0297 -0.0297 0.0297 -0.0220 0.0576
#> D 0.0955 0.105 0.0874 0.0178 0.0777 0.0584 0.1003
#> U -0.0020 -0.002 -0.0003 -0.0017 -0.0003 -0.0036 0.0027
#> Q 0.0975 0.107 0.0877 0.0195 0.0780 0.0561 0.1033
#> B 0.2247 0.222 0.2269 -0.0045 0.2292 0.2236 0.2346
#> g 0.7562 0.795 0.7176 0.0776 0.6786 0.5970 0.7755
#> gp 0.1742 0.181 0.1662 0.0153 0.1589 0.1430 0.1781
#> n
#> Dxy 10
#> R2 10
#> Intercept 10
#> Slope 10
#> Emax 10
#> D 10
#> U 10
#> Q 10
#> B 10
#> g 10
#> gp 10
validate(f, B=10, group=y)
#> index.orig training test optimism index.corrected Lower Upper n
#> Dxy 0.3549 0.373 0.3311 0.0419 0.3130 0.2229 0.3751 10
#> R2 0.1230 0.138 0.1080 0.0296 0.0933 0.0451 0.1349 10
#> Intercept 0.0000 0.000 0.0210 -0.0210 0.0210 -0.0339 0.0891 10
#> Slope 1.0000 1.000 0.8725 0.1275 0.8725 0.6686 1.1346 10
#> Emax 0.0000 0.000 0.0359 -0.0359 0.0359 -0.0312 0.0973 10
#> D 0.0955 0.108 0.0832 0.0245 0.0710 0.0309 0.1054 10
#> U -0.0020 -0.002 0.0007 -0.0027 0.0007 -0.0039 0.0158 10
#> Q 0.0975 0.110 0.0825 0.0272 0.0702 0.0204 0.1101 10
#> B 0.2247 0.222 0.2281 -0.0062 0.2309 0.2216 0.2420 10
#> g 0.7562 0.808 0.6997 0.1081 0.6481 0.4723 0.7995 10
#> gp 0.1742 0.183 0.1623 0.0211 0.1531 0.1186 0.1821 10
# two-sample validation: make resamples have same numbers of
# successes and failures as original sample
#Validate stepwise model with typical (not so good) stopping rule
validate(f, B=10, bw=TRUE, rule="p", sls=.1, type="individual")
#>
#> Backwards Step-down - Original Model
#>
#> Deleted Chi-Sq d.f. P Residual d.f. P AIC
#> blood.pressure 0.24 1 0.6249 0.24 1 0.6249 -1.76
#>
#> Approximate Estimates after Deleting Factors
#>
#> Coef S.E. Wald Z P
#> Intercept 0.8702433 2.8457736 0.3058 0.75976
#> sex=male -6.8502929 4.0911144 -1.6744 0.09405
#> cholesterol -0.0179958 0.0155344 -1.1584 0.24668
#> cholesterol' 0.0534096 0.0889959 0.6001 0.54842
#> cholesterol'' -0.3488630 0.4851078 -0.7191 0.47205
#> cholesterol''' 0.5401488 0.7134443 0.7571 0.44899
#> age 0.0543458 0.0514153 1.0570 0.29051
#> age^2 -0.0001235 0.0005229 -0.2362 0.81329
#> sex=male * cholesterol 0.0405070 0.0245927 1.6471 0.09953
#> sex=male * cholesterol' -0.0211647 0.1346974 -0.1571 0.87514
#> sex=male * cholesterol'' -0.1769058 0.7182590 -0.2463 0.80545
#> sex=male * cholesterol''' 0.5391529 1.0405621 0.5181 0.60436
#>
#> Factors in Final Model
#>
#> [1] sex cholesterol age sex * cholesterol
#> index.orig training test optimism index.corrected Lower Upper n
#> Dxy 0.3534 0.3547 0.3390 0.0157 0.3377 0.2520 0.3939 10
#> R2 0.1227 0.1240 0.1111 0.0129 0.1098 0.0501 0.1476 10
#> Intercept 0.0000 0.0000 0.0693 -0.0693 0.0693 -0.1022 0.2296 10
#> Slope 1.0000 1.0000 0.9405 0.0595 0.9405 0.6450 1.1698 10
#> Emax 0.0000 0.0000 0.0360 -0.0360 0.0360 -0.0037 0.0738 10
#> D 0.0952 0.0966 0.0857 0.0109 0.0844 0.0349 0.1155 10
#> U -0.0020 -0.0020 0.0013 -0.0033 0.0013 -0.0044 0.0067 10
#> Q 0.0972 0.0986 0.0844 0.0142 0.0831 0.0297 0.1161 10
#> B 0.2248 0.2254 0.2275 -0.0022 0.2269 0.2188 0.2400 10
#> g 0.7551 0.7582 0.7126 0.0456 0.7095 0.4913 0.8503 10
#> gp 0.1740 0.1736 0.1644 0.0092 0.1647 0.1235 0.1913 10
#>
#> Factors Retained in Backwards Elimination
#>
#> sex cholesterol age blood.pressure sex * cholesterol
#> * * * *
#> * * * *
#> * * *
#> * * * *
#> * * * *
#> * * * *
#> * * *
#> * * * *
#> * * *
#> * * *
#>
#> Frequencies of Numbers of Factors Retained
#>
#> 3 4
#> 4 6
if (FALSE) { # \dontrun{
#Fit a continuation ratio model and validate it for the predicted
#probability that y=0
u <- cr.setup(y)
Y <- u$y
cohort <- u$cohort
attach(mydataframe[u$subs,])
f <- lrm(Y ~ cohort+rcs(age,4)*sex, penalty=list(interaction=2))
validate(f, cluster=u$subs, subset=cohort=='all')
#see predab.resample for cluster and subset
} # }