nlpembeds

Co-occurrence and PMI-SVD

Background

Pointwise Mutual Information (PMI) and Singular Value Decomposition (SVD) are methods to analyse co-occurrence matrices. Following word2vec's gain in popularity in the year 2010s, it has been proven that word2vec produces an approximation of the following pipeline: 1) computation of co-occurrences in a sliding window 2) computation of single shift PMI (SSPMI) 3) SVD.

Word2vec enables the approximation to be computed in a very parallelized way, however it may lack stability, i.e. the same computation will produce significantly different results. To improve the stability, one can increase the number of iterations, but it may be preferable to compute co-occurrence, SSPMI, and SVD.

The nlpembeds package aims to provide efficient methods to compute co-occurrence matrices and PMI-SVD. In the biomedical and clinical settings, one challenge is the huge size of databases, e.g. when analyzing data of millions of patients over tens of years. To address this, the nlpembeds package provides functions to efficiently compute monthly co-occurrence matrices, which is the computational bottleneck of the analysis, by using the RcppAlgos package and sparse matrices. Furthermore, the functions can be called on SQL databases, enabling the computation of co-occurrence matrices from tens of GBs of data, representing millions of patients over tens of years.

Co-occurrence matrix

Consider the following data object:

library(nlpembeds)
df_ehr = data.frame(Month = c(1, 1, 1, 2, 2, 3, 3, 4, 4),
                    Patient = c(1, 1, 2, 1, 2, 1, 1, 3, 4),
                    Parent_Code = c('C1', 'C2', 'C2', 'C1', 'C1', 'C1', 'C2',
                                    'C3', 'C4'),
                    Count = 1:9)
df_ehr

This represents our patient-level data: Patient is the patient identifier, Parent_Code is the diagnosis that was performed (or medication prescribed, etc.), Month is the year/month of the diagnosis, and Count is the number of times the diagnosis was performed during that month.

We can compute a monthly co-occurrence which gives us:

spm_cooc = build_df_cooc(df_ehr)
spm_cooc

What the data represents is:

• At month 1, patient 1 had the code C1 once and C2 twice, patient 2 had the code C2 thrice.
• At month 2, patients 1 and 2 had the code C1 4 and 5 times, respectively.
• At month 3, patient 1 had the codes C1 and C2 6 and 7 times, respectively.
• At month 4, patients 3 and 4 had the codes C3 and C4 8 and 9 times, respectively.

The monthly co-occurrence between two codes is defined as the minimum of their two counts within a month. The co-occurrences are first computed for each patient for each month, and sum aggregated.

Here, let's consider the co-occurrence between codes C1 and C2. We decompose the computation to make it easily understandable. The codes co-occurred in months 1 and 3 in patient 1. We can decompose the co-occurrence matrices as follows:

• The codes co-occurred min(1, 2) = 1, in Month 1

  cooc_1 = build_df_cooc(subset(df_ehr, Patient == 1 & Month == 1), min_term_freq = 0)
  cooc_1

• The codes co-occurred min(6, 7) = 6, in Month 3

  cooc_2 = build_df_cooc(subset(df_ehr, Patient == 1 & Month == 3))
  cooc_2

When summed up, we get a co-occurrence of 7 between C1 and C2

  cooc_1 + cooc_2

Note that although codes C3 and C4 were both present during month 4, they were not for the same patient, so their co-occurrence is 0.

PMI

We can then call PMI on that object

m_pmi = get_pmi(spm_cooc)
m_pmi

Here,

• PMI(C1, C2) is log(7 * 59 / ((16 + 7) * (12 + 7)))
• PMI(C1, C1) is log(16 * 59 / (16 + 7) ^ 2)

SVD

We then call SVD on the PMI matrix, and then feed the m_svd object (i.e. the embedding matrix) to downstream analyses, such as knowledge graphs.

m_svd = get_svd(m_pmi, svd_rank = 2)
m_svd

