n <- 1000 # define sample size
set.seed(17) # so can reproduce the results
age <- rnorm(n, 50, 10)
blood.pressure <- rnorm(n, 120, 15)
cholesterol <- rnorm(n, 200, 25)
sex <- factor(sample(c('female','male'), n,TRUE))
label(age) <- 'Age' # label is in Hmisc
label(cholesterol) <- 'Total Cholesterol'
label(blood.pressure) <- 'Systolic Blood Pressure'
label(sex) <- 'Sex'
units(cholesterol) <- 'mg/dl' # uses units.default in Hmisc
units(blood.pressure) <- 'mmHg'
# 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)
ddist <- datadist(age, blood.pressure, cholesterol, sex)
options(datadist='ddist')
fit <- lrm(y ~ blood.pressure + sex * (age + rcs(cholesterol,4)),
x=TRUE, y=TRUE)
p <- Predict(fit, age=., cholesterol=., np=100)
bplot(p) # image plot for age, cholesterol with color
# coming from yhat; use default ranges for
# both continuous predictors
bplot(p, method='persp', theta=30, ticktype='detailed')
# 3-d perspective plot
bplot(p, method='contour') # contour plot
boundaries <- perimeter(age, cholesterol, lowess=TRUE)
plot(age, cholesterol) # show bivariate data density
lines(boundaries) # and perimeter that will be used for 3-D plot
p <- Predict(fit, age=., cholesterol=.)
par(mgp=c(1.7, .35, 0))
bplot(p, perim=boundaries)
# draws image() plot
# don't show estimates where data are sparse
# doesn't make sense here since vars don't interact
iLegend(p, x=c(30,40), y=c(230, 245)) # original logit scale
iLegend(p, x=c(65,75), y=c(230, 245), fun=plogis,
at=qlogis(c(.1,.25,.5,.75,.9)),
zlab='Probability') # probability scale
options(datadist=NULL)
Run the code above in your browser using DataCamp Workspace