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.

[aalexandersson]

Thank you, Jack. The R code snippet sorts the input data with N=15,000. But it seems to me that it is the output data with N=1,000 that needs to be sorted by id in ascending order to get comparable results. It would be very helpful to have a Minimal, Complete, and Verifiable Example (MVCE) of icd that produces the three mismatches that you claim are the only mismatches. If that is too complicated, then perhaps simply mention the three problematic IDs that you found? Then at least we can look at the actual ICD-9 codes for those three IDs and compare directly.

One part I really like with comorbidity for reproducibility is that it keeps the id variable. If icd9_comorbid_quan_elix() similarly could keep the id variable somehow, then the direct programmatic comparison I did with cf in Stata would work and be trivial.

[jackwasey]

Will do. FYI, icd returns a matrix of logicals by default because it is much more compact and faster with big data sets. If you add the argument return_df = TRUE you get a data.frame with an id column. Otherwise, with the matrix result, rownames(icd9_comorbid_quan_elix(x)) gives you the IDs, and the colnames are the names of the comorbidities.

[ellessenne]

Hi guys,
I modified my example as before to obtain IDs of individuals whose score differ:

library(comorbidity)                                                                                         
library(icd)                                                                                                 
#> Welcome to the 'icd' package for finding comorbidities and interpretation of ICD-9 and ICD-10 codes.
#> ?icd to get started, then see the vignettes and help for details and examples.
#> 
#> Suggestions and contributions are welcome at https://github.com/jackwasey/icd . Please cite this package if you find it useful in your published work citation(package = "icd")
library(tidyverse)                                                                                           
#> ── Attaching packages ──────────────────────────────────────────────── tidyverse 1.2.1 ──
#> ✔ ggplot2 2.2.1     ✔ purrr   0.2.4
#> ✔ tibble  1.4.2     ✔ dplyr   0.7.4
#> ✔ tidyr   0.8.0     ✔ stringr 1.3.1
#> ✔ readr   1.1.1     ✔ forcats 0.3.0
#> ── Conflicts ─────────────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::explain() masks icd::explain()
#> ✖ dplyr::filter()  masks stats::filter()
#> ✖ dplyr::lag()     masks stats::lag()
                                                                                                             
# 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"), # the codes simulated via sample_diag are only valid codes     
stringsAsFactors = FALSE                                                                                     
)                                                                                                            
head(x)                                                                                                      
#>    id  code
#> 1 266  3089
#> 2 373  V554
#> 3 573 85254
#> 4 909 45384
#> 5 202  2191
#> 6 899 20970
x <- x[order(x[["id"]]), ]                                                                                   
row.names(x) <- NULL                                                                                         
head(x)                                                                                                      
#>   id  code
#> 1  1  7966
#> 2  1 01226
#> 3  1 64260
#> 4  1 23877
#> 5  1 11509
#> 6  1 74710

# Package icd                                                                                                
icd_quan <- icd9_comorbid_quan_elix(x, hierarchy = FALSE)                                                    
icd_quan <- data.frame(id.icd = row.names(icd_quan), score.icd = rowSums(icd_quan), stringsAsFactors = FALSE)

# Package comorbidity                                                                                        
comorbidity9 <- comorbidity(x = x, id = "id", code = "code", score = "elixhauser_icd9") %>%                  
mutate(id.comorbidity = as.character(id)) %>%                                                                
rename(score.comorbidity = score) %>%                                                                        
select(id.comorbidity, score.comorbidity)                                                                    

# Merge and id those that differ                                                                             
singledf <- full_join(icd_quan, comorbidity9, by = c("id.icd" = "id.comorbidity")) %>%                       
mutate(same = as.numeric(score.icd == score.comorbidity)) %>%                                                
rename(id = id.icd)                                                                                          
x[["id"]] <- as.character(x[["id"]])                                                                         
x <- left_join(x, singledf, "id")                                                                            