The SVD is computed with "random SVD", an efficient approximation of truncated SVD, in which only the first principal components are returned. It is computed with the rsvd package, and the author suggests that the number of dimensions requested k should follow: k < n / 4, where n is the number of features, for it to be efficient, and that otherwise one should rather use either SVD or truncated SVD.

I recommend a SVD rank between 1000 and 2000 when analyzing codified + NLP (15k - 25k features), and less if only analyzing codified (300 - 800 dimensions for 5k - 10k codified features). One can use known pairs to get the AUC and select the best performing rank, e.g. with the kgraph app or the kgaph R package. The aim is to identify the "shoulder" after which increasing the number of dimensions does not significantly increase the AUC anymore. For generalization and reproducibility within other cohorts, it may also be better to have a slightly smaller AUC and a lower number of dimensions, as a large number of dimensions can lead to overfitting and a lower number of dimensions will act as a regularization.

In my opinion, the "edge cases" would be if for example the AUC is 0.732 at 1000 dimensions, and 0.75 at 2000 dimensions. In that case either way would be fine and perhaps a more thorough investigation should be performed to choose between 1000 or 2000 (or just go with 1500). This is just an example on how one can balance the 2 objectives: higher AUC and lower number of dimensions.

Out-of-memory SQL databases

To efficiently perform co-occurrence matrices of large databases of tens of GBs, two important optimizations are possible:

• Batching by patients
• Code subsets based on dictionaries

Batching by patients

To perform batching by patients, we want to have our input file as a SQL file, which we will read by batches of patients and aggregate the results. This is performed with the sql_cooc function.

To demonstrate how it is used, we first write the example object above as a SQL database. It is important to name the columns and the SQL table as here (i.e. four columns Patient, Month, Parent_Code, Count and table name df_monthly). Also, we index the table by patients, and write as a second table df_uniq_codes the unique codes (this is optional, it is done automatically in sql_cooc if the table df_uniq_codes is not found in the SQL file). It is a good idea to also order your table by patients before writing it, as I believe it will make read accesses more efficient, though it is not required and should not impact performance too much I believe.

library(RSQLite)

test_db_path = tempfile()
test_db = dbConnect(SQLite(), test_db_path)
dbWriteTable(test_db, 'df_monthly', df_ehr, overwrite = TRUE)

###
# optional, done automatically by sql_cooc if table df_uniq_codes not found
df_uniq_codes = unique(df_ehr['Parent_Code'])
dbWriteTable(test_db, 'df_uniq_codes', df_uniq_codes, overwrite = TRUE)

dbExecute(test_db, "CREATE INDEX patient_idx ON df_monthly (Patient)")
###

dbDisconnect(test_db)

We can then call the sql_cooc function on that filename, and specify a new filename for the output.

For each batch, information about the current memory usage is printed, which can be helpful for debugging memory usage and getting an idea of the computation progress from the log file.

output_db_path = tempfile()
sql_cooc(dbpath = test_db_path, outpath = output_db_path, select_cuis = 'all',
         n_batch = 2)

Once the computation is complete, we read the co-occurrence sparse matrix from the output SQL file.

test_db = dbConnect(SQLite(), output_db_path)
spm_cooc = dbGetQuery(test_db, 'select * from df_monthly;')
dbDisconnect(test_db)

spm_cooc

As previously, we can then feed it to get_pmi.

  m_pmi = get_pmi(spm_cooc)

Or transform it as a classic matrix.

  spm_cooc = build_spm_cooc_sym(spm_cooc)                                       
  m_cooc = as.matrix(spm_cooc)                                       

That is the minimum call. Two important parameters come into play here:

• n_batch which is the number of patients to include by batch
• n_cores which is the number of cores on which the computation will be parallelized

For the two parameters, the higher the value, the faster the computation and the more RAM memory required. The user can fine-tune these parameters to fit machine specifications and to optimize. As an example, n_batch = 300 and n_cores = 6 should be able to run on 16GB machines if the number of codes is not too large (e.g. if only considering rolled-up codified data).

Code subsets based on dictionaries

When the number of unique codes is large, e.g. when analyzing natural language processing (NLP) concepts, one may want to perform a first subset of the codes by providing a dictionary of codes to include.

