Unexpected behaviour in rms::lrtest

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

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: [hidden email]  Tel: 416.864.5776  Fax: 416.864.3016

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.