12  Logistic Model Case Study: Survival of Titanic Passengers

Data source: The Titanic Passenger List edited by Michael A. Findlay, originally published in Eaton & Haas (1994) Titanic: Triumph and Tragedy, Patrick Stephens Ltd, and expanded with the help of the Internet community. The original html files were obtained from Philip Hind (1999). The dataset was compiled and interpreted by Thomas Cason. It is available in R and spreadsheet formats from hbiostat.org/data under the name titanic3.

12.1 Descriptive Statistics

Code
require(rms)
options(prType='html')   # for print, summary, anova
getHdata(titanic3)        # get dataset from web site
# List of names of variables to analyze
v <- c('pclass','survived','age','sex','sibsp','parch')
t3 <- titanic3[, v]
units(t3$age) <- 'years'
describe(t3)
t3 Descriptives
t3

6 Variables   1309 Observations

pclass
image
nmissingdistinct
130903
 Value        1st   2nd   3rd
 Frequency    323   277   709
 Proportion 0.247 0.212 0.542 

survived: Survived
nmissingdistinctInfoSumMeanGmd
1309020.7085000.3820.4725

age: Age years
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
1046263980.99929.8816.06 5142128395057
lowest : 0.1667 0.3333 0.4167 0.6667 0.75 , highest: 70.5 71 74 76 80
sex
nmissingdistinct
130902
 Value      female   male
 Frequency     466    843
 Proportion  0.356  0.644 

sibsp: Number of Siblings/Spouses Aboard
image
nmissingdistinctInfoMeanGmd
1309070.670.49890.777
 Value          0     1     2     3     4     5     8
 Frequency    891   319    42    20    22     6     9
 Proportion 0.681 0.244 0.032 0.015 0.017 0.005 0.007 
For the frequency table, variable is rounded to the nearest 0
parch: Number of Parents/Children Aboard
image
nmissingdistinctInfoMeanGmd
1309080.5490.3850.6375
 Value          0     1     2     3     4     5     6     9
 Frequency   1002   170   113     8     6     6     2     2
 Proportion 0.765 0.130 0.086 0.006 0.005 0.005 0.002 0.002 
For the frequency table, variable is rounded to the nearest 0
Code
spar(ps=6,rt=3)
dd <- datadist(t3)
# describe distributions of variables to rms
options(datadist='dd')
s <- summary(survived ~ age + sex + pclass +
             cut2(sibsp,0:3) + cut2(parch,0:3), data=t3)
plot(s, main='', subtitles=FALSE)

Figure 12.1: Univariable summaries of Titanic survival

Show 4-way relationships after collapsing levels. Suppress estimates based on \(<25\) passengers.

A
Code
require(ggplot2)
tn <- transform(t3,
  agec = ifelse(age < 21, 'child', 'adult'),
  sibsp= ifelse(sibsp == 0, 'no sib/sp', 'sib/sp'),
  parch= ifelse(parch == 0, 'no par/child', 'par/child'))
g <- function(y) if(length(y) < 25) NA else mean(y)
s <- with(tn, summarize(survived,
           llist(agec, sex, pclass, sibsp, parch), g))
# llist, summarize in Hmisc package
ggplot(subset(s, agec != 'NA'),
  aes(x=survived, y=pclass, shape=sex)) +
  geom_point() + facet_grid(agec ~ sibsp * parch) +
  xlab('Proportion Surviving') + ylab('Passenger Class') +
  scale_x_continuous(breaks=c(0, .5, 1))

Figure 12.2: Multi-way summary of Titanic survival

12.3 Binary Logistic Model with Casewise Deletion of Missing Values

  • First fit a model that is saturated with respect to age, sex, pclass
  • Insufficient variation in sibsp, parch to fit complex interactions or nonlinearities.
  • With age appearing in so many terms, giving too many parameters to age creates instabilities and makes many bootstrap repetitions fail to converge or to yield singular covariance matrices
  • Use AIC to determine the global number of knots for age that is “best for the money” in terms of being the most likely to cross-validate well
Code
for(k in 3 : 5) {
  f <- lrm(survived ~ sex*pclass*rcs(age, k) +
           rcs(age, k)*(sibsp + parch), data=t3)
  cat('k=', k, '  AIC=', AIC(f), '\n')
}
k= 3   AIC= 922.9147 
k= 4   AIC= 916.6481 
k= 5   AIC= 921.2103 
  • 4 knots has best (lowest) AIC and we’ll use that going forward
  • Refit that model with x=TRUE, y=TRUE so can do likelihood ratio (LR) tests
  • But start with Wald tests
Code
f1 <- lrm(survived ~ sex*pclass*rcs(age,4) +
          rcs(age,4)*(sibsp + parch), data=t3, x=TRUE, y=TRUE)
print(f1, r2=1:4)   # print all 4 R^2 measures that use only the global LR chi-square

Logistic Regression Model

lrm(formula = survived ~ sex * pclass * rcs(age, 4) + rcs(age, 
    4) * (sibsp + parch), data = t3, x = TRUE, y = TRUE)
Frequencies of Missing Values Due to Each Variable
survived      sex   pclass      age    sibsp    parch 
       0        0        0      263        0        0 
Model Likelihood
Ratio Test
Discrimination
Indexes
Rank Discrim.
Indexes
Obs 1046 LR χ2 561.97 R21046 0.416 C 0.876
0 619 d.f. 31 R231,1046 0.398 Dxy 0.752
1 427 Pr(>χ2) <0.0001 R2758.1 0.524 γ 0.753
max |∂log L/∂β| 0.0007 R231,758.1 0.504 τa 0.363
Brier 0.129
β S.E. Wald Z Pr(>|Z|)
Intercept   -2.2942  3.4139 -0.67 0.5016
sex=male   6.3349  4.2247 1.50 0.1337
pclass=2nd   14.3540  8.4673 1.70 0.0900
pclass=3rd   3.5271  3.2329 1.09 0.2753
age   0.3671  0.2187 1.68 0.0932
age'   -0.8270  0.5684 -1.45 0.1457
age''   2.9159  2.3083 1.26 0.2065
sibsp   -0.8241  0.3173 -2.60 0.0094
parch   0.2397  0.7406 0.32 0.7462
sex=male × pclass=2nd  -13.7215  9.0533 -1.52 0.1296
sex=male × pclass=3rd   -6.3991  4.3000 -1.49 0.1367
sex=male × age   -0.5937  0.2582 -2.30 0.0215
sex=male × age'   1.2395  0.6406 1.93 0.0530
sex=male × age''   -4.3891  2.5546 -1.72 0.0858
pclass=2nd × age   -0.9460  0.4793 -1.97 0.0484
pclass=3rd × age   -0.4106  0.2097 -1.96 0.0502
pclass=2nd × age'   2.2111  1.0827 2.04 0.0411
pclass=3rd × age'   0.7450  0.5632 1.32 0.1859
pclass=2nd × age''   -8.5916  4.1621 -2.06 0.0390
pclass=3rd × age''   -2.0708  2.3726 -0.87 0.3828
age × sibsp   0.0035  0.0277 0.13 0.9005
age' × sibsp   0.1309  0.1076 1.22 0.2237
age'' × sibsp   -0.7549  0.5438 -1.39 0.1651
age × parch   0.0145  0.0468 0.31 0.7558
age' × parch   -0.1092  0.1262 -0.87 0.3869
age'' × parch   0.5123  0.5365 0.95 0.3396
sex=male × pclass=2nd × age   0.7993  0.5140 1.56 0.1199
sex=male × pclass=3rd × age   0.4755  0.2641 1.80 0.0718
sex=male × pclass=2nd × age'   -1.9165  1.1705 -1.64 0.1016
sex=male × pclass=3rd × age'   -0.7422  0.6754 -1.10 0.2719
sex=male × pclass=2nd × age''   7.6430  4.5357 1.69 0.0920
sex=male × pclass=3rd × age''   1.1688  2.8864 0.40 0.6855
Code
anova(f1)
Wald Statistics for survived
χ2 d.f. P
sex (Factor+Higher Order Factors) 187.59 12 <0.0001
All Interactions 60.55 11 <0.0001
pclass (Factor+Higher Order Factors) 100.33 16 <0.0001
All Interactions 47.44 14 <0.0001
age (Factor+Higher Order Factors) 61.35 24 <0.0001
All Interactions 37.51 21 0.0147
Nonlinear (Factor+Higher Order Factors) 28.15 16 0.0303
sibsp (Factor+Higher Order Factors) 20.38 4 0.0004
All Interactions 11.84 3 0.0080
parch (Factor+Higher Order Factors) 3.79 4 0.4349
All Interactions 3.79 3 0.2848
sex × pclass (Factor+Higher Order Factors) 43.72 8 <0.0001
sex × age (Factor+Higher Order Factors) 14.39 9 0.1093
Nonlinear (Factor+Higher Order Factors) 12.54 6 0.0510
Nonlinear Interaction : f(A,B) vs. AB 4.95 2 0.0843
pclass × age (Factor+Higher Order Factors) 18.59 12 0.0989
Nonlinear (Factor+Higher Order Factors) 15.56 8 0.0492
Nonlinear Interaction : f(A,B) vs. AB 9.22 4 0.0559
age × sibsp (Factor+Higher Order Factors) 11.84 3 0.0080
Nonlinear 2.22 2 0.3302
Nonlinear Interaction : f(A,B) vs. AB 2.22 2 0.3302
age × parch (Factor+Higher Order Factors) 3.79 3 0.2848
Nonlinear 1.02 2 0.5994
Nonlinear Interaction : f(A,B) vs. AB 1.02 2 0.5994
sex × pclass × age (Factor+Higher Order Factors) 11.24 6 0.0813
Nonlinear 10.12 4 0.0385
TOTAL NONLINEAR 28.15 16 0.0303
TOTAL INTERACTION 77.40 23 <0.0001
TOTAL NONLINEAR + INTERACTION 80.04 25 <0.0001
TOTAL 243.00 31 <0.0001

