incorrect comorbidity calculations for pulmonary and cancer codes, maybe more

[jackwasey]

Hi, I'm the author of the R package 'icd', and I'm glad to see that several of us have worked on solving the comorbidity computation problem. Just noticed your package today. Also, glad to see you live in my home country!

cc @patrickmdnet who is the author of 'medicalrisk'

I took some time this morning to compare the comorbidity computations between our packages both in speed and content. I was distressed to see we all differed from each other, particularly in the COPD/chronic lung disease and cancer/tumor categories. I dug into your source code and noticed you grep for descendents of a non-existent top-level code (498), giving a false positive for chronic lung disease with a random test code 498.82 . It is an open question what we should all be doing if potentially valid, but utterly non-existent codes appear, particularly as different annual revisions may gain or lose codes, and we would probably want to sweep them all up when looking for comorbidities.

In 'icd' I took the view that I would count non-existent descendents of extant codes, but I would exclude codes which had no parent with any association with a comorbidity.

I didn't look into the cancer side, but there were many more discrepancies which I suspect to be of the same origin.

Can I suggest randomly generating strings for testing. You can see I do this in the icd source code. I see you generated test data by sampling only valid codes.

One way for us all to work together might be for you to continue to implement comorbidities how you wish, and consider importing the 'icd' package for validation and explanation of actual codes, which I've put a lot of time into. I'm open to considering other ways for us to collaborate.

Best wishes,
Jack

[ellessenne]

Hello Jack!
Thanks for the heads up.
I choose to grep for all the codes assuming that grepping for non-existent codes should never results in hits - that code should never be encountered in practice. I use should as we all know that mistakes do happen... but that would affect all aspects of comorbidity scores.
What are your thoughts on this?
Also, about the testing for random codes: could you elaborate a bit on the benefits of doing that when it is not possible to determine what a code would be? I am not sure I completely follow 😃
Thanks!

[salmasian]

@jackwasey Excellent observation!

We need to standardize our approaches across all these packages. To that end, we need to first analyze how ICD codes change over time (eg. was 498.82 ever a valid ICD code? in this case, no) and specifically determine if any ICD9/ICD10 codes that have been added or removed over time were part of any of the categories for Charlson or Elixhauser.

We also need to discuss a higher level question: if ICD codes change over time (which they do) but algorithms like Charlson and Elixhauser are not updated periodically, who is accountable for making those algorithms "current"? Programmers or researchers?

I generally am a fan of ellessenne's approach: the fact that a data set might contain invalid ICD codes is the data curator's fault and should not be addressed by comorbidity package. However, I think having a function in comorbidity (and icd and medicalrisk) package called validate() which just checks your data to see if it contains invalid ICD codes is a great idea. And in the documentation, the user should be encouraged to runvalidate() before running comorbidity() or any other function.

PS: I usually used my own programs to extract comorbidities. Recently, for a specific project, I decided to use R; in comparing icd and comorbidity I found the latter to be much faster (in fact icd crashed with my data set of millions of rows). So while it makes sense to say that the validate() function should only be part of icd package and other packages should use that, until icd is improved in its ability to handle large data sets I think we should not introduce a dependency from comorbidity to icd.

[jackwasey]

Yes, agreed that nobody ever talks about the changing codes year to year, although in reality, these are relatively slight, and may hardly change the comorbidities results at all.

My approach with icd is as you describe: I generously include non-existent (but potentially existent) sub-codes when matching comorbidities. However, I thought adding whole 3-digit top level codes was a step too far, since these can differ significantly in character from neighboring codes. I have the sense it is unlikely new three digit codes will spring up or disappear from year to year in ICD-10, or did so with ICD-9.

Regarding your questions: researchers (like me) are pretty much forced to do what has been done before, which is to gloss over the annual changes. Ideally, we could at least audit year-to-year changes, which I did make possible with icd.

icd does have lots of detailed, carefully crafted and thoroughly tested validation code already, written largely after I had to deal with a lot of potentially invalid codes on a research project. It's GPL-3, so anyone can take the relevant code and/or tests out and put it in their own GPL project, with appropriate credit!

I'm sorry you experienced a crash with icd. I'd really appreciate bug reports if this ever happens in a reproducible way. I have a lot of test code and it has been used successfully on huge data sets, e.g. by a Propublica journalist who used it on some huge data sets to look at national health disparities. In any case, I've developed an algorithm for significantly reducing the problem for both CD-9 and ICD-10 codes which appears to be orders of magnitude faster, e.g. doing tens of millions or rows in a few seconds on a three-year old laptop.

Looking again for common ground for us to collaborate, it would be nice to share test data or test code. This can be used to compare results, and we will only grow stronger from finding how our different implementations give different results. Speed is important for big data sets, and this is something I've been working hard on, so we can also compare speed with some common benchmark code. Happy to have this discussion.

[salmasian]