# Show                                                                                                       
filter(x, same == 0)                                                                                         
#>      id  code score.icd score.comorbidity same
#> 1    43  1639         2                 3    0
#> 2    43  6114         2                 3    0
#> 3    43  9658         2                 3    0
#> 4    43 39890         2                 3    0
#> 5    43 52525         2                 3    0
#> 6    43 83119         2                 3    0
#> 7    43 40211         2                 3    0
#> 8    43 66901         2                 3    0
#> 9    43 40211         2                 3    0
#> 10   43 33901         2                 3    0
#> 11   43 94502         2                 3    0
#> 12   43 37312         2                 3    0
#> 13   43 71680         2                 3    0
#> 14   43  3545         2                 3    0
#> 15   43 V1221         2                 3    0
#> [...]

I kept only the output for the individual with ID == 43. I need some spare time to check the codes one by one for individuals that differ, I hope to get this done soon... it's pretty busy over here these days.

Anyway, I think there is no real, unique solution for the "keeping the ID" problem: it is more efficient to use matrices (as evident by your benchmarks, @jackwasey), but on the other side they can only contain single-type data. I designed comorbidity to return the ID column as it becomes easy to join the result with a "main" data set!

[aalexandersson]

Thanks! So, there are 35 mismatched ids? That's worse than 3 but much better than 727 :-)

> table(singledf$same)

  0   1 
 35 965

[jackwasey]

Alright, thank you! I still find three minor renal codes are different. You found a minor error in the renal section of the transcribed Quan Elixhauser ICD-9 map: I had accidentally written 588 instead of 588.0 . The following code shows the output is now identical:

library(comorbidity)
library(icd)
requireNamespace("testthat")
library(tidyverse)
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
)
x <- x[order(x[["id"]]), ]
row.names(x) <- NULL
# Package icd
icd_quan <- icd9_comorbid_quan_elix(x, hierarchy = FALSE, return_df = TRUE)
# Package comorbidity
cmb_quan <- comorbidity(x = x, id = "id", code = "code", score = "elixhauser_icd9")[1:32]
# change T/F to 1/0. (This will be done with new argument return_binary in next icd release.)
cmb_quan[2:32] <- icd:::binary_to_logical(cmb_quan[2:32])
# make id column and column names identical
colnames(cmb_quan) <- colnames(icd_quan)
# make sure the content is identical before converting:
stopifnot(identical(as.integer(icd_quan$id), cmb_quan$id))
cmb_quan$id <- icd_quan$id
testthat::compare(icd_quan, cmb_quan)
# now identical
# http://www.icd9data.com/2015/Volume1/580-629/580-589/588/default.htm
# link shows the few renal ICD codes which icd included inadvertently under renal category in Quan Elix map

There is still a question of the scores, which I'll get to. Most people use Van Walraven system for scoring Elixhauser comorbidities, not just the sum of the number of flags. What scoring system did you implement?

[ellessenne]

comorbidity returns both a simple sum score and a weighted score as well; weights for the Charlson score are according to the original formulation by Charlson et al. in 1987, while weights for the Elixhauser score are based on work by Moore et al. and van Walraven et al.
A categorisation of the scores is returned as well, based on work by Menendez et al.

References:

• Charlson ME, Pompei P, Ales KL, et al. A new method of classifying prognostic comorbidity in longitudinal studies: development and validation. Journal of Chronic Diseases 1987; 40:373-383.
• Moore BJ, White S, Washington R, Coenen N, and Elixhauser A. Identifying increased risk of readmission and in-hospital mortality using hospital administrative data: the AHRQ Elixhauser comorbidity index. Medical Care 2017; 55(7):698-705.
• van Walraven C, Austin PC, Jennings A, Quan H and Forster AJ. A modification of the Elixhauser comorbidity measures into a point system for hospital death using administrative data. Medical Care 2009; 47(6):626-633.
• Menendez ME, Neuhaus V, van Dijk CN, Ring D. The Elixhauser comorbidity method outperforms the Charlson index in predicting inpatient death after orthopaedic surgery. Clinical Orthopaedics and Related Research 2014; 472(9):2878-2886.

I have an unrelated question though: in my experience as an applied researcher, I never used the actual weighted comorbidity score for anything; I always ended up using single comorbidities or even adding additional comorbidities based on ICD codes, as I was dealing with large health-care data and registries. What is your experience on this point? I reckon this is probably just my biased experience, so that's why I thought about asking! I am genuinely curious here 😄