Compute the slightly more time-consuming LR tests

Code
af1 <- anova(f1, test='LR')
print(af1, which='subscripts')
Likelihood Ratio Statistics for survived
χ2 d.f. P Tested
sex (Factor+Higher Order Factors) 339.48 12 <0.0001 1,9-13,26-31
All Interactions 76.17 11 <0.0001 9-13,26-31
pclass (Factor+Higher Order Factors) 154.71 16 <0.0001 2-3,9-10,14-19,26-31
All Interactions 64.95 14 <0.0001 9-10,14-19,26-31
age (Factor+Higher Order Factors) 109.11 24 <0.0001 4-6,11-31
All Interactions 53.85 21 0.0001 11-31
Nonlinear (Factor+Higher Order Factors) 37.75 16 0.0016 5-6,12-13,16-19,21-22,24-25,28-31
sibsp (Factor+Higher Order Factors) 26.75 4 <0.0001 7,20-22
All Interactions 12.10 3 0.0070 20-22
parch (Factor+Higher Order Factors) 3.96 4 0.4109 8,23-25
All Interactions 3.95 3 0.2666 23-25
sex × pclass (Factor+Higher Order Factors) 54.58 8 <0.0001 9-10,26-31
sex × age (Factor+Higher Order Factors) 19.68 9 0.0200 11-13,26-31
Nonlinear (Factor+Higher Order Factors) 16.43 6 0.0116 12-13,28-31
Nonlinear Interaction : f(A,B) vs. AB 7.76 2 0.0206 12-13
pclass × age (Factor+Higher Order Factors) 27.45 12 0.0066 14-19,26-31
Nonlinear (Factor+Higher Order Factors) 22.59 8 0.0039 16-19,28-31
Nonlinear Interaction : f(A,B) vs. AB 12.97 4 0.0114 16-19
age × sibsp (Factor+Higher Order Factors) 12.10 3 0.0070 20-22
Nonlinear 2.26 2 0.3224 21-22
Nonlinear Interaction : f(A,B) vs. AB 2.26 2 0.3224 21-22
age × parch (Factor+Higher Order Factors) 3.95 3 0.2666 23-25
Nonlinear 1.03 2 0.5990 24-25
Nonlinear Interaction : f(A,B) vs. AB 1.03 2 0.5990 24-25
sex × pclass × age (Factor+Higher Order Factors) 14.94 6 0.0207 26-31
Nonlinear 14.00 4 0.0073 28-31
TOTAL NONLINEAR 37.75 16 0.0016 5-6,12-13,16-19,21-22,24-25,28-31
TOTAL INTERACTION 107.47 23 <0.0001 9-31
TOTAL NONLINEAR + INTERACTION 117.47 25 <0.0001 5-6,9-31
TOTAL 561.97 31 <0.0001 1-31
  • In the RMS text, 5 knots were used for age and only Wald tests were performed
  • Large \(p\)-value for the 3rd order interaction was used to justify exclusion of these highest-order interactions from the model (and one other term)
  • More evidence for 3rd order interaction from the more accurate LR test
  • Keep this model

Show the many effects of predictors.

B
Code
p <- Predict(f1, age, sex, pclass, sibsp=0, parch=0, fun=plogis)
ggplot(p)

Figure 12.5: Effects of predictors on probability of survival of Titanic passengers, estimated for zero siblings/spouses and zero parents/children
Code
ggplot(Predict(f1, sibsp, age=c(10,15,20,50), conf.int=FALSE))
#

Figure 12.6: Effect of number of siblings and spouses on the log odds of surviving, for third class males

Note that children having many siblings apparently had lower survival. Married adults had slightly higher survival than unmarried ones.

C

But moderate problem with missing data must be dealt with

12.4 Examining Missing Data Patterns

Code
spar(mfrow=c(2,2), top=1, ps=11)
na.patterns <- naclus(titanic3)
require(rpart)      # Recursive partitioning package
who.na <- rpart(is.na(age) ~ sex + pclass + survived +
                sibsp + parch, data=titanic3, minbucket=15)
naplot(na.patterns, 'na per var')
plot(who.na, margin=.1); text(who.na)
plot(na.patterns)

Figure 12.7: Patterns of missing data. Upper left panel shows the fraction of observations missing on each predictor. Lower panel depicts a hierarchical cluster analysis of missingness combinations. The similarity measure shown on the \(Y\)-axis is the fraction of observations for which both variables are missing. Right panel shows the result of recursive partitioning for predicting is.na(age). The rpart function found only strong patterns according to passenger class.
Code
spar(ps=7, rt=3)
plot(summary(is.na(age) ~ sex + pclass + survived +
             sibsp + parch, data=t3))

Figure 12.8: Univariable descriptions of proportion of passengers with missing age

But models almost always provide better descriptive statistics

Code
m <- lrm(is.na(age) ~ sex * pclass + survived + sibsp + parch,
         data=t3)
m

Logistic Regression Model

lrm(formula = is.na(age) ~ sex * pclass + survived + sibsp + 
    parch, data = t3)
Model Likelihood
Ratio Test
Discrimination
Indexes
Rank Discrim.
Indexes
Obs 1309 LR χ2 114.99 R2 0.133 C 0.703
FALSE 1046 d.f. 8 R28,1309 0.078 Dxy 0.406
TRUE 263 Pr(>χ2) <0.0001 R28,630.5 0.156 γ 0.451
max |∂log L/∂β| 5×10-6 Brier 0.148 τa 0.131
β S.E. Wald Z Pr(>|Z|)
Intercept  -2.2030  0.3641 -6.05 <0.0001
sex=male   0.6440  0.3953 1.63 0.1033
pclass=2nd  -1.0079  0.6658 -1.51 0.1300
pclass=3rd   1.6124  0.3596 4.48 <0.0001
survived  -0.1806  0.1828 -0.99 0.3232
sibsp   0.0435  0.0737 0.59 0.5548
parch  -0.3526  0.1253 -2.81 0.0049
sex=male × pclass=2nd   0.1347  0.7545 0.18 0.8583
sex=male × pclass=3rd  -0.8563  0.4214 -2.03 0.0422
Code
anova(m)
Wald Statistics for is.na(age)
χ2 d.f. P
sex (Factor+Higher Order Factors) 5.61 3 0.1324
All Interactions 5.58 2 0.0614
pclass (Factor+Higher Order Factors) 68.43 4 <0.0001
All Interactions 5.58 2 0.0614
survived 0.98 1 0.3232
sibsp 0.35 1 0.5548
parch 7.92 1 0.0049
sex × pclass (Factor+Higher Order Factors) 5.58 2 0.0614
TOTAL 82.90 8 <0.0001

pclass and parch are the important predictors of missing age.

12.5 Single Conditional Mean Imputation

Single imputation is not the preferred approach here. Click below to see this section.

First try: conditional mean imputation
Default spline transformation for age caused distribution of imputed values to be much different from non-imputed ones; constrain to linear. Also force discrete numeric variables to be linear because knots are hard to determine for them.