Sharing test data and test code would be a great idea.

[ellessenne]

Thanks @salmasian and @jackwasey for the interesting discussion!
I like the idea of having a validate() function that returns all potential invalid codes - I plan to include more ICD revisions in comorbidity, so this would be straightforward to implement and check.
Regarding your initial comments @jackwasey: I was thinking, if you simulate valid ICD codes (e.g. using comorbidity::sample_diag()), does the output of comorbidity and icd differ?
On a final note, I agree with both @salmasian and @jackwasey comments about "accountability": I reckon we should do with what we have, but assessing the quality of the input data is a whole different topic and I am not sure this should be on us!

[jackwasey]

Hi Folks,

Please take a look at this PDF paper which I've submitted to JSS. I roughly benchmarked the three packages and the results are in the paper. The benchmark code is in the github repo if you want to run it.

I also included a more detailed discussion about whether to assign incorrect three-digit codes to a neighboring comorbidity or not.

Still interested in sharing ideas on code validation. At some point, there are diminishing returns, and it seems most researchers are just not interested in getting microscopic detail of ICD code accuracy when calculating comorbidities.

Jack

[aalexandersson]

Both your packages are great. But it would be very helpful if you can better describe or explain the differences. My question is:

Why does the regular sum score differ between the two packages for ICD-9 "Quan" Elixhauser?

I restrict the question to "Quan" Elixhauser because I prefer Elixhauser over Charlson, and it seems to me that Alessandro's package comorbidity cannot calculate "original" Elixhauser or "AHRQ" Elixhauser. I additionally restrict the question to a regular sum score and ICD-9 to keep it simple.

To illustrate the problematic difference with a reproducible example, I use the Vermont sample data that is shipped with Jack's package icd.

> # Vermont sample data
> library(icd)
> library(haven)
> library(magrittr)
> head(vermont_dx[1:10])
  visit_id   age_group    sex death DRG   DX1   DX2   DX3   DX4   DX5
1        7       40-44   male  TRUE 640 27801 03842 51881 41519 99591
2       10 75 and over female FALSE 470 71526 25000 42830  4280  4019
3       13 75 and over female FALSE 470 71535 59651 78052 27800 V8537
4       16       55-59 female FALSE 470 71535 49390 53081 27800  V140
5       37       70-74   male FALSE 462 71536  4241  2859  2720  4414
6       41       70-74   male FALSE 462 71536 V1259 V1582  V160  V171
> v <- icd_wide_to_long(vermont_dx)[c("visit_id", "icd_code")]
> head(v)
  visit_id icd_code
1        7    27801
2        7    03842
3        7    51881
4        7    41519
5        7    99591
6        7    42842

Here, I calculate a regular sum score for icd:

> # Package icd
> icd_quan <- icd9_comorbid_quan_elix(v) %>%
+   apply(2, as.integer) # convert logical to integer
> dim(icd_quan)
[1] 1000   30
> elixsum <- rowSums(icd_quan, na.rm = TRUE)
> table(elixsum)
elixsum
  0   1   2   3   4   5   6   7   8   9 
221 189 156 134 118  77  59  26  14   6 

Here, I calculate a regular sum score for comorbidity. It has an extra first variable visit_id, so the wanted variable range for the calculation is 2-32:

> # Package comorbidity
> library(comorbidity)
> comorbidity9 <- comorbidity(x = v, id = "visit_id", code = "icd_code", 
+     score = "elixhauser_icd9") 
> dim(comorbidity9)
[1] 1000   36	
> elixsum <- rowSums(comorbidity9[2:32], na.rm = TRUE)
> table(elixsum)
elixsum
  0   1   2   3   4   5   6   7   8   9  10 
221 186 156 130 117  75  59  34  12   9   1 

Why does the sum score variable elixsum differ slightly? Since the data and the calculation of the sum score are the same in my example, the comorbidity coding must differ. But how and/or why does the coding differ?

Anders

[ellessenne]

This is interesting, thank you @aalexandersson!

I replicated your example using simulated ICD9 data (only valid codes):

library(comorbidity)
library(icd)
library(tidyverse)

Simulating real ICD-9 codes:

set.seed(1)
N <- 15000
x <- data.frame(
   id = sample(1:1000, size = N, replace = TRUE),
   code = sample_diag(N, version = "ICD9_2015"),
   stringsAsFactors = FALSE)
head(x)
#>    id  code
#> 1 266  3089
#> 2 373  V554
#> 3 573 85254
#> 4 909 45384
#> 5 202  2191
#> 6 899 20970

Package icd:

icd_quan <- icd9_comorbid_quan_elix(x) %>%
   apply(2, as.integer)
elixsum <- rowSums(icd_quan, na.rm = TRUE)
table(elixsum)
#> elixsum
#>   0   1   2   3   4   5   6
#> 244 366 235 114  34   6   1