[jackwasey]

Interesting to hear that. I also use the comorbidities directly, usually for matching. icd has icd::van_walraven score function, which was contributed, and based on van Walraven C, Austin PC, Jennings A, Quan H and Forster AJ. Let's compare!

[aalexandersson]

My experience is the same. However, applied researchers differ widely in skills and needs. I work for the Florida cancer data registry.

Before discussing scoring other than sum scores, it would help me if we first could resolve the mismatches so far: Do we agree that there are 35 mismatched IDs in the example data in terms of sum scores? Jack wrote:

The following code shows the output is now identical

But I got 29 "element mismatches":

testthat::compare(icd_quan, cmb_quan)
Component “HTNcx”: 8 element mismatches
Component “Pulmonary”: 1 element mismatch
Component “Renal”: 14 element mismatches
Component “Alcohol”: 6 element mismatches
Component “Depression”: 1 element mismatch> # now identical

[jackwasey]

Would you please make sure icd is up-to-date and also try devtools::install_github("jackwasey/icd") which also fixes a recent regression related to duplicate codes in comorbidities which might explain those few discrepancies.

[aalexandersson]

I had up-to-date icd stable version 3.1.2. I have now updated to development version 3.2.0 (fixed typo!).

packageVersion("icd")
[1] ‘3.2.0’
packageVersion("comorbidity")
[1] ‘0.1.1’

As a result, when I re-run Jack's code the 29 element mismatches are resolved:

> testthat::compare(icd_quan, cmb_quan)
Equal

However, there are still these 9 mismatched ids when running Alessandro's code:

> filter(x, same == 0)
     id  code score.icd score.comorbidity same