Code
xtrans <- transcan(~ I(age) + sex + pclass + I(sibsp) + I(parch),
                   imputed=TRUE, pl=FALSE, pr=FALSE, data=t3)
summary(xtrans)
transcan(x = ~I(age) + sex + pclass + I(sibsp) + I(parch), imputed = TRUE, 
    pr = FALSE, pl = FALSE, data = t3)

Iterations: 4 

R-squared achieved in predicting each variable:

   age    sex pclass  sibsp  parch 
 0.236  0.075  0.232  0.200  0.173 

Adjusted R-squared:

   age    sex pclass  sibsp  parch 
 0.233  0.072  0.229  0.197  0.170 

Coefficients of canonical variates for predicting each (row) variable

       age   sex   pclass sibsp parch
age           1.33  5.98  -3.16 -0.85
sex     0.04       -0.67  -0.04 -0.80
pclass  0.08 -0.32         0.14  0.02
sibsp  -0.02 -0.01  0.08         0.39
parch   0.00 -0.15  0.01   0.28      

Summary of imputed values

Starting estimates for imputed values:

   age    sex pclass  sibsp  parch 
    28      2      3      0      0 
Code
# Look at mean imputed values by sex,pclass and observed means
# age.i is age, filled in with conditional mean estimates
age.i <- with(t3, impute(xtrans, age, data=t3))
i <- is.imputed(age.i)
with(t3, tapply(age.i[i], list(sex[i],pclass[i]), mean))
            1st      2nd      3rd
female 37.64677 29.78567 21.67031
male   42.21854 32.55474 26.19231
Code
with(t3, tapply(age, list(sex,pclass), mean, na.rm=TRUE))
            1st      2nd      3rd
female 37.03759 27.49919 22.18531
male   41.02925 30.81540 25.96227
Code
dd   <- datadist(dd, age.i)
f.si <- lrm(survived ~ sex * pclass * rcs(age.i, 4) +
            rcs(age.i, 4) * (sibsp + parch), data=t3, x=TRUE, y=TRUE)
print(f.si, coefs=FALSE)

Logistic Regression Model

lrm(formula = survived ~ sex * pclass * rcs(age.i, 4) + rcs(age.i, 
    4) * (sibsp + parch), data = t3, x = TRUE, y = TRUE)
Model Likelihood
Ratio Test
Discrimination
Indexes
Rank Discrim.
Indexes
Obs 1309 LR χ2 649.29 R2 0.532 C 0.864
0 809 d.f. 31 R231,1309 0.376 Dxy 0.728
1 500 Pr(>χ2) <0.0001 R231,927 0.487 γ 0.732
max |∂log L/∂β| 0.0006 Brier 0.132 τa 0.344
Code
spar(ps=12)
p1 <- Predict(f1,   age,   pclass, sex, sibsp=0, fun=plogis)
p2 <- Predict(f.si, age.i, pclass, sex, sibsp=0, fun=plogis)
p  <- rbind('Casewise Deletion'=p1, 'Single Imputation'=p2,
            rename=c(age.i='age'))   # creates .set. variable
ggplot(p, groups='sex', ylab='Probability of Surviving')
anova(f.si, test='LR')
Likelihood Ratio Statistics for survived
χ2 d.f. P
sex (Factor+Higher Order Factors) 399.94 12 <0.0001
All Interactions 74.26 11 <0.0001
pclass (Factor+Higher Order Factors) 163.16 16 <0.0001
All Interactions 61.31 14 <0.0001
age.i (Factor+Higher Order Factors) 109.88 24 <0.0001
All Interactions 55.34 21 <0.0001
Nonlinear (Factor+Higher Order Factors) 40.70 16 0.0006
sibsp (Factor+Higher Order Factors) 28.84 4 <0.0001
All Interactions 12.81 3 0.0051
parch (Factor+Higher Order Factors) 1.55 4 0.8177
All Interactions 0.26 3 0.9681
sex × pclass (Factor+Higher Order Factors) 50.28 8 <0.0001
sex × age.i (Factor+Higher Order Factors) 19.61 9 0.0205
Nonlinear (Factor+Higher Order Factors) 15.35 6 0.0177
Nonlinear Interaction : f(A,B) vs. AB 8.33 2 0.0156
pclass × age.i (Factor+Higher Order Factors) 23.86 12 0.0213
Nonlinear (Factor+Higher Order Factors) 19.67 8 0.0117
Nonlinear Interaction : f(A,B) vs. AB 11.63 4 0.0203
age.i × sibsp (Factor+Higher Order Factors) 12.81 3 0.0051
Nonlinear 1.50 2 0.4718
Nonlinear Interaction : f(A,B) vs. AB 1.50 2 0.4718
age.i × parch (Factor+Higher Order Factors) 0.26 3 0.9681
Nonlinear 0.02 2 0.9876
Nonlinear Interaction : f(A,B) vs. AB 0.02 2 0.9876
sex × pclass × age.i (Factor+Higher Order Factors) 11.88 6 0.0647
Nonlinear 10.57 4 0.0318
TOTAL NONLINEAR 40.70 16 0.0006
TOTAL INTERACTION 108.27 23 <0.0001
TOTAL NONLINEAR + INTERACTION 117.26 25 <0.0001
TOTAL 649.29 31 <0.0001
Figure 12.9: Predicted probability of survival for males from fit using casewise deletion (bottom) and single conditional mean imputation (top). sibsp is set to zero for these predicted values.

Figure 12.10: Predicted probability of survival for males from fit using casewise deletion (bottom) and single conditional mean imputation (top). sibsp is set to zero for these predicted values.
D

12.6 Multiple Imputation

The following uses aregImpute with predictive mean matching. By default, aregImpute does not transform age when it is being predicted from the other variables. Four knots are used to transform age when used to impute other variables (not needed here as no other missings were present). Since the fraction of observations with missing age is \(\frac{263}{1309} = 0.2\) we use 20 imputations.

Force sibsp and parch to be linear for imputation, because their highly discrete distributions make it difficult to choose knots for splines.
Code
set.seed(17)         # so can reproduce random aspects
mi <- aregImpute(~ age + sex + pclass +
                 I(sibsp) + I(parch) + survived,
                 data=t3, n.impute=20, nk=4, pr=FALSE)
mi

Multiple Imputation using Bootstrap and PMM

aregImpute(formula = ~age + sex + pclass + I(sibsp) + I(parch) + 
    survived, data = t3, n.impute = 20, nk = 4, pr = FALSE)

n: 1309     p: 6    Imputations: 20     nk: 4 

Number of NAs:
     age      sex   pclass    sibsp    parch survived 
     263        0        0        0        0        0 

         type d.f.
age         s    1
sex         c    1
pclass      c    2
sibsp       l    1
parch       l    1
survived    l    1

Transformation of Target Variables Forced to be Linear

R-squares for Predicting Non-Missing Values for Each Variable
Using Last Imputations of Predictors
  age 
0.294 
Code
# Print the first 10 imputations for the first 10 passengers
#  having missing age
mi$imputed$age[1:10, 1:10]
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
16    29 71.0   62   41   24   71 48.0   30   28    33
38    42 58.0   58   64   62   28 51.0   36   29    29
41    42 32.5   55   24   58   60 54.0   47   23    54
47    31 28.5   48   37   60   50 28.5   38   42    47
60    28 42.0   38   31   58   21 45.0    2   61    42
70    38 58.0   30   17   43   39 64.0   52   33    30
71    37 46.0   30   47   30   36 47.0   65   30    40
75    62 46.0   47   70   65   54 21.0   47   46    56
81    24 25.0   17   28   36   29 42.0   56   48    41
107   42 23.0   60   41   46   58 21.0   61   33    62

Show the distribution of imputed (black) and actual ages (gray).

E
Code
plot(mi)
Ecdf(t3$age, add=TRUE, col='gray', lwd=2,
     subtitles=FALSE)

Figure 12.11: Distributions of imputed and actual ages for the Titanic dataset. Imputed values are in black and actual ages in gray.
  • Fit logistic models for 20 completed datasets and print the ratio of imputation-corrected variances to average ordinary variances.
  • Use method of Chan & Meng to get LR tests
  • This method takes final \(\hat{\beta}\) from a single model fit on 20 stacked completed datasets
  • But standard errors come from the usual Rubin’s rule and the 20 fits
  • rms::processMI computes the LR statistics from special information saved by fit.mult.impute triggered by lrt=TRUE
  • The Hmisc package runifChanged function is used to save the result and not spend 1m running it again until an input changes
  • The rms LRupdate function is run to fix likelihood ratio-related statistics (LR test, its \(p\)-value, various \(R^2\) measures) using the overall Chan & Meng model LR \(\chi^2\) computed by processMI
  • Two of the \(R^2\) printed use an effective sample size of 927 for the unbalanced binary survived variable
