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.
# 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, ...)a fit derived by lrm or orm. The options x=TRUE and
y=TRUE must have been specified.
see validate and predab.resample
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.
deprecated and ignored. lrm through lrm.fit computes exact rank correlation coefficients as of version 6.9-0.
range of predicted probabilities over which to compute the maximum error. Default is entire range.
other arguments to pass to lrm.fit and to predab.resample (note especially
the group, cluster, and subset parameters)
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.
prints a summary, and optionally statistics for each re-fit
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.
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.
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 n
#> Dxy 0.3549 0.371 0.3391 0.0320 0.3230 10
#> R2 0.1230 0.135 0.1131 0.0215 0.1014 10
#> Intercept 0.0000 0.000 0.0165 -0.0165 0.0165 10
#> Slope 1.0000 1.000 0.8952 0.1048 0.8952 10
#> Emax 0.0000 0.000 0.0274 0.0274 0.0274 10
#> D 0.0955 0.105 0.0874 0.0178 0.0777 10
#> U -0.0020 -0.002 -0.0003 -0.0017 -0.0003 10
#> Q 0.0975 0.107 0.0877 0.0195 0.0780 10
#> B 0.2247 0.222 0.2269 -0.0045 0.2292 10
#> g 0.7562 0.795 0.7176 0.0776 0.6786 10
#> gp 0.1742 0.181 0.1662 0.0153 0.1589 10
validate(f, B=10, group=y)
#> index.orig training test optimism index.corrected n
#> Dxy 0.3549 0.373 0.3311 0.0419 0.3130 10
#> R2 0.1230 0.138 0.1080 0.0296 0.0933 10
#> Intercept 0.0000 0.000 0.0210 -0.0210 0.0210 10
#> Slope 1.0000 1.000 0.8725 0.1275 0.8725 10
#> Emax 0.0000 0.000 0.0339 0.0339 0.0339 10
#> D 0.0955 0.108 0.0832 0.0245 0.0710 10
#> U -0.0020 -0.002 0.0007 -0.0027 0.0007 10
#> Q 0.0975 0.110 0.0825 0.0272 0.0702 10
#> B 0.2247 0.222 0.2281 -0.0062 0.2309 10
#> g 0.7562 0.808 0.6997 0.1081 0.6481 10
#> gp 0.1742 0.183 0.1623 0.0211 0.1531 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 n
#> Dxy 0.3534 0.3547 0.3390 0.0157 0.3377 10
#> R2 0.1227 0.1240 0.1111 0.0129 0.1098 10
#> Intercept 0.0000 0.0000 0.0693 -0.0693 0.0693 10
#> Slope 1.0000 1.0000 0.9405 0.0595 0.9405 10
#> Emax 0.0000 0.0000 0.0261 0.0261 0.0261 10
#> D 0.0952 0.0966 0.0857 0.0109 0.0844 10
#> U -0.0020 -0.0020 0.0013 -0.0033 0.0013 10
#> Q 0.0972 0.0986 0.0844 0.0142 0.0831 10
#> B 0.2248 0.2254 0.2275 -0.0022 0.2269 10
#> g 0.7551 0.7582 0.7126 0.0456 0.7095 10
#> gp 0.1740 0.1736 0.1644 0.0092 0.1647 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
} # }