1    90  V626         2                 1    0
2    90 E8355         2                 1    0
3    90  1907         2                 1    0
4    90 70719         2                 1    0
5    90 66610         2                 1    0
6    90 94239         2                 1    0
7    90 E8375         2                 1    0
8    90 E9909         2                 1    0
9    90  1537         2                 1    0
10   90 62131         2                 1    0
11   90 77084         2                 1    0
12   90 66881         2                 1    0
13   90 71210         2                 1    0
14   90 E8694         2                 1    0
15   90 36463         2                 1    0
16   90  1978         2                 1    0
17  116  V734         4                 3    0
18  116  0213         4                 3    0
19  116 01340         4                 3    0
20  116 79021         4                 3    0
21  116  9839         4                 3    0
22  116 01895         4                 3    0
23  116 20720         4                 3    0
24  116  1100         4                 3    0
25  116 33909         4                 3    0
26  116 25063         4                 3    0
27  116 35981         4                 3    0
28  116 65841         4                 3    0
29  116 77439         4                 3    0
30  116  1319         4                 3    0
31  116   470         4                 3    0
32  116 20246         4                 3    0
33  116 25030         4                 3    0
34  116 36833         4                 3    0
35  116  1762         4                 3    0
36  116 E9093         4                 3    0
37  287 05882         5                 4    0
38  287  2329         5                 4    0
39  287  7905         5                 4    0
40  287 64812         5                 4    0
41  287 E8717         5                 4    0
42  287 01092         5                 4    0
43  287  6145         5                 4    0
44  287  1985         5                 4    0
45  287 74600         5                 4    0
46  287  1412         5                 4    0
47  287  9791         5                 4    0
48  287 66811         5                 4    0
49  287 66710         5                 4    0
50  287 29630         5                 4    0
51  287 42611         5                 4    0
52  287 36275         5                 4    0
53  287 71640         5                 4    0
54  287 33510         5                 4    0
55  287  5716         5                 4    0
56  589  1649         3                 2    0
57  589  1968         3                 2    0
58  589 71914         3                 2    0
59  589 V0261         3                 2    0
60  589  2134         3                 2    0
61  589 65283         3                 2    0
62  589  9931         3                 2    0
63  589 V9122         3                 2    0
64  589 E8860         3                 2    0
65  589 13100         3                 2    0
66  589 34831         3                 2    0
67  589 71885         3                 2    0
68  589 V5041         3                 2    0
69  589 71902         3                 2    0
70  589 94400         3                 2    0
71  589 36600         3                 2    0
72  589 87379         3                 2    0
73  589 73719         3                 2    0
74  589  9107         3                 2    0
75  604  1611         2                 1    0
76  604  0088         2                 1    0
77  604 80496         2                 1    0
78  604 E8780         2                 1    0
79  604 94892         2                 1    0
80  604  0213         2                 1    0
81  604 71805         2                 1    0
82  604  2110         2                 1    0
83  604  1966         2                 1    0
84  604 80043         2                 1    0
85  604 29910         2                 1    0
86  604  7901         2                 1    0
87  604  6310         2                 1    0
88  604  9783         2                 1    0
89  604 43884         2                 1    0
90  604 09381         2                 1    0
91  659  2863         6                 5    0
92  659 29535         6                 5    0
93  659 01580         6                 5    0
94  659 74357         6                 5    0
95  659 37255         6                 5    0
96  659  4572         6                 5    0
97  659  9999         6                 5    0
98  659 80844         6                 5    0
99  659 E0111         6                 5    0
100 659 80330         6                 5    0
101 659 53391         6                 5    0
102 659 E8300         6                 5    0
103 659  7964         6                 5    0
104 659 25080         6                 5    0
105 659  1978         6                 5    0
106 659 E9559         6                 5    0
107 659  6185         6                 5    0
108 659 31500         6                 5    0
109 659  1580         6                 5    0
110 659  2290         6                 5    0
111 659 36053         6                 5    0
112 659 83819         6                 5    0
113 659 V8904         6                 5    0
114 659 65263         6                 5    0
115 659 65920         6                 5    0
116 659 67684         6                 5    0
117 659 V5863         6                 5    0
118 659  6982         6                 5    0
119 659 51661         6                 5    0
120 659  7241         6                 5    0
121 752 40402         5                 4    0
122 752 78901         5                 4    0
123 752  2588         5                 4    0
124 752 77081         5                 4    0
125 752  1977         5                 4    0
126 752 01162         5                 4    0
127 752 53510         5                 4    0
128 752  V616         5                 4    0
129 752 20203         5                 4    0
130 752 71892         5                 4    0
131 752 V1046         5                 4    0
132 752  1915         5                 4    0
133 752 20143         5                 4    0
134 752 41032         5                 4    0
135 752  9925         5                 4    0
136 752  6861         5                 4    0
137 752  V442         5                 4    0
138 752 E9315         5                 4    0
139 752 66410         5                 4    0
140 762  1338         2                 1    0
141 762 56989         2                 1    0
142 762 73382         2                 1    0
143 762  0782         2                 1    0
144 762 46411         2                 1    0
145 762  6948         2                 1    0
146 762  8940         2                 1    0
147 762 V5302         2                 1    0
148 762 90289         2                 1    0
149 762  2700         2                 1    0
150 762  1490         2                 1    0
151 762  9739         2                 1    0
152 762  1972         2                 1    0
153 762 53641         2                 1    0
154 996 30189         2                 1    0
155 996  1480         2                 1    0
156 996 53431         2                 1    0
157 996 V1048         2                 1    0
158 996 01565         2                 1    0
159 996  1981         2                 1    0
160 996 20975         2                 1    0
161 996 36459         2                 1    0
162 996 76418         2                 1    0
163 996 87332         2                 1    0
164 996 E8710         2                 1    0
165 996 71190         2                 1    0
166 996 59652         2                 1    0
167 996 45919         2                 1    0
168 996 63612         2                 1    0
169 996 75249         2                 1    0
170 996 37714         2                 1    0
171 996 36308         2                 1    0
172 996 E8341         2                 1    0
> 
> 
> table(singledf$same)

  0   1 
  9 991 

I do not understand what the "element mismatches" were but that's okay. More important, do you still have the 9 mismatched ids listed above?

[aalexandersson]

I agree with Jack that it is now time to focus on the scoring:

In Alessandro's code, I do not understand the score.comorbidity for the remaining 9 problematic ids. For example, id 90 has "Metastatic cancer" (metacanc = 1) and "Solid tumor without metastasis" (solidtum = 1) in comorbidity9 which implies that score.comorbidity should be 2, but it is 1.