F
Code
runmi <- function()
  fit.mult.impute(survived ~ sex * pclass * rcs(age, 4) + rcs(age, 4) * (sibsp + parch),
                  lrm, mi, data=t3, pr=FALSE, lrt=TRUE)  # lrt implies x=TRUE y=TRUE + more
seed <- 17
f.mi <- runifChanged(runmi, seed, mi, t3)

Re-run because of changes in the following objects: mi 
Code
afmi <- processMI(f.mi, 'anova')
# Print imputation penalty indexes
prmiInfo(afmi)
Imputation penalties
Test Missing
Information
Fraction
Denominator
d.f.
χ2 Discount
sex (Factor+Higher Order Factors) 0.131 13387.9 0.869
All Interactions 0.180 6455.1 0.820
pclass (Factor+Higher Order Factors) 0.106 27217.2 0.894
All Interactions 0.154 11285.5 0.846
age (Factor+Higher Order Factors) 0.179 14281.1 0.821
All Interactions 0.175 12960.7 0.825
Nonlinear (Factor+Higher Order Factors) 0.160 11937.3 0.840
sibsp (Factor+Higher Order Factors) 0.209 1744.4 0.791
All Interactions 0.215 1235.9 0.785
parch (Factor+Higher Order Factors) 0.179 2362.9 0.821
All Interactions 0.219 1183.5 0.781
sex × pclass (Factor+Higher Order Factors) 0.153 6502.3 0.847
sex × age (Factor+Higher Order Factors) 0.210 3875.9 0.790
Nonlinear (Factor+Higher Order Factors) 0.223 2293.9 0.777
Nonlinear Interaction : f(A,B) vs. AB 0.000 Inf 1.000
pclass × age (Factor+Higher Order Factors) 0.169 7940.7 0.831
Nonlinear (Factor+Higher Order Factors) 0.186 4413.0 0.814
Nonlinear Interaction : f(A,B) vs. AB 0.181 2330.0 0.819
age × sibsp (Factor+Higher Order Factors) 0.215 1235.9 0.785
Nonlinear 0.147 1765.7 0.853
Nonlinear Interaction : f(A,B) vs. AB 0.147 1765.7 0.853
age × parch (Factor+Higher Order Factors) 0.219 1183.5 0.781
Nonlinear 0.213 837.2 0.787
Nonlinear Interaction : f(A,B) vs. AB 0.213 837.2 0.787
sex × pclass × age (Factor+Higher Order Factors) 0.215 2476.2 0.785
Nonlinear 0.260 1123.0 0.740
TOTAL NONLINEAR 0.160 11937.3 0.840
TOTAL INTERACTION 0.167 15608.7 0.833
TOTAL NONLINEAR + INTERACTION 0.165 17345.0 0.835
TOTAL 0.144 28342.6 0.856
  • None of the denominator d.f. is small enough for us to worry about the \(\chi^2\) approximation
  • Take the ratio of selected LR statistics after multiple imputation to that from casewise deletion
Code
afmi
Likelihood Ratio Statistics for survived
χ2 d.f. P
sex (Factor+Higher Order Factors) 345.17 12 <0.0001
All Interactions 59.41 11 <0.0001
pclass (Factor+Higher Order Factors) 161.47 16 <0.0001
All Interactions 50.55 14 <0.0001
age (Factor+Higher Order Factors) 101.66 24 <0.0001
All Interactions 43.61 21 0.0026
Nonlinear (Factor+Higher Order Factors) 39.97 16 0.0008
sibsp (Factor+Higher Order Factors) 24.23 4 <0.0001
All Interactions 8.94 3 0.0300
parch (Factor+Higher Order Factors) 3.19 4 0.5272
All Interactions 1.72 3 0.6329
sex × pclass (Factor+Higher Order Factors) 42.26 8 <0.0001
sex × age (Factor+Higher Order Factors) 14.42 9 0.1081
Nonlinear (Factor+Higher Order Factors) 11.47 6 0.0748
Nonlinear Interaction : f(A,B) vs. AB 7.94 2 0.0189
pclass × age (Factor+Higher Order Factors) 19.68 12 0.0734
Nonlinear (Factor+Higher Order Factors) 14.76 8 0.0639
Nonlinear Interaction : f(A,B) vs. AB 8.93 4 0.0629
age × sibsp (Factor+Higher Order Factors) 8.94 3 0.0300
Nonlinear 1.26 2 0.5313
Nonlinear Interaction : f(A,B) vs. AB 1.26 2 0.5313
age × parch (Factor+Higher Order Factors) 1.72 3 0.6329
Nonlinear 1.73 2 0.4214
Nonlinear Interaction : f(A,B) vs. AB 1.73 2 0.4214
sex × pclass × age (Factor+Higher Order Factors) 9.11 6 0.1676
Nonlinear 7.66 4 0.1050
TOTAL NONLINEAR 39.97 16 0.0008
TOTAL INTERACTION 87.90 23 <0.0001
TOTAL NONLINEAR + INTERACTION 100.00 25 <0.0001
TOTAL 567.58 31 <0.0001
Code
f.mi <- LRupdate(f.mi, afmi)
print(f.mi, r2=1:4)   # print all 4 imputation-adjusted R^2

Logistic Regression Model

fit.mult.impute(formula = survived ~ sex * pclass * rcs(age, 
    4) + rcs(age, 4) * (sibsp + parch), fitter = lrm, xtrans = mi, 
    data = t3, lrt = TRUE, pr = FALSE)
Model Likelihood
Ratio Test
Discrimination
Indexes
Rank Discrim.
Indexes
Obs 1309 LR χ2 567.58 R21309 0.352 C 0.868
0 809 d.f. 31 R231,1309 0.336 Dxy 0.736
1 500 Pr(>χ2) <0.0001 R2927 0.458 γ 0.737
max |∂log L/∂β| 3×10-8 R231,927 0.439 τa 0.347
Brier 0.130
β S.E. Wald Z Pr(>|Z|)
Intercept   -0.3199  3.2655 -0.10 0.9220
sex=male   5.8145  4.1248 1.41 0.1586
pclass=2nd   11.5383  8.2719 1.39 0.1631
pclass=3rd   2.3785  3.1614 0.75 0.4518
age   0.2701  0.2149 1.26 0.2087
age'   -0.6430  0.5367 -1.20 0.2309
age''   2.0278  2.2600 0.90 0.3696
sibsp   -0.7625  0.3165 -2.41 0.0160
parch   -0.4562  0.5576 -0.82 0.4133
sex=male × pclass=2nd  -11.5679  8.8617 -1.31 0.1918
sex=male × pclass=3rd   -6.0402  4.1905 -1.44 0.1495
sex=male × age   -0.5758  0.2578 -2.23 0.0255
sex=male × age'   1.2105  0.6099 1.98 0.0472
sex=male × age''   -3.8105  2.5114 -1.52 0.1292
pclass=2nd × age   -0.8021  0.4775 -1.68 0.0930
pclass=3rd × age   -0.3556  0.2096 -1.70 0.0898
pclass=2nd × age'   1.9084  1.0268 1.86 0.0631
pclass=3rd × age'   0.6770  0.5353 1.26 0.2059
pclass=2nd × age''   -6.6070  4.0713 -1.62 0.1046
pclass=3rd × age''   -1.8293  2.3224 -0.79 0.4309
age × sibsp   0.0070  0.0275 0.26 0.7981
age' × sibsp   0.0987  0.0986 1.00 0.3169
age'' × sibsp   -0.4979  0.5199 -0.96 0.3382
age × parch   0.0362  0.0396 0.91 0.3607
age' × parch   -0.1208  0.1115 -1.08 0.2783
age'' × parch   0.4435  0.5094 0.87 0.3839
sex=male × pclass=2nd × age   0.6870  0.5140 1.34 0.1813
sex=male × pclass=3rd × age   0.4564  0.2625 1.74 0.0821
sex=male × pclass=2nd × age'   -1.6435  1.1151 -1.47 0.1405
sex=male × pclass=3rd × age'   -0.7801  0.6367 -1.23 0.2205
sex=male × pclass=2nd × age''   5.7658  4.4553 1.29 0.1956
sex=male × pclass=3rd × age''   1.7728  2.7888 0.64 0.5250
Code
round(afmi[c(1,3,5,30), 'Chi-Square'] / af1[c(1,3,5,30), 'Chi-Square'], 3)
   sex  (Factor+Higher Order Factors) pclass  (Factor+Higher Order Factors) 
                                1.017                                 1.044 
   age  (Factor+Higher Order Factors)                                 TOTAL 
                                0.932                                 1.010 

