When I perform a logistic regression using lrm and specify weights for the observations, I get the following warning message:
Warning message: In lrm(Tag ~ DLL, weights = W, data = tagdata, x = TRUE, y = TRUE) : currently weights are ignored in model validation and bootstrapping lrm fits
My interpretation is that everything that the rms package will tell me regarding goodness-of-fit, notably using the residuals.lrm tool, is wrong. Is this correct?
To be more specific, I have working example. All the code and output can be found in this GitHub repository. I have two CSV tables of data, toystudy.csv and realstudy.csv. There are three columns in each:
- The binomial response $y$ (0 or 1) [called Tag in code]
- The predictor $x$ [called DLL in code]
- The weight for the observation [called W in code]
The former is simulated data, where all the weights are unity and where a logistic regression $log(\pi) = \theta_0 + \theta_1 x$ should fit the data perfectly. The latter is real data from my analysis, where the validity of this simple model is in question. The real data has weighted observations. (Some of the weights are negative, but there is a well-defined reason for this). The analysis code in contained completely in regressionTest.R; the meat of the code is
library(rms) fit <- lrm(Tag ~ DLL, weights = W, data = tagdata, x=TRUE, y=TRUE) residuals(fit,"gof")
Here are the results for the two tables of data.
Case 1: Toy data
The goodness-of-fit claimed by lrm (which is something called the le Cessie-van Houwelingen-Copas-Hosmer test, I understand?) is very good:
This is confirmed by grouping the data into 20 quantiles of the predictor and overlaying the predicted success rate over the average actual success rate:
Case 2: Real data
In this case, the goodness-of-fit reported by lrm is horrendous:
However, I don't think it should be that bad. Again grouping the data into quantiles, and taking into account the weights when computing the average values in each bin:
Comparing the prediction to the observed values and their standard errors, I don't think this is that bad (the error bars here depend on how the standard error on a weighted mean is computed, so they might not be 100% right, but should at least be close). On the other hand, if I produce the same plot while ignoring the weights:
I can definitely imagine this fit being as poor as the goodness-of-fit test says.
So, is residuals.rm simply ignoring the weights when it calculates its goodness-of-fit statistic? And if so, is there any R package that will do this correctly?