This can be verified by running Alessandro's code until this line:

select(id.comorbidity, score.comorbidity) .

There, run this code:

> filter(comorbidity9, id == 90)
  id chf carit valv pcd pvd hypunc hypc para ond cpd diabunc diabc hypothy rf ld
1 90   0     0    0   0   0      0    0    0   0   0       0     0       0  0  0
  pud aids lymph metacanc solidtum rheumd coag obes wloss fed blane dane alcohol
1   0    0     0        1        1      0    0    0     0   0     0    0       0
  drug psycho depre score.comorbidity index wscore windex id.comorbidity
1    0      0     0                 1   1-4     14    >=5             90

[aalexandersson]

Here is a list of all 9 observations for comorbidity that have problematic scoring to me:

> filter(comorbidity9, id %in% c(90,116,287,589,604,659,752,762,996)) %>%
+    select(id, hypc,ond,diabunc,diabc,rf,ld,pud,lymph,metacanc,solidtum,
+             coag,psycho,depre, score.comorbidity)
   id hypc ond diabunc diabc rf ld pud lymph metacanc solidtum coag psycho depre
1  90    0   0       0     0  0  0   0     0        1        1    0      0     0
2 116    0   0       1     1  0  0   0     1        0        1    0      0     0
3 287    0   1       0     0  0  1   0     0        1        1    0      0     1
4 589    0   1       0     0  0  0   0     0        1        1    0      0     0
5 604    0   0       0     0  0  0   0     0        1        1    0      0     0
6 659    0   0       0     1  0  0   1     0        1        1    1      1     0
7 752    1   0       0     0  1  0   0     1        1        1    0      0     0
8 762    0   0       0     0  0  0   0     0        1        1    0      0     0
9 996    0   0       0     0  0  0   0     0        1        1    0      0     0
  score.comorbidity
1                 1
2                 3
3                 4
4                 2
5                 1
6                 5
7                 4
8                 1
9                 1

[ellessenne]

Hello @aalexandersson!

The difference in score is due to the fact that comorbidity by default does not count the same comorbidity with different degree of severity twice when computing the score. This is controlled by the assign0 argument, which is set to TRUE by default. In the example above, for instance, the individual with id == 116 has both diabetes and diabetes with complications. If we consider both, the score would be 2; however, we would be counting diabetes twice!

See below, calling comorbidity with assign0 = FALSE yields exactly the same result!

library(comorbidity)
library(icd)
#> Welcome to the 'icd' package for finding comorbidities and interpretation of ICD-9 and ICD-10 codes.
#> ?icd to get started, then see the vignettes and help for details and examples.
#> 
#> Suggestions and contributions are welcome at https://github.com/jackwasey/icd . Please cite this package if you find it useful in your published work citation(package = "icd")
requireNamespace("testthat")
#> Loading required namespace: testthat
library(tidyverse)
#> -- Attaching packages --------------------------------- tidyverse 1.2.1 --
#> v ggplot2 2.2.1     v purrr   0.2.4
#> v tibble  1.4.2     v dplyr   0.7.4
#> v tidyr   0.8.0     v stringr 1.3.1
#> v readr   1.1.1     v forcats 0.3.0
#> -- Conflicts ------------------------------------ tidyverse_conflicts() --
#> x dplyr::explain() masks icd::explain()
#> x dplyr::filter()  masks stats::filter()
#> x dplyr::lag()     masks stats::lag()
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)
x <- x[order(x[["id"]]), ]
row.names(x) <- NULL

# Package icd
icd_quan <- icd9_comorbid_quan_elix(x, hierarchy = FALSE)
icd_quan <- data.frame(id.icd = row.names(icd_quan), score.icd = rowSums(icd_quan), 
  stringsAsFactors = FALSE)

# Package comorbidity
comorbidity9 <- comorbidity(x = x, id = "id", code = "code", score = "elixhauser_icd9", assign0 = FALSE) %>% 
  mutate(id.comorbidity = as.character(id)) %>% 
  rename(score.comorbidity = score) %>% 
  select(id.comorbidity, score.comorbidity)