G
  • Using all available data resulted in increases in predictive information for sex, pclass and strangely a reduction for age

For each completed dataset run bootstrap validation of model performance indexes and the nonparametric calibration curve. Because the 20 analyses of completed datasets help to average out some of the noise in bootstrap estimates we can use fewer bootstrap repetitions (100) than usual (300 or so).

Code
val <- function(fit)
  list(validate  = validate (fit, B=100),
       calibrate = calibrate(fit, B=100) )

runmi <- function()
  fit.mult.impute(       # 1m
    survived ~ sex * pclass * rcs(age,4) +
    rcs(age,4) * (sibsp + parch),
    lrm, mi, data=t3, pr=FALSE,
    fun=val, fitargs=list(x=TRUE, y=TRUE))
seed <- 19
f <- runifChanged(runmi, seed, mi, t3, val)

Re-run because of changes in the following objects: mi 


Divergence or singularity in 29 samples

Divergence or singularity in 24 samples

Divergence or singularity in 18 samples

Divergence or singularity in 23 samples

Divergence or singularity in 22 samples

Divergence or singularity in 20 samples

Divergence or singularity in 27 samples

Divergence or singularity in 20 samples

Divergence or singularity in 26 samples

Divergence or singularity in 18 samples

Divergence or singularity in 26 samples

Divergence or singularity in 24 samples

Divergence or singularity in 23 samples

Divergence or singularity in 23 samples

Divergence or singularity in 25 samples

Divergence or singularity in 26 samples

Divergence or singularity in 26 samples

Divergence or singularity in 23 samples

Divergence or singularity in 29 samples

Divergence or singularity in 25 samples

Divergence or singularity in 25 samples

Divergence or singularity in 19 samples

Divergence or singularity in 18 samples

Divergence or singularity in 25 samples

Divergence or singularity in 27 samples

Divergence or singularity in 20 samples

Divergence or singularity in 31 samples

Divergence or singularity in 22 samples

Divergence or singularity in 38 samples

Divergence or singularity in 27 samples

Divergence or singularity in 23 samples

Divergence or singularity in 32 samples

Divergence or singularity in 23 samples

Divergence or singularity in 16 samples

Divergence or singularity in 23 samples

Divergence or singularity in 20 samples

Divergence or singularity in 24 samples

Divergence or singularity in 19 samples

Divergence or singularity in 21 samples

Divergence or singularity in 22 samples

  • Display the 20 bootstrap internal validations averaged over the multiple imputations.
  • Show the 20 individual calibration curves then the first 3 in more detail followed by the overall calibration estimate
Code
val <- processMI(f, 'validate')
print(val, digits=3)
Index Original
Sample
Training
Sample
Test
Sample
Optimism Corrected
Index
Successful
Resamples
Dxy 0.739 0.754 0.728 0.026 0.713 1496
R2 0.543 0.561 0.503 0.058 0.486 1496
Intercept 0 0 -0.099 0.099 -0.099 1496
Slope 1 1 0.846 0.154 0.846 1496
Emax 0 0 0.055 0.055 0.055 1496
D 0.509 0.531 0.462 0.069 0.44 1496
U -0.002 -0.002 0.014 -0.015 0.014 1496
Q 0.511 0.532 0.448 0.085 0.426 1496
B 0.129 0.126 0.133 -0.007 0.136 1496
g 2.392 3.135 2.604 0.531 1.861 1496
gp 0.352 0.357 0.334 0.023 0.329 1496
Code
spar(mfrow=c(2,2), top=1, bot=2)
cal <- processMI(f, 'calibrate', nind=3)

n=1309   Mean absolute error=0.008   Mean squared error=0.00012
0.9 Quantile of absolute error=0.018

n=1309   Mean absolute error=0.008   Mean squared error=1e-04
0.9 Quantile of absolute error=0.016

n=1309   Mean absolute error=0.009   Mean squared error=0.00017
0.9 Quantile of absolute error=0.023

n=1309   Mean absolute error=0.009   Mean squared error=0.00017
0.9 Quantile of absolute error=0.022
Code
# plot(cal) for full-size final calibration curve

Figure 12.12: Estimated calibration curves for the Titanic risk model, accounting for multiple imputation

Figure 12.13: Estimated calibration curves for the Titanic risk model, accounting for multiple imputation

Return to the stacked fit and compare it to the fit from single imputation

Code
p1 <- Predict(f.si,  age.i, pclass, sex, sibsp=0, fun=plogis)
p2 <- Predict(f.mi,  age,   pclass, sex, sibsp=0, fun=plogis)
p  <- rbind('Single Imputation'=p1, 'Multiple Imputation'=p2,
            rename=c(age.i='age'))
ggplot(p, groups='sex', ylab='Probability of Surviving')

Figure 12.14: Predicted probability of survival for males from fit using single conditional mean imputation again (top) and multiple random draw imputation (bottom). Both sets of predictions are for sibsp=0.

12.7 Summarizing the Fitted Model

Show odds ratios for changes in predictor values

H
Code
spar(bot=1, top=0.5, ps=8)
# Get predicted values for certain types of passengers
s <- summary(f.mi, age=c(1,30), sibsp=0:1)
# override default ranges for 3 variables
plot(s, log=TRUE, main='')

Figure 12.15: Odds ratios for some predictor settings
Code
phat <- predict(f.mi,
                combos <-
         expand.grid(age=c(2,21,50),sex=levels(t3$sex),
                     pclass=levels(t3$pclass),
                     sibsp=0, parch=0), type='fitted')
# Can also use Predict(f.mi, age=c(2,21,50), sex, pclass,
#                      sibsp=0, fun=plogis)$yhat
options(digits=1)
data.frame(combos, phat)
   age    sex pclass sibsp parch phat
1    2 female    1st     0     0 0.55
2   21 female    1st     0     0 0.99
3   50 female    1st     0     0 0.96
4    2   male    1st     0     0 0.99
5   21   male    1st     0     0 0.49
6   50   male    1st     0     0 0.28
7    2 female    2nd     0     0 1.00
8   21 female    2nd     0     0 0.88
9   50 female    2nd     0     0 0.80
10   2   male    2nd     0     0 0.99
11  21   male    2nd     0     0 0.11
12  50   male    2nd     0     0 0.07
13   2 female    3rd     0     0 0.87
14  21 female    3rd     0     0 0.58
15  50 female    3rd     0     0 0.45
16   2   male    3rd     0     0 0.81
17  21   male    3rd     0     0 0.15
18  50   male    3rd     0     0 0.05
Code
options(digits=5)

We can also get predicted values by creating an R function that will evaluate the model on demand, but that only works if there are no 3rd-order interactions.

I
Code
pred.logit <- Function(f.mi)
# Note: if don't define sibsp to pred.logit, defaults to 0
plogis(pred.logit(age=c(2,21,50), sex='male', pclass='3rd'))

A nomogram could be used to obtain predicted values manually, but this is not feasible when so many interaction terms are present.

J

12.8 Bayesian Analysis

  • Repeat the multiple imputation-based approach but using a Bayesian binary logistic model
  • Using default blrm function normal priors on regression coefficients with zero mean and large SD making the priors almost flat
  • blrm uses the rcmdstan and rstan packages that provides the full power of Stan to R
  • Here we use cmdstan with rcmdstan
  • rmsb has its own caching mechanism that efficiently stores the model fit object (and all its posterior draws) and reads it back from disk install of running it again, until one of the inputs change
  • See this for more about the rmsb package
  • Could use smaller prior SDs to get penalized estimates
  • Using 4 independent Markov chain Hamiltonion posterior sampling procedures each with 1000 burn-in iterations that are discarded, and 1000 “real” iterations for a total of 4000 posterior sample draws
  • Use the first 10 multiple imputations already developed above (object mi), running the Bayesian procedure separately for 10 completed datasets
  • Merely have to stack the posterior draws into one giant sample to account for imputation and get correct posterior distribution
