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.tho...@utoronto.ca Tel: 416.864.5776 Fax: 416.864.3016
______________________________________________
R-help@r-project.org 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.