# Merge and id those that differ
singledf <- full_join(icd_quan, comorbidity9, by = c(id.icd = "id.comorbidity")) %>% 
  mutate(same = as.numeric(score.icd == score.comorbidity)) %>% 
  rename(id = id.icd)
x[["id"]] <- as.character(x[["id"]])
x <- left_join(x, singledf, "id")

# Show
filter(x, same == 0)
#> [1] id                code              score.icd         score.comorbidity
#> [5] same             
#> <0 rows> (or 0-length row.names)

[jackwasey]

Hello, as I mentioned before: icd has functions charlson and van_walraven to calculate scores which take these things into account and weight them in validated ways. If you just sum the columns of the icd results, yes, you will count related things twice. In fact, different comorbidity schemes have different expectations when thinking about this. icd has the hierarchy argument to the comorbidity functions so that each set of rules can be applied. See the internal functions named like apply_hier_ahrq in R/comorbid.R in the icd package. If you want to apply your own rules, which do merge e.g. hypertension with and wihout complications, then just or the HTN columns before row summing. I have occasionally seen people do what you are doing, and just sum the number of comorbidities, but most researchers who even use scores will use Charlson or Van Walraven in my experience. Or, of course, calculate a propensity score from other features, too.

I'm glad at last we have verified the calculations!

[jackwasey]

This has been a really helpful discussion, so can we go back to the original problem now? Let me know when you've read the discussion about interpretation of incorrect three-digit codes in my JSS article.

[ellessenne]

@jackwasey comorbidity computes the simple sum of comorbidities (with hierarchy depending on the assign0 argument), and also the Charlson and van Walraven scores (with weights and hierarchy via assign0 as well)!

I also read the draft of your JSS paper (very well written by the way, good job): summarising your main points,

1. three-digits codes that are not valid are discarded by default;

2. children codes (with 4 digits or more) are assigned, even if invalid.

Please correct me if I am mistaken.
While I do agree with what you have implemented (and in fact comorbidity does point (2) as well), I think this is still not that big of a deal: assuming only valid ICD codes are supplied, no difference should be observed. Digging into the quality of the input data is a whole different thing, with great implications: what if mistyped codes actually exist, e.g. in a different domain?

I have some ideas on validation (well, sort of) to play around and potentially implement in comorbidity. For instance, what if we throw warnings about invalid codes instead?

Alessandro

P.s.: it was very interesting to read about the performance of icd! I had similar ideas that I am planning to implement to make comorbidity faster, I am looking forward to a new comparison (hopefully soon)!

[aalexandersson]

Thank you both. The official Stata command icd might give you some more ideas for validation.

[ellessenne]

Thanks @aalexandersson, I was not aware of Stata's icd!

[aalexandersson]

@jackwasey @ellessenne
The new version 0.3.0 of comorbidity can calculate the 29 "AHRQ comorbidities". Like we did for the 31 Quan comorbidities, I would like to be able to calculate the same unweighted sum scores using your packages icd version 3.3 and comorbidity version 0.3.0, respectively for AHRQ comorbidities. I get different AHRQ results when using regular unweighted sum scores for icd and comorbidity, how come?

Here is the same test data as before:

library(comorbidity)
library(icd)
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)
x <- x[order(x[["id"]]), ]
row.names(x) <- NULL
save(x,file="simulated.Rda")

Here are the AHRQ sum scores from icd:

load("simulated.Rda")
v_cmb <- icd9_comorbid_ahrq(x) %>%
   apply(2, as.integer) # convert logical to integer
elixsum <- rowSums(v_cmb, na.rm = TRUE)
table(elixsum)

This is the output:

elixsum
  0   1   2   3   4   5   6   7 
211 363 251 123  46   4   1   1

Here are the AHRQ sum scores from comorbidity:

elixhauser9 <- comorbidity(x = x, id = "id", code = "code", score = "elixhauser", 
  icd = "icd9", assign0 = TRUE, factorise = TRUE)  %>%
  mutate(id.comorbidity = as.character(id)) %>% 
  rename(score.comorbidity = score) %>% 
  rename(wscore.comorbidity = wscore_ahrq) %>% 
  rename(vw.comorbidity = wscore_vw) %>% 
  select(id.comorbidity, score.comorbidity, wscore.comorbidity, vw.comorbidity)