K
Code
# Use all available CPU cores less 1.  Each chain will be run on its
# own core.
require(rmsb)
options(mc.cores=parallel::detectCores() - 1, rmsb.backend='cmdstan')
# 10 Bayesian analyses took 3m on 11 cores
set.seed(21)
bt <- stackMI(survived ~ sex * pclass * rcs(age, 4) +
          rcs(age, 4) * (sibsp + parch),
          blrm, mi, data=t3, n.impute=10, refresh=25,
          file='bt.rds')
Initial log joint probability = -1835.57 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -695.969      0.325784      0.511413           1           1      123    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     199      -695.889     0.0173782     0.0272125      0.1729     0.01729      260    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     251      -695.888      0.014204    0.00568045           1           1      346    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1065.1 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -712.959     0.0176801     0.0049831      0.0834      0.9743      151    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     107      -712.959     0.0155748    0.00165401      0.9995      0.9995      163    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1398.32 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -702.384      0.153329     0.0616481      0.5674      0.5674      116    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     161      -702.345     0.0238888    0.00780664      0.9805      0.9805      207    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1131.31 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -697.097     0.0461989     0.0261213           1           1      140    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     148      -697.081     0.0160222    0.00829009      0.7163      0.7163      217    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1298.98 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -709.813    0.00734109     0.0102637    0.009593       0.772      141    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     113      -709.813    0.00565778    0.00330685           1           1      162    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1098.27 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -707.254     0.0212415     0.0491055      0.2018      0.8095      135    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     129      -707.237     0.0239951    0.00142576           1           1      174    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1856.31 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -711.472     0.0635463      0.403605           1           1      167    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     177      -711.469     0.0133301    0.00288863           1           1      268    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1039.36 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -712.284      0.365486     0.0279831      0.2897    0.002897      168    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     170      -712.259     0.0237121    0.00336129           1           1      268    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1062.29 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -706.193     0.0594262     0.0328537           1           1      128    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     130      -706.174     0.0329229    0.00113446           1           1      184    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Initial log joint probability = -1236.08 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
      99      -702.868     0.0179831    0.00544577           1           1      137    
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
     124      -702.867      0.013636    0.00165081           1           1      168    
Optimization terminated normally:  
  Convergence detected: relative gradient magnitude is below tolerance 
Finished in  0.1 seconds.
Code
bt

Bayesian Logistic Model

Dirichlet Priors With Concentration Parameter 0.541 for Intercepts

stackMI(formula = survived ~ sex * pclass * rcs(age, 4) + rcs(age, 
    4) * (sibsp + parch), fitter = blrm, xtrans = mi, data = t3, 
    n.impute = 10, refresh = 25, file = "bt.rds")
Mixed Calibration/
Discrimination Indexes
Discrimination
Indexes
Rank Discrim.
Indexes
Obs 1309 B 0.132 [0.129, 0.134] g 2.755 [2.352, 3.262] C 0.866 [0.862, 0.871]
0 809 gp 0.36 [0.341, 0.374] Dxy 0.733 [0.724, 0.742]
1 500 EV 0.467 [0.419, 0.506]
Draws 40000 v 8 [4.52, 14.104]
Chains 4 vp 0.111 [0.101, 0.12]
Time 19.1s
Imputations 10
p 31
Mean β Median β S.E. Lower Upper Pr(β>0) Symmetry
Intercept   -2.9180   -1.9279   5.0946  -13.8307   5.4465  0.3110  0.60
sex=male   9.7219   8.9155   5.7709   -0.2478  21.4845  0.9843  1.48
pclass=2nd   20.6712   19.4793   9.6672   4.0582  40.0654  0.9992  1.45
pclass=3rd   5.3431   4.3585   5.0030   -2.8204  16.1250  0.9042  1.69
age   0.4730   0.4161   0.3352   -0.0866   1.1784  0.9656  1.65
age'   -1.0841   -0.9802   0.7893   -2.7174   0.3109  0.0516  0.67
age''   4.1287   3.7805   3.1780   -1.7305  10.5189  0.9270  1.39
sibsp   -0.9452   -0.9300   0.3188   -1.5863  -0.3411  0.0006  0.87
parch   -0.5069   -0.5802   0.6982   -1.7930   1.1419  0.1654  1.58
sex=male × pclass=2nd  -20.8879  -19.8074  10.3886  -41.6335  -2.2696  0.0048  0.74
sex=male × pclass=3rd   -9.7203   -8.9451   5.8375  -21.7590   0.3773  0.0186  0.68
sex=male × age   -0.8588   -0.8101   0.3674   -1.6112  -0.2241  0.0008  0.68
sex=male × age'   1.7918   1.7045   0.8437   0.3132   3.5473  0.9960  1.38
sex=male × age''   -6.7443   -6.4441   3.3672  -13.5387  -0.5338  0.0091  0.76
pclass=2nd × age   -1.3651   -1.3013   0.5635   -2.4805  -0.3633  0.0002  0.71
pclass=3rd × age   -0.5800   -0.5230   0.3304   -1.2831  -0.0286  0.0086  0.60
pclass=2nd × age'   2.9890   2.8748   1.2062   0.7540   5.3614  0.9995  1.31
pclass=3rd × age'   1.1631   1.0588   0.7839   -0.2334   2.7849  0.9629  1.47
pclass=2nd × age''  -11.8192  -11.4281   4.7308  -21.3021  -3.1176  0.0011  0.80
pclass=3rd × age''   -4.0531   -3.7331   3.2041  -10.6001   1.8523  0.0810  0.74
age × sibsp   0.0171   0.0167   0.0273   -0.0358   0.0713  0.7324  1.04
age' × sibsp   0.0694   0.0688   0.0972   -0.1181   0.2632  0.7620  1.02
age'' × sibsp   -0.4731   -0.4697   0.5184   -1.4783   0.5511  0.1808  0.97
age × parch   0.0414   0.0465   0.0477   -0.0656   0.1307  0.8441  0.67
age' × parch   -0.1313   -0.1396   0.1268   -0.3666   0.1398  0.1415  1.26
age'' × parch   0.5631   0.5845   0.5641   -0.6044   1.6249  0.8461  0.87
sex=male × pclass=2nd × age   1.2556   1.1966   0.6072   0.1356   2.4435  0.9952  1.31
sex=male × pclass=3rd × age   0.7208   0.6718   0.3716   0.0648   1.4778  0.9933  1.44
sex=male × pclass=2nd × age'   -2.7456   -2.6436   1.3071   -5.3674  -0.3413  0.0060  0.80
sex=male × pclass=3rd × age'   -1.3313   -1.2436   0.8607   -3.0676   0.2594  0.0381  0.75
sex=male × pclass=2nd × age''   10.8604   10.5240   5.1564   1.0237  20.9899  0.9920  1.21
sex=male × pclass=3rd × age''   3.9649   3.6945   3.5433   -2.7075  11.1409  0.8812  1.23
  • Note that fit indexes have HPD uncertainty intervals
  • Everthing above accounts for imputation
  • Look at diagnostics
L
Code
stanDx(bt)
Diagnostics for each of 10 imputations

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)


Imputation 1 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 1.015 1.006 0.895 1.029 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.004     1322     1660
2    beta[1] 1.002     1023     1664
3    beta[2] 1.004      908      964
4    beta[3] 1.001     1531     1641
5    beta[4] 1.002     1167     1676
6    beta[5] 1.004     1035     1069
7    beta[6] 1.002     1532     1978
8    beta[7] 1.001     2822     3085
9    beta[8] 1.000     2791     2717
10   beta[9] 1.004      795      885
11  beta[10] 1.001     1700     1701
12  beta[11] 1.004     1027     1520
13  beta[12] 1.002      998     1311
14  beta[13] 1.004     1628     1681
15  beta[14] 1.003      880      928
16  beta[15] 1.002     1689     1912
17  beta[16] 1.003      836      886
18  beta[17] 1.001     2052     2338
19  beta[18] 1.004     1096     1519
20  beta[19] 1.003     2006     2752
21  beta[20] 1.002     4008     2827
22  beta[21] 1.000     4098     3179
23  beta[22] 1.002     4679     2963
24  beta[23] 1.001     3423     3170
25  beta[24] 1.001     4173     2884
26  beta[25] 1.001     4581     2916
27  beta[26] 1.005      767      891
28  beta[27] 1.001     1343     2266
29  beta[28] 1.005      789      861
30  beta[29] 1.002     1805     2457
31  beta[30] 1.003      928     1064
32  beta[31] 1.001     2203     1999