Two parameters are available here:

• select_cuis which has three possible values: none (to discard all NLP concepts), mh_5k to use the file nlp_dict_cuis_5188.csv, and mh_20k to use the union of previous file and the file general_mental_health_nlp_dict_15k_cuis.tsv.
• mh_cuis_dirpath which is the folder where the NLP dictionaries are located

The behavior of these parameters will be improved and generalized later on to enable the user to specify personal dictionaries. For the time being the user can just reuse those filenames and change the files.

To demonstrate the behavior, we will rename the two first codes in the object above to two NLP concepts in those dictionaries.

df_ehr$Parent_Code %<>% ifelse(. == 'C1', 'C0000545', .)
df_ehr$Parent_Code %<>% ifelse(. == 'C2', 'C0000578', .)

df_ehr

We write it as a SQL file.

test_db = dbConnect(SQLite(), test_db_path)
dbWriteTable(test_db, 'df_monthly', df_ehr, overwrite = TRUE)

df_uniq_codes = unique(df_ehr['Parent_Code'])
dbWriteTable(test_db, 'df_uniq_codes', df_uniq_codes, overwrite = TRUE)
dbDisconnect(test_db)

The code is then called like this:

sql_cooc(dbpath = test_db_path, outpath = output_db_path,
         select_cuis = 'mh_20k', n_batch = 2)

We then read the co-occurrence sparse matrix from that SQL file.

test_db = dbConnect(SQLite(), output_db_path)
spm_cooc = dbGetQuery(test_db, 'select * from df_monthly;')
dbDisconnect(test_db)

spm_cooc

Running on HPC servers as O2

Installation

You need to first install v2.4.0 of RcppAlgos, then the nlpembeds package.

remotes::install_git('https://github.com/jwood000/RcppAlgos@v2.4.0.git')
remotes::install_git('https://gitlab.com/thomaschln/nlpembeds.git')

Parameters tuning

Assuming 350 GB are available, I would recommend n_batch = 75 and n_cores = 25, and estimate the co-occurrence matrix for ~60k patients should take 3 hours (assuming 20k CUIs considered). You should be able to more or less multiply the parameters by the amount of RAM, meaning if you have 700 GB available you could use n_batch = 150 and n_cores = 50 and it should run in ~1h30.

The number of CUIs have a quadratic impact on memory / computation time, so if not considering any CUIs you can significantly increase n_batch (several thousands should work). If only considering 5k CUIs, you can most likely have n_batch = 600. If considering more than 20k CUIs, you will need to reduce n_batch and/or n_cores.

The n_batch patients are parallelized on the n_cores available, so it is of no use to give more n_cores than n_batch, and it is best to have at least 2-3 patients per core. The aggregated sparse matrix will logarithmically increase in size, so if you are close to the RAM limits your job may fail only after several hours.

Ontologies

Bioportal datasets

Populate the 'ontologies' folder in 'inst'

wget -O ontologies.7z "https://www.dropbox.com/scl/fi/10xu9nvr5pb52tvklz9ct/ontologies.7z?rlkey=w3a2n4uitrrfhuwktt5hx89zw&st=og0mjxvm&dl=0"

Downloaded from https://bioportal.bioontology.org/

Raw and propagated graphs

Copy to 'data' folder, visualize with shiny app cui_onto

wget -O onto_graph1.rds "https://www.dropbox.com/scl/fi/wmb58g6vmjecmg5ay089j/onto_graph1.rds?rlkey=44azfaun7sg5tazbnz0333vjw&st=1qh5rxve&dl=0"

wget -O onto_graph2.rds "https://www.dropbox.com/scl/fi/od12ghd0xto3o2quo2j1b/onto_graph2.rds?rlkey=qklrzzrmubbu8nilu7nef45kd&st=ifo6pmrd&dl=0"

EPMC test files

wget https://europepmc.org/ftp/oa/PMC1126903_PMC1129089.xml.gz
wget https://europepmc.org/ftp/oa/PMC1743054_PMC1747532.xml.gz