table(elixhauser9$score.comorbidity)

This is the output:

  0   1   2   3   4   5   6 
245 359 235 113  36  10   2 

Why do the two sets of sum scores differ for AHRQ? How can the two programs calculate the same sum scores for AHRQ?

[ellessenne]

Hi @aalexandersson, could you try running the same code with assign0 = FALSE - the new default for comorbidity?

[aalexandersson]

Here is the output when instead using syntax assign0 = FALSE. There is still a huge difference:

  0   1   2   3   4   5   6 
245 355 238 113  35  11   3 

[ellessenne]

I see. Let me have a deeper look, I'll get back to you ASAP!

[aalexandersson]

For regular (unweighted) sum scores, I think we should use assign0 = TRUE as we discussed and did for sum scores of Quan comorbidities.

[ellessenne]

@aalexandersson I had a quick look at some of the differences and I think it may depend on which codes have been used to detect the comorbidities; I am overall using the "Enhanced ICD9-CM" from Quan et al., not sure what @jackwasey is using though.
Regarding the hierarchy of assign0, I had a discussion via e-mail with some researchers in California using comorbidity scores routinely and I was convinced of switching the default to FALSE. I thought about it a bit more now and I actually think I will change again to no default, "forcing" the user to decide whether to apply a hierarchy or not. In this way, every user is aware of what the algorithm is doing behind the scenes. What do you think?

[aalexandersson]

I like your idea of having no default for assign0 since the choice is fundamental. Whatever you decide, good documentation is a must. Thanks!

[aalexandersson]

For AHRQ comorbidities, it seems to me that @jackwasey in icd instead is using the "official AHRQ minus DRG", not the "Enhanced ICD9-CM" by Quan.

I base this conclusion from reading jackwasey/icd#112 and https://cran.r-project.org/web/packages/icd/vignettes/compare-maps.html . I guess that explains the differences in sum scores between comorbidity and icd for AHRQ comorbidities.

[aalexandersson]

I just remembered that for icd sumscores, we need the option hierarchy = FALSE. Here is the output:

  0   1   2   3   4   5   6   7 
211 357 255 124  42   8   1   2

[jackwasey]

Thanks for bringing this up and you both looking into it. I will have some time soon to dig into this.

[aalexandersson]

I used icd9_comorbid_ahrq(); I think that was wrong. My bad.
I now think this is correct about AHRQ scoring:

• The (31) Quan Elixhauser comorbidities: The packages comorbidity and icd can both calculate it. The sum scores were the same in the example. For ICD-9, the equivalent function in package icd is icd9_comorbid_quan_elix(). The new version 0.3.0 of comorbidity adds AHRQ scoring. A comparison of scoring between the two packages is needed.

• The (29) AHRQ Elixhauser comorbidities: The package icd can calculate it; the package comorbidity cannot. Therefore, a comparison of scoring between the two packages is not possible.

I think this is the remaining issue:
Are the AHRQ scoring of Quan Elixhauser comorbidities the same or different in the two packages?

I will try to look into this next week. Any feedback is welcome. After all, you both are the experts! At least, I am learning :-)

[aalexandersson]

Here is yet another attempt to be more precise about "AHRQ Elixhauser" comorbidities and its AHRQ scoring. Again, this is simply my understanding of what you both developed:

The (29) AHRQ Elixhauser comorbidities:
-The package icd: It uses the current definition of "AHRQ-Web" with modifications such as not including DRGs and "other modifications" listed on the icd website. It would be helpful to have a concise coding definition similar to "AHRQ-Web" in Quan et al. (2005, Table 2). Does it exist, or what are the modifications other than not using DRGs? What AHRQ scoring is used?

-The package comorbidity: It uses 29 of the 31 "Enhanced ICD-9-CM" comorbidities in Quan et al. (2005, Table 2) . Cardiac arrythmias is not included in AHRQ and therefore it has scoring weight 0. There is only one hypertension comorbidity (combined uncomplicated and complicated), and it has scoring weight -1. The used AHRQ scoring is defined in Moore et al. (2017). There is a helpful vignette: vignette("comorbidityscores", package = "comorbidity").