Imputation 2 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 0.967 0.916 0.911 1.053 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     1430     1768
2    beta[1] 1.001     1296     1914
3    beta[2] 1.002     1079     1126
4    beta[3] 1.001     1757     1949
5    beta[4] 1.001     1299     1924
6    beta[5] 1.002     1121     1616
7    beta[6] 1.001     1423     1725
8    beta[7] 1.000     2756     2857
9    beta[8] 1.002     4060     3356
10   beta[9] 1.003      911     1081
11  beta[10] 1.001     1922     1963
12  beta[11] 1.001     1231     1634
13  beta[12] 1.002     1147     1496
14  beta[13] 1.001     1506     2013
15  beta[14] 1.002      996     1009
16  beta[15] 1.000     1950     2339
17  beta[16] 1.002      963     1197
18  beta[17] 1.001     2957     3139
19  beta[18] 1.001     1287     1704
20  beta[19] 1.002     2370     2302
21  beta[20] 1.003     4398     2704
22  beta[21] 1.000     4826     2652
23  beta[22] 1.003     4957     2903
24  beta[23] 1.003     4259     2739
25  beta[24] 1.000     4545     2872
26  beta[25] 1.000     4594     3069
27  beta[26] 1.004      909     1254
28  beta[27] 1.001     1820     2218
29  beta[28] 1.003      891     1202
30  beta[29] 1.001     2397     2499
31  beta[30] 1.002     1077     1363
32  beta[31] 1.000     2229     2243

Imputation 3 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 0.911 0.945 0.958 0.959 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.003     1100     1705
2    beta[1] 1.001      967     1311
3    beta[2] 1.004      847      886
4    beta[3] 1.000     1943     2209
5    beta[4] 1.002     1083     1153
6    beta[5] 1.003      969     1175
7    beta[6] 1.001     1277     1615
8    beta[7] 1.001     2444     2764
9    beta[8] 1.000     2319     2324
10   beta[9] 1.004      748      750
11  beta[10] 1.001     1578     1895
12  beta[11] 1.002     1013     1426
13  beta[12] 1.002      907      826
14  beta[13] 1.001     1383     1692
15  beta[14] 1.004      866     1035
16  beta[15] 1.002     1283     1325
17  beta[16] 1.004      806      983
18  beta[17] 1.001     1490     1595
19  beta[18] 1.003     1031     1263
20  beta[19] 1.001     1815     2427
21  beta[20] 1.001     2952     2740
22  beta[21] 1.000     3275     2649
23  beta[22] 1.003     4003     2894
24  beta[23] 1.001     2101     2650
25  beta[24] 1.001     2379     1706
26  beta[25] 1.001     3106     2706
27  beta[26] 1.005      685      785
28  beta[27] 1.002     1301     1453
29  beta[28] 1.004      734      793
30  beta[29] 1.001     1604     2275
31  beta[30] 1.003      841      959
32  beta[31] 1.001     2066     2079

Imputation 4 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 0.99 1.073 0.922 0.949 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.002     1199     1900
2    beta[1] 1.003     1072     1657
3    beta[2] 1.006      837     1187
4    beta[3] 1.002     1891     2141
5    beta[4] 1.003     1105     1539
6    beta[5] 1.007      953     1323
7    beta[6] 1.004     1242     1923
8    beta[7] 1.001     2830     3034
9    beta[8] 1.000     2526     2587
10   beta[9] 1.005      696      955
11  beta[10] 1.002     1968     2144
12  beta[11] 1.004      997     1364
13  beta[12] 1.005      793     1267
14  beta[13] 1.004     1406     2156
15  beta[14] 1.006      710     1143
16  beta[15] 1.003     1495     2111
17  beta[16] 1.007      757     1097
18  beta[17] 1.001     2076     2131
19  beta[18] 1.004      865     1524
20  beta[19] 1.001     2125     2441
21  beta[20] 1.001     3053     3112
22  beta[21] 1.000     3822     3160
23  beta[22] 1.001     3843     2911
24  beta[23] 1.001     2720     2893
25  beta[24] 1.002     3316     2743
26  beta[25] 1.000     3396     2895
27  beta[26] 1.007      674      902
28  beta[27] 1.002     1561     1719
29  beta[28] 1.007      666      992
30  beta[29] 1.003     1974     2175
31  beta[30] 1.005      793     1155
32  beta[31] 1.001     2443     2130

Imputation 5 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 1.089 0.972 0.909 0.932 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.002     1504     1873
2    beta[1] 1.001     1147     1572
3    beta[2] 1.007     1092     1499
4    beta[3] 1.003     2141     2458
5    beta[4] 1.002     1214     1493
6    beta[5] 1.006     1421     1833
7    beta[6] 1.006     1640     1968
8    beta[7] 1.002     2563     2733
9    beta[8] 1.001     4016     3347
10   beta[9] 1.007     1004     1114
11  beta[10] 1.003     1929     2280
12  beta[11] 1.002     1153     1090
13  beta[12] 1.007     1226     1485
14  beta[13] 1.004     1826     2389
15  beta[14] 1.008     1039     1255
16  beta[15] 1.000     1668     2211
17  beta[16] 1.006      990     1288
18  beta[17] 1.001     1915     2474
19  beta[18] 1.003     1332     1868
20  beta[19] 1.003     2232     2344
21  beta[20] 1.001     3441     2951
22  beta[21] 1.000     4203     3293
23  beta[22] 1.001     3182     2879
24  beta[23] 1.002     3437     3048
25  beta[24] 1.002     3945     2956
26  beta[25] 1.000     4525     3068
27  beta[26] 1.009      920      872
28  beta[27] 1.002     1413     2149
29  beta[28] 1.008      906      837
30  beta[29] 1.003     2047     2294
31  beta[30] 1.003     1036     1440
32  beta[31] 1.004     2037     1802

Imputation 6 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 1.005 0.944 1.046 1.025 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.004     1048     1501
2    beta[1] 1.006      951     1125
3    beta[2] 1.003      831      968
4    beta[3] 1.001     1883     2003
5    beta[4] 1.001     1058     1433
6    beta[5] 1.003      933     1214
7    beta[6] 1.004     1361     1502
8    beta[7] 1.001     2265     2754
9    beta[8] 1.001     3647     2955
10   beta[9] 1.003      798      883
11  beta[10] 1.003     1822     2280
12  beta[11] 1.003      900     1109
13  beta[12] 1.006      943     1015
14  beta[13] 1.001     1603     1894
15  beta[14] 1.004      838      907
16  beta[15] 1.002     1475     1681
17  beta[16] 1.005      773      876
18  beta[17] 1.000     2447     2098
19  beta[18] 1.003      940     1307
20  beta[19] 1.001     2089     2136
21  beta[20] 1.001     3604     2893
22  beta[21] 1.002     3737     2413
23  beta[22] 1.001     4124     3131
24  beta[23] 1.001     4331     2585
25  beta[24] 1.002     5082     3089
26  beta[25] 1.000     4677     2875
27  beta[26] 1.005      696      820
28  beta[27] 1.002     1263     1502
29  beta[28] 1.004      726      870
30  beta[29] 1.002     1746     2371
31  beta[30] 1.003      800     1009
32  beta[31] 1.001     1899     2197

Imputation 7 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 0.977 1.006 0.96 1.002 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     1323     1968
2    beta[1] 1.001     1167     1615
3    beta[2] 1.007      784      873
4    beta[3] 1.001     2052     2219
5    beta[4] 1.003     1024     1358
6    beta[5] 1.006      886     1164
7    beta[6] 1.003     1175     1439
8    beta[7] 1.003     1571     2054
9    beta[8] 1.002     4445     2714
10   beta[9] 1.006      692      783
11  beta[10] 1.000     2012     2217
12  beta[11] 1.003      936     1253
13  beta[12] 1.005      940     1278
14  beta[13] 1.003     1340     1942
15  beta[14] 1.005      700      805
16  beta[15] 1.002     1674     2356
17  beta[16] 1.006      752      868
18  beta[17] 1.002     2118     2813
19  beta[18] 1.004      948     1329
20  beta[19] 1.002     2018     2376
21  beta[20] 1.000     4006     2861
22  beta[21] 1.001     4986     3245
23  beta[22] 1.001     5129     2789
24  beta[23] 1.001     4555     2628
25  beta[24] 1.001     4994     2997
26  beta[25] 1.000     5199     2897
27  beta[26] 1.006      649      718
28  beta[27] 1.004     1328     1736
29  beta[28] 1.005      694      826
30  beta[29] 1.003     1841     2318
31  beta[30] 1.003      871      992
32  beta[31] 1.001     2094     2049