Package comorbidity:

comorbidity9 <- comorbidity(x = x, id = "id", code = "code", score = "elixhauser_icd9")
elixsum <- rowSums(comorbidity9[2:32], na.rm = TRUE)
table(elixsum)
#> elixsum
#>   0   1   2   3   4   5   6 
#> 245 355 238 113  35  11   3

These are real ICD9 codes, so there should be no difference in principle if the only difference between icd and comorbidity was that comorbidity greps codes that do not exist - all codes simulated here are real codes! There must be something else going on here, any idea @jackwasey? I will start having a closer look as soon as I can.

P.s.: @jackwasey thanks for the draft of your JSS paper. I am a bit busy at the moment with other stuff but I will definitely read it sooner rather than later! I have some opinions on assign incorrect three-digit codes to a neighboring comorbidity or not, but I want to read your thoughts more carefully before commenting on it.

Edit: I noticed that icd9_comorbid_quan_elix() returns only one column with hypertension: shouldn't Quan's version of the Elixhauser score (Quan et al., 2005) have hypertension (uncomplicated) and hypertension (complicated)? Are we actually comparing the same score here? comorbidity computes the "Enhanced ICD-9-CM" version according to the Quan et al. (2005) paper.

[jackwasey]

One immediate issue is that comorbidity doesn't account for the fact that Hypertension and Hypertension with Complications are combined in the Elixhauser schemes, although they are reported together and separately in Quan's report. E.g. Van Walraven scores use combined hypertension. The default in icd is to be consistent with all other Elixhauser comorbidity groups, but this can be changed using the argument hierarchy = FALSE. This may help.

[aalexandersson]

@jackwasey @ellessenne
The user-written Stata command elixhauser (SSC) agrees with comorbidity. Here I show the Stata output for the simulated data.

# Stata code:
elixhauser code, idvar(id) index(e)
tabm yn*
tab elixsum
	
       ELIX |
COMORBIDITY |
        SUM |      Freq.     Percent        Cum.
------------+-----------------------------------
          0 |        245       24.50       24.50
          1 |        355       35.50       60.00
          2 |        238       23.80       83.80
          3 |        113       11.30       95.10
          4 |         35        3.50       98.60
          5 |         11        1.10       99.70
          6 |          3        0.30      100.00
------------+-----------------------------------
      Total |      1,000      100.00

icd's icd9_comorbid_quan_elix() even with argument hierarchy = FALSE does NOT produce the same "Quan" Elixhauser result:

> icd_quan <- icd9_comorbid_quan_elix(x, hierarchy = FALSE) %>%
+   apply(2, as.integer)
> elixsum <- rowSums(icd_quan, na.rm = TRUE)
> table(elixsum)
> #> elixsum
> #>   0   1   2   3   4   5   6
> #> 244 362 238 113  35   6   2

Is the remaining difference because icd assigns incorrect three-digit codes to a neighboring comorbidity or because of something else?

[aalexandersson]

Which mismatched codes does icd9_comorbid_quan_elix(x, hierarchy = FALSE) count and why?

Using the Stata command cf for comparing two datasets, I notice that here are 727 mismatches of 1,000 possible:

. cf elixsum using comorbidity9, all verbose
         elixsum:  727 mismatches
                   obs    1. 1 in master; 0 in using
                   obs    3. 2 in master; 0 in using
                   obs    4. 2 in master; 1 in using
                   obs    6. 0 in master; 2 in using

That is too many to review manually, so for now I only look at the first four mismatches (as noted above).

For id 1, icd counts tumor. Why?
For id 3, icd counts PVD and Tumor. Why?
For id 4, both count Tumor (code 1478) but only icd counts Psychoses, presumably code 299.1. Why?
For id 6, icd counts nothing. Why not? comorbidity counts pcd (presumably 417.8) and lymph (presumably 201118).

[jackwasey]

If you compare the results with the IDs in the same order, there are only three mismatches. The first one I look at: icd counts renal failure for the code "Unspecified disorder resulting from impaired renal function", which comorbidity misses. icd returns the patients in the order the first appearance of each ID is presented.

I haven't had time to look any deeper, or at why the sums don't match. I'm glad for this validation, though.

[aalexandersson]

How do I sort the icd data so that the two sets of data have the same sort order? I wrongly assumed icd would not change the sort order. I am mostly a Stata user, so it is not clear to me what Jack did:

icd returns the patients in the order the first appearance of each ID is presented.

[jackwasey]

I use something like:

x_sorted_by_id <- x[order(x$id), ]

icd doesn't change the order! It uses the order it finds the IDs in the input data. The way this sample data is generated, the IDs are random and not grouped as would be expected in most health data sets. In this case, the order of patients is the exactly the same as the input data; but when there are out-of-order IDs, there is no 'right' order, so icd orders by the order of first occurrence.

Beware row name vs id column distinction in R, too.