<jupyter_start><jupyter_text>Latent Dirichelet Allocation on Clinical Notes Author: <NAME>This script builds the Latent Dirichelet Allocation (LDA) models on the MIMIC-III clinical notes. In order to manage the large number of clinical notes (over 2 Million) we choose to subsample these notes. We began with modestly sized training and test samples of 1000 clinical notes, however determined that a larger sample was necessary to appropriately capture the variation in this dataset.We selected a subsample size of 100,000 clinical notes to develop the Topic Model. We used perplexity analysis (described below) to select our hyperparameters. The Date-Time ShiftThe dates and times have been shifted to preserve patient anonymity. The relative dates for a given patient are conserved, however the dates between patients or relative to the chronological ordering cannot be inferred. In order to appropriately generate the training and test sets, we need to perform a within-patient hold-out set. That is to say, we hold-out the data for the LAST visit to validate our methods. [MIMIC Date-Time Shifting](https://mimic.physionet.org/mimicdata/time/)<jupyter_code>from nltk.tokenize import RegexpTokenizer
from stop_words import get_stop_words
from nltk.stem.porter import PorterStemmer
from gensim import corpora, models
import gensim
import random
import time

# MIMIC-III Connection Imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import psycopg2
# ----------------------------

# Setup some Global Variables
VERBOSE = True
PICKLE_RICK = True


# Set the database connection information
sqluser = 'postgres'
dbname = 'mimic'
schema_name = 'mimiciii'

# Connect to postgres with a copy of the MIMIC-III database
con = psycopg2.connect(dbname=dbname, user=sqluser)

# the below statement is prepended to queries to ensure they select from the right schema
query_schema = 'set search_path to ' + schema_name + ';'

if VERBOSE: print "Setup Complete!"
tokenizer = RegexpTokenizer(r'\w+')

# create English stop words list
en_stop = get_stop_words('en')
if VERBOSE: print en_stop
# Create p_stemmer of class PorterStemmer
p_stemmer = PorterStemmer()
# Determine the average number of clinical notes per patient...
# This will help us determine the number of clinical notes to subsample when creating the LDA topic model!
query = query_schema + """
SELECT subject_id, text
FROM noteevents
"""
df = pd.read_sql_query(query, con)
total_notes = len(df)

query = query_schema + """
SELECT subject_id, COUNT(subject_id) AS occurences
FROM noteevents
GROUP BY subject_id
"""
df = pd.read_sql_query(query, con)
avg_notes = df["occurences"].mean()
subsample_size = total_notes / avg_notes
print 'Average number of Clinical Notes per Patient: ' + str(avg_notes)
print 'Approximation of Minimum Subsample Size: ' + str(subsample_size)
df.head()<jupyter_output><empty_output><jupyter_text>Selected Subsample Size: 100,000The average number of clinical notes per patient approximates the number of clinical notes. By setting the subsample size to 100,000 we can use ~2 clinical notes per patient which is useful when separating into training and testing subsamples.<jupyter_code># === RANDOMLY SUBSAMPLE CLINICAL NOTES ===
query = query_schema + """
SELECT subject_id, text
FROM noteevents
ORDER BY random()
LIMIT 100000
"""
df = pd.read_sql_query(query, con)

print 'Subsampled Size of Clinical Notes: ' + str(len(df))
# Compile the clinical notes into a list
doc_set = []
print 'DataFrame Size: ' + str(len(df))
for i in range(0, len(df)):
    doc_set.append(df.at[i,'text'])

print 'Length of Doc_set: ' + str(len(doc_set))
#print doc_set[2] # Look at the Second Entry to Validate
#doc_set = [doc_a, doc_b, doc_c, doc_d, doc_e]<jupyter_output>DataFrame Size: 100000
Length of Doc_set: 100000<jupyter_text>Pre-Processing of Clinical NotesPrior to dumping the files into the LDA, we need to perform a number of pre-processing steps.1) **Tokenization**: The text must be chopped into individual units ("tokens") for independent processing2) **Stop Word Removal**: Stop words are words which are filtered out before or after processing of natural language data and comprise the most common words in the English language.3) **Stemmitization**: The process of reducing inflected words to their word stem, base, or root form. (*e.g.* stemming --> stem)---An important additional step to be incorporated in future work is the removal of negated terms:[NegEx](https://github.com/mongoose54/negex/tree/master/negex.python)<jupyter_code># list for tokenized documents in loop
texts = []

# Variables to Track Summary Stats
token_numbers = []
stop_numbers  = []
stem_numbers  = []

# Now loop through document list
for i in doc_set:
    
    # clean and tokenize document string
    raw = i.lower()
    tokens = tokenizer.tokenize(raw)
    token_numbers.append(len(tokens))
    
    # remove stop words from tokens
    stopped_tokens = [i for i in tokens if not i in en_stop]
    stop_numbers.append(len(stopped_tokens))
    
    # stem tokens
    stemmed_tokens = [p_stemmer.stem(i) for i in stopped_tokens]
    stem_numbers.append(len(stemmed_tokens))
    
    # add tokens to list
    texts.append(stemmed_tokens)
    
# Get a sense of sizes:
print 'Average Number of Tokens per Document:  ' + str(sum(token_numbers) / float(len(token_numbers)))
print 'Average Number after Stop & Stem per Document: ' + str(sum(stem_numbers) / float(len(stem_numbers)))
print 'Number of Tokens over all Documents: ' + str(len(texts))
# === OPTIONAL ===
# Saving the list of tokens to file for relaoding later (previous computation time was several hours)
f = open('./saved_tokens/token_list_100000_Subset.txt', 'w')
for t in texts:
    f.write(str(t) + '\n')
f.close()
# === OPTIONAL ===
# Convert the tokenized documents into a Key:Value Dictionary where Key = Id; Value = Term
dictionary = corpora.Dictionary(texts)
#if VERBOSE: print dictionary
# Convert tokenized documents into a document-term matrix
# This is a matrix where the 
bow_corpus = [dictionary.doc2bow(text) for text in texts]
#if VERBOSE: print corpus
    
# === HERE WE ORGANIZE THE TRAIN AND TEST CORPI ===
# split into train and test - random sample, but preserving order
import random
train_size = int(round(len(bow_corpus)*0.8))
train_index = sorted(random.sample(xrange(len(bow_corpus)), train_size))
test_index = sorted(set(xrange(len(bow_corpus)))-set(train_index))
train_corpus = [bow_corpus[i] for i in train_index]
test_corpus = [bow_corpus[j] for j in test_index]

print 'Training Corpus Size :: ' + str(len(train_corpus))
print 'Testing Corpus Size :: ' + str(len(test_corpus))<jupyter_output>Training Corpus Size :: 80000
Testing Corpus Size :: 20000<jupyter_text>Perplexity AnalysisWith our subsampled set of clinical notes, we run perplexity analysis to determine the best number of topics to appropriately model our patients.The most common way to evaluate a probabilistic model is by training a the model on a corpus and measuring the log-likelihood of the held-out test set. The test set comprises a collection of unseen documents, in this case, the clinical notes of an independent patient. Given the naive subsampling approach, there is a small risk that the test corpus contains a clinical note pertaining to a patient with another note in the training set. This does not constitute duplication, however is a less stringent form of model generation and is suitable to our purpose given that the occureent is negligible. The test set is a collection of unseen documents, $w_{d}$.The model is described by the topic matrix $\Phi$ and the hyper parameter $\alpha$ for topic distribution of documents.The log-likelihood is evaluated as:$$ \mathcal{L}(w) = log p(w|\Phi, \alpha) = \sum_{d}{}log p(w_{d} | \Phi, \alpha) $$The higher the likelihood value, the better the model is considered. We can compare different topic models by considering their likelihood of unseen documents. This convention for measuring the suitability of a topic model is known as the *perplexity* of the held-out documents:$$ perplexity(w) = e^{-\frac{\mathcal{L(w)}}{||tokens||}}$$This is a decreasing function of the log-likelihood $\mathcal{L}(w)$ over the unseen documents $w_{d}$ and a lower perplexity value is indicative of a better model since it indicates less **surprise** when evaluatating the unseen documents. A larger perplexity value indicates a larger degree of **surprise**.---  Parameter Tuning over the Training and Test SetHere, we prepare the training and test set using the prepared bag-of-words corpus of all subsampled clinical notes. We will iterate over a parameter range from 10 to 300 in steps of 10, for 30 evaluations, each with 10 iterations. We will then plot the results and identify our hyperparameters.<jupyter_code># Prepare the Variables for looping...
data_path = '/home/kevindick/MIMIC-III/code/latent_dirichelet_allocation/saved_LDA_models/'
num_words = sum(cnt for document in test_corpus for _, cnt in document)
parameter_range = range(10, 310, 10)
results = {}

for num_topics in parameter_range:
    print "Beginning pass for number of topics :: " + str(num_topics)
    start_time = time.time()
    model = models.LdaMulticore(corpus=train_corpus, workers=None, id2word=dictionary, num_topics=num_topics, iterations=10)
    print 'Model Build Time: ' + str(time.time() - start_time)

    start_time = time.time()
    perplexity = model.bound(test_corpus) # The Topic Model perplexity, not the per word perplexity
    print "Topic Model Perplexity: " + str(perplexity)
    print 'Model Test Time: ' + str(time.time() - start_time)
    #grid.get(str(num_topics), []).append(perplexity)
    #print grid[str(num_topics)]
    
    per_word_perplex = np.exp2(-perplexity / num_words)
    print "Per-word Perplexity: %s" % per_word_perplex
    #grid.get(num_topics).append(per_word_perplex)
    model.save(data_path + 'ldaMulticore_NumTopics_' + str(num_topics) + '_training_corpus.lda')
    print 'Saved Model to File!'

    # Append to dictionary of results
    results["numTopics" + str(num_topics)] = [perplexity, per_word_perplex]

# Generate the LDA model
#ldamodel = gensim.models.ldamodel.LdaModel(corpus, num_topics=300, id2word = dictionary, passes=5)
# Here we plot the results...
import cPickle as pickle
import plotly.graph_objs as go
import plotly.plotly as py

# Save the results
#pickle.dump(results, open(data_path + "perplexity_results.p", "wb"))

perplexies = []
for num_topics in parameter_range:
    perplexies.append(results["numTopics" + str(num_topics)][1])
    #print num_topics, '\t',  results["numTopics" + str(num_topics)]

trace = go.Scatter(
    x = parameter_range,
    y = perplexies,
    mode = 'lines+markers',
    name = 'lines+markers'
)    

data = [trace]
layout = dict(title = 'Perplexity Analysis Results',
              yaxis = dict(title = 'Perplexity'),
              xaxis = dict(title = 'Number of Topics'))
fig = dict(data=data, layout=layout)
py.iplot(fig, filename='perplexities-plot')

#df.to_pickle(data_path + 'gensim_multicore_topic_perplexity.df')<jupyter_output><empty_output><jupyter_text>Post-Perplexity Analysis ThoughtsFrom the resulting distribution of perplexity values the following conclusions are drawn:* **Minimum Perplexity Difficult to Define**: We anticipated using a minimum in the distribution as our operating point. No clear minium presented itself in this range and over this corpus. Small deviances of the general trend are seen at 220 and 270.* **Instability at Higher Topic Ranges**: Clearly a large amount of deviance from the original trend is percieved in the upper topic ranges. It is difficult to determine whether trends within these ranges account for meaningful representation of topics. Perhaps a higher resolution of model comparison might elucidate trends in these regions and establish a veritable "dip" indicative of a more appropriate topic value.* **Dip at 290**: The only point in the curve that resulted in a lower perplexity value than the previous iteration was at 290 topics.**Conclusion**: The operating threshold selected at this time is 270 given the small dip prior to the dramatic rise in perplexity and instability of higher topic ranges.<jupyter_code># Re-Load the 270 Topic Model to move forward!
NUM_TOPICS = 270
model = model.load(data_path + 'ldaMulticore_NumTopics_' + str(NUM_TOPICS) + '_training_corpus.lda')<jupyter_output><empty_output><jupyter_text>Inferring Topic Distributions of All Clinical Notes for our Patient Population We select the operational point for the LDA, reload the model, and iterate over all notesWe first segregate our patient population. We consider patients with at least two clinical notes to ensure they have a minimum number of notes to be represented in both the training and test sets.The results of each inference will be saved relative to the patient via **subject_id**. Once all notes have been processed, we need to average the topic distributions of all notes for a given patient to produce a single representative topic distribution. When possible, we withhold one note for each patient for our test set in the subsequent classification phase.<jupyter_code># Make fancy query to pull the clinical notes of only patients with two or more admissions

# Count number of times subject_id appears in the noteevents table (indicates number of unique notes)
# 
query = query_schema + """
SELECT fusion_table.subject_id, fusion_table.chartdate, fusion_table.hadm_id, fusion_table.text
FROM (
    SELECT subject_id, num_visits
    FROM (SELECT subject_id, COUNT(subject_id) AS num_visits
        FROM admissions
        GROUP BY subject_id) AS count_table
    WHERE num_visits > 1) AS visit_table
JOIN noteevents AS fusion_table ON visit_table.subject_id = fusion_table.subject_id 
"""
df = pd.read_sql_query(query, con)
df.head(20)
# Sort the subject_id column to group together the date, and then internally sort the dates for each patient...
# We will then post-process this further to ensure the data relevant to the last visit is excluded!
sorted_df = df.sort_values(['subject_id', 'chartdate'], ascending=[True, True])
sorted_df.head(20)<jupyter_output><empty_output><jupyter_text>Experimental DesignWe have conserved the clinical notes of only those patients with two or more visits. Given the sorted dataframe of subject_id and datetimes relative to the clinical notes, we can appropriately divide this dataset into a training and test set on a per-patient basis. We use the *hadm_id* in the last record to allocate all notes related to that visit to the test set.---We  will now divide the dataset into a training and testing dataset. We differentitate the test set from the training set as the LAST record in the EHR. This  ensure that the clinical note is the most *future* record in the dataset. This same approach is employed accross this study. We use the latest temporal record based in *hadm_id* on a within-subject basis to ensure we predict the future diagnoses of that patient rather than performing simple recall.<jupyter_code># Generate the Training set by removing the last record from each subject within their group of notes
train_df = pd.DataFrame(columns=['subject_id', 'chartdate', 'hadm_id', 'text'])
test_df = pd.DataFrame(columns=['subject_id', 'chartdate', 'hadm_id', 'text'])

# Create array of all unique subject_ids
subjects = sorted_df.subject_id.unique()
print 'Number of Subjects: ' + str(len(subjects))

# Iterate over all subjects and get their grouped data in the sorted_df
for subject in subjects:
    group = sorted_df.loc[sorted_df['subject_id'] == subject]
    
    last_hadm = group.tail(1)['hadm_id']
    
    # List of items we want removed from the 
    to_remove = [last_hadm.item()]#, 'nan', 'NaN']
        
    not_last  = group.hadm_id.unique().tolist()
    not_last  = [item for item in not_last if str(item) not in to_remove]
    
    # Put the data related to the LAST VISIT into the test set
    test_df = pd.concat([test_df, group[group.hadm_id.isin(last_hadm)]])
    
    # Remaining get appended to the training set
    #train_df = pd.concat([train_df, group.loc[(group['hadm_id'] != last_hadm)]])
    train_df = pd.concat([train_df, group[group.hadm_id.isin(not_last)]])

    
print "Training Set Size: " + str(len(train_df))
print "Testing Set Size: " + str(len(test_df))

if PICKLE_RICK:
    pickle.dump(train_df, open(data_path + "train_df.p", "wb"))
    pickle.dump(test_df, open(data_path + "test_df.p", "wb"))
# Now submit each of the training notes to be computed... 
# We leave out the testing notes since these will be used to create the representation of the patient in the testing phase
topic_distributions = {}
row = -1
# Submit each clinical note to the model to obtain the corresponding topic distribution
for note in train_df["text"]:
    start_time = time.time()
    row += 1
    texts = [] # Only want one note at a time so refresh for each note in the loop...
    patient_id = train_df['subject_id'][row]
    patient_hadm = train_df['hadm_id'][row]
    
    # Pre-Process
    raw = note.lower()
    tokens = tokenizer.tokenize(raw)
    stopped_tokens = [i for i in tokens if not i in en_stop]
    stemmed_tokens = [p_stemmer.stem(i) for i in stopped_tokens]
    texts.append(stemmed_tokens)
    dictionary = corpora.Dictionary(texts)
    doc_bow = dictionary.doc2bow(text)
    
    # Now obtain the topic distribution!
    doc_lda = model[doc_bow] 
    
    # Save the results to a new df with the subject_id conserved...
    topic_distributions[str(row)] = [patient_id, patient_hadm, doc_lda]
    if VERBOSE and row < 100: # Only print the first couple to get an idea!
        print "Patient Id   :: " + str(patient_id)
        print "Patient HADM :: " + str(patient_hadm)
        print doc_lda
        print "Note " + str(row) + " :: Topic Compute Time = " + str(time.time() - start_time)

pickle.dump(topic_distribution, open(data_path + "train_df_topic_distributions.p", "wb"))
print "Execution Complete!~"<jupyter_output>Patient Id   :: 17
Patient HADM :: nan
[]
Note 0 :: Topic Compute Time = 0.0442798137665
Patient Id   :: 17
Patient HADM :: 194023.0
[(14, 0.012824726650240037), (116, 0.70188837837473772), (134, 0.019360043611712931), (211, 0.1855474973422605), (241, 0.050430253867526385)]
Note 1 :: Topic Compute Time = 0.0260610580444
Patient Id   :: 17
Patient HADM :: 194023.0
[(251, 0.75092581786602886)]
Note 2 :: Topic Compute Time = 0.00909900665283
Patient Id   :: 17
Patient HADM :: 194023.0
[(14, 0.010761039994930197), (94, 0.029200275790082983), (117, 0.20038510267512843), (128, 0.11103096035874453), (132, 0.20490983783465841), (137, 0.21364957242226087), (142, 0.090140396543946871), (211, 0.1216139394675925)]
Note 3 :: Topic Compute Time = 0.0219399929047
Patient Id   :: 17
Patient HADM :: 194023.0
[(94, 0.4540388551536565), (116, 0.34744262632781858)]
Note 4 :: Topic Compute Time = 0.00867795944214
Patient Id   :: 17
Patient HADM :: 194023.0
[(14, 0.41408765006485787), (130, 0.14701499021372[...]
{'fraction_non_alphanumeric': 0.08689256632519043, 'fraction_numerical': 0.032991857105332284, 'mean_word_length': 4.695990424895272, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 29, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '338790', 'license': 'MIT', 'n_tokens_mistral': 6158, 'n_tokens_neox': 5568, 'n_words': 2228, 'path': 'LDA/Latent_Dirichelet_Allocation_on_Clinical_Notes.ipynb', 'repo_name': 'chazingtheinfinite/panacea', 'size': '19039'}