Imputation 8 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 0.922 1.024 1.066 1.062 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.002     1059     1517
2    beta[1] 1.002      938     1441
3    beta[2] 1.002      735      911
4    beta[3] 1.001     1851     2049
5    beta[4] 1.002      877     1634
6    beta[5] 1.002      849     1275
7    beta[6] 1.002     1174     2212
8    beta[7] 1.000     2965     2598
9    beta[8] 1.000     4289     2685
10   beta[9] 1.004      614      773
11  beta[10] 1.002     1728     1744
12  beta[11] 1.002      815     1195
13  beta[12] 1.004      774     1149
14  beta[13] 1.003     1245     2388
15  beta[14] 1.003      716     1194
16  beta[15] 1.000     1893     2384
17  beta[16] 1.004      620      766
18  beta[17] 1.001     3421     3048
19  beta[18] 1.002      793     1726
20  beta[19] 1.002     2049     2007
21  beta[20] 1.001     4143     2723
22  beta[21] 1.003     4477     3378
23  beta[22] 1.002     4844     2754
24  beta[23] 1.001     4116     2843
25  beta[24] 1.001     5294     2702
26  beta[25] 0.999     5061     3312
27  beta[26] 1.004      588      717
28  beta[27] 1.001      987     1342
29  beta[28] 1.003      572      593
30  beta[29] 1.001     1506     2496
31  beta[30] 1.001      681     1154
32  beta[31] 1.002     2708     2124

Imputation 9 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 1.034 0.982 0.985 0.966 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.001     1288     1808
2    beta[1] 1.001     1308     1509
3    beta[2] 1.002      988     1110
4    beta[3] 1.002     2077     2371
5    beta[4] 1.002     1162     1756
6    beta[5] 1.002     1231     1641
7    beta[6] 1.001     1745     2317
8    beta[7] 1.001     2327     2540
9    beta[8] 1.000     3415     2936
10   beta[9] 1.002      978     1091
11  beta[10] 1.001     1627     1936
12  beta[11] 1.001     1202     1657
13  beta[12] 1.002      946     1337
14  beta[13] 1.000     1801     2142
15  beta[14] 1.001     1072     1315
16  beta[15] 1.000     1561     1551
17  beta[16] 1.001      917     1142
18  beta[17] 1.000     1576     1503
19  beta[18] 1.002     1203     1675
20  beta[19] 1.001     2156     2086
21  beta[20] 1.000     4455     2757
22  beta[21] 1.000     4832     3132
23  beta[22] 1.000     4471     2712
24  beta[23] 1.001     3473     2821
25  beta[24] 1.002     4676     2666
26  beta[25] 1.000     4105     2817
27  beta[26] 1.002      832      944
28  beta[27] 1.001     1381     2351
29  beta[28] 1.001      963      929
30  beta[29] 1.000     2081     2423
31  beta[30] 1.001     1153     1240
32  beta[31] 1.001     2425     2135

Imputation 10 

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.

EBFMI: 0.981 1.003 0.999 0.946 

   Parameter  Rhat ESS bulk ESS tail
1   alpha[1] 1.004      858     1047
2    beta[1] 1.005      822     1446
3    beta[2] 1.010      653      778
4    beta[3] 1.003     1925     2051
5    beta[4] 1.007      814     1253
6    beta[5] 1.006      687     1072
7    beta[6] 1.006     1223     2188
8    beta[7] 1.001     2145     2863
9    beta[8] 1.001     3759     2907
10   beta[9] 1.009      575      715
11  beta[10] 1.003     1631     1796
12  beta[11] 1.006      611      955
13  beta[12] 1.006      711     1031
14  beta[13] 1.007      933     1348
15  beta[14] 1.009      515      766
16  beta[15] 1.001     1348     1935
17  beta[16] 1.009      540      780
18  beta[17] 1.000     2889     2673
19  beta[18] 1.007      852     1281
20  beta[19] 1.004     1685     1747
21  beta[20] 1.001     4361     2503
22  beta[21] 1.001     4937     3326
23  beta[22] 0.999     4179     3156
24  beta[23] 1.000     3177     2595
25  beta[24] 1.000     4262     2767
26  beta[25] 1.000     3767     2631
27  beta[26] 1.008      497      667
28  beta[27] 1.005     1154     1827
29  beta[28] 1.009      460      648
30  beta[29] 1.003     1384     2463
31  beta[30] 1.004      643      932
32  beta[31] 1.002     1930     1870
Code
# Look at convergence of only 2 parameters
stanDxplot(bt, c('sex=male', 'pclass=3rd', 'age'), rev=TRUE)

  • Difficult to see but there are 40 traces (10 imputations \(\times\) 4 chains)
  • Diagnostics look good; posterior samples can be trusted
  • Plot posterior densities for select parameters
  • Also shows the 10 densities before stacking
M
Code
plot(bt, c('sex=male', 'pclass=3rd', 'age'), nrow=2)

  • Plot partial effect plots with 0.95 highest posterior density intervals
N
Code
p <- Predict(bt, age, sex, pclass, sibsp=0, fun=plogis, funint=FALSE)
ggplot(p)

  • Compute approximate measure of explained outcome variation for predictors
O
Code
plot(anova(bt))

  • Contrast second class males and females, both at 5 years and 30 years of age, all other things being equal
  • Compute 0.95 HPD interval for the contrast and a joint uncertainty region
  • Compute P(both contrasts < 0), both < -2, and P(either one < 0)
P
Code
k <- contrast(bt, list(sex='male',   age=c(5, 30), pclass='2nd'),
                  list(sex='female', age=c(5, 30), pclass='2nd'),
              cnames = c('age 5 M-F', 'age 30 M-F'))
k
            age Contrast    S.E.    Lower   Upper Pr(Contrast>0)
1age 5 M-F    5  -9.1826 6.22740 -21.7300  1.0987         0.0272
2age 30 M-F  30  -4.8996 0.61722  -6.1015 -3.6944         0.0000

Intervals are 0.95 highest posterior density intervals
Contrast is the posterior mean 
Code
plot(k)

Code
plot(k, bivar=TRUE)                        # assumes an ellipse
plot(k, bivar=TRUE, bivarmethod='kernel')  # doesn't
P <- PostF(k, pr=TRUE)
Contrast names: age 5 M-F, age 30 M-F 
Code
P(`age 5 M-F` <  0 & `age 30 M-F` <  0)    # note backticks
[1] 0.97285
Code
P(`age 5 M-F` < -2 & `age 30 M-F` < -2)
[1] 0.9019
Code
P(`age 5 M-F` <  0 | `age 30 M-F` <  0)
[1] 1

  • Show posterior distribution of predicted survival probability for a 21 year old male in third class with sibsp=0
  • Predict summarizes with a posterior mean (set posterior.summary='median' to use posterior median)
  • Frequentist multiple imputation estimate was 0.1342
Code
pmean <- Predict(bt, age=21, sex='male', pclass='3rd', sibsp=0, parch=0,
                 fun=plogis, funint=FALSE)
pmean
  age  sex pclass sibsp parch    yhat    lower   upper
1  21 male    3rd     0     0 0.14656 0.097635 0.19579

Response variable (y):  

Limits are 0.95 confidence limits
Code
p <- predict(bt,
             data.frame(age=21, sex='male', pclass='3rd', sibsp=0, parch=0),
             posterior.summary='all', fun=plogis, funint=FALSE)
plot(density(p), main='',
     xlab='Pr(survival) For One Covariate Combination')
abline(v=with(pmean, c(yhat, lower, upper)), col=alpha('blue', 0.5))

  • Compute Pr(survival probability > 0.2) for this man
Code
mean(p > 0.2)
[1] 0.023925
R software used
Package Purpose Functions
Hmisc Miscellaneous functions summary,plsmo,naclus,llist,latex, summarize,Dotplot,describe
Hmisc Imputation transcan,impute,fit.mult.impute,aregImpute,stackMI
rms Modeling datadist,lrm,rcs
Accounting for imputation processMI, LRupdate
Model presentation plot,summary,nomogram,Function,anova
Estimation Predict,summary,contrast
Model validation validate,calibrate
rmsb Misc. Bayesian blrm, stanDx,stanDxplot,plot
rpart1 Recursive partitioning rpart
  • 1 Written by Atkinson and Therneau