[R] Unexpected behaviour in rms::lrtest

[Kevin E. Thorpe]

Hello.

One of my teaching assistants was experimenting and encountered 
unexpected behaviour with the lrtest function in the rms package. It 
appears that when you have a pair of non-nested models that employ an 
RCS, the error checking for non-nested models appears not to work.

Here is a reproducible example.

 > library(rms)
Loading required package: Hmisc
Loading required package: lattice
Loading required package: survival
Loading required package: Formula
Loading required package: ggplot2

Attaching package: ‘Hmisc’

The following objects are masked from ‘package:base’:

     format.pval, units

Loading required package: SparseM

Attaching package: ‘SparseM’

The following object is masked from ‘package:base’:

     backsolve

 >
 > ### Code below that generates the data taken from the
 > ### help page for lrm()
 >
 > 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)
 >
 > exdata <- 
data.frame(y=y,age=age,blood.pressure=blood.pressure,cholesterol=cholesterol,sex=sex)
 >
 > ###
 > ### Example
 > ###
 >
 > fit1 <- lrm(y ~ blood.pressure + sex * age + cholesterol, data = exdata)
 > fit2 <- lrm(y ~ blood.pressure + age + sex * cholesterol, data = exdata)
 > lrtest(fit1, fit2) # error as expected
Error in lrtest(fit1, fit2) : models are not nested
 > fit3 <- lrm(y ~ blood.pressure + sex * age + rcs(cholesterol, 4), 
data = exdata)
 > fit4 <- lrm(y ~ blood.pressure + age + sex * rcs(cholesterol, 4), 
data = exdata)
 > lrtest(fit3,fit4) # gives result for apparently non-nested models

Model 1: y ~ blood.pressure + sex * age + rcs(cholesterol, 4)
Model 2: y ~ blood.pressure + age + sex * rcs(cholesterol, 4)

   L.R. Chisq         d.f.            P
2.043630e+01 2.000000e+00 3.650168e-05

For reference, here is my sessionInfo() although my TA got the same 
results and I do not know what his sessionInfo() was.

 > sessionInfo()
R version 3.4.3 Patched (2017-12-12 r73903)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Slackware 14.2

Matrix products: default
BLAS: /usr/local/lib64/R/lib/libRblas.so
LAPACK: /usr/local/lib64/R/lib/libRlapack.so

locale:
  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C
  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
  [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rms_5.1-2       SparseM_1.77    Hmisc_4.1-1     ggplot2_2.2.1
[5] Formula_1.2-2   survival_2.41-3 lattice_0.20-35 knitr_1.18

loaded via a namespace (and not attached):
  [1] Rcpp_0.12.14        pillar_1.0.1        compiler_3.4.3
  [4] RColorBrewer_1.1-2  plyr_1.8.4          base64enc_0.1-3
  [7] tools_3.4.3         rpart_4.1-11        digest_0.6.13
[10] polspline_1.1.12    nlme_3.1-131        tibble_1.4.1
[13] gtable_0.2.0        htmlTable_1.11.1    checkmate_1.8.5
[16] rlang_0.1.6         Matrix_1.2-12       rstudioapi_0.7
[19] mvtnorm_1.0-6       gridExtra_2.3       stringr_1.2.0
[22] cluster_2.0.6       MatrixModels_0.4-1  htmlwidgets_0.9
[25] grid_3.4.3          nnet_7.3-12         data.table_1.10.4-3
[28] foreign_0.8-69      multcomp_1.4-8      TH.data_1.0-8
[31] latticeExtra_0.6-28 magrittr_1.5        codetools_0.2-15
[34] MASS_7.3-48         scales_0.5.0        backports_1.1.2
[37] htmltools_0.3.6     splines_3.4.3       colorspace_1.3-2
[40] quantreg_5.34       sandwich_2.4-0      stringi_1.1.6
[43] acepack_1.4.1       lazyeval_0.2.1      munsell_0.4.3
[46] zoo_1.8-1


-- 
Kevin E. Thorpe
Head of Biostatistics,  Applied Health Research Centre (AHRC)
Li Ka Shing Knowledge Institute of St. Michael's Hospital
Assistant Professor, Dalla Lana School of Public Health
University of Toronto
email: kevin.thorpe at utoronto.ca  Tel: 416.864.5776  Fax: 416.864.3016