# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

from helpers.utilities import *
# %run helpers/notebook_setup.ipynb

# + tags=["parameters"]
raw_clinical_path = 'data/raw/PatientClinicalData.xlsx'
glossary_path = 'data/clean/clinical/glossary.csv'
samples_path = 'data/clean/samples_list.csv'
# output
clean_clinical_path = 'data/clean/clinical/data.csv'
# -

# Clinical data did not require a separate extraction step, thus I will clean the spreadsheet and explore the data in a single notebook (this one).

# # Data Extraction

clinical_spreadsheets = read_excel(raw_clinical_path, sheet_name=None)

clinical_spreadsheets.keys()

# ## Glossary

clinical_glossary = clinical_spreadsheets['Parameter Glossary']

clinical_glossary.head()

len(clinical_glossary)

# Perfect! Nice and computer readable format. Not much to see there.
#
# The first raw was taken as a column name, I will mitigate it reloading the data with proper headers:

clinical_glossary.to_csv(glossary_path, index=False)
clinical_glossary = read_csv(glossary_path, names=['variable', 'description'])

clinical_glossary.head(n=2)

# And save the clean copy:

clinical_glossary.to_csv(glossary_path)

# ## Clinical data

clinical_data = clinical_spreadsheets['Clinical MetaData']

clinical_data

from helpers.presentation import show_list, compare_sets

# ### Do we have description for each column?

compare_sets(clinical_data.columns, clinical_glossary.variable)

# Okay, so there is a typo the glossary, I will fix that so to make data merges easier in the future:

clinical_glossary.variable = clinical_glossary.variable.replace({'OnTBTReat': 'OnTBTreat'})
clinical_glossary.to_csv(glossary_path, index=False)

# Should work now:

assert set(clinical_data.columns) == set(clinical_glossary.variable)

# # Data exploration and reformatting

# I do not see an atomic column with condition; for simplicity I will create it now:

clinical_data['condition'] = clinical_data.PatientID.str.split('.').str[1]

clinical_data['condition'].value_counts()

# ## Getting to know the variables

# There are only ~~37~~ 39 variables so I will explore each of these, one-by-one.
#
# The descriptions from glossary will come in handy:

glossary = clinical_glossary.set_index('variable').description

# ### 0. PatientID

# We should have data for all patients:

samples_list = read_csv(samples_path)

compare_sets(samples_list.sample_id, clinical_data.PatientID)

# Great!

len(samples_list.sample_id)

# ### 1. The admission date

clinical_data.AdmissionDate.hist();

assert not clinical_data.AdmissionDate.isnull().any()

# ### 2. The date of birth

clinical_data.Birthday.hist();

clinical_data.Birthday.describe()

# The value ranges are plausible (as the has data from ~2015, the oldest participant had about 80 years), no suspicious duplicates.

# ### 3. History of TB

# PrevTB is expected to be a binary column:

set(clinical_data.PrevTB)

clinical_data[~clinical_data.PrevTB.isin(['N', 'Y'])]

# I assume that the 'Unknown' and missing values convey the same information of the history of TB not being known.

from numpy import nan


def yes_no_unknown_to_binary(column):
    mapping = {'N': False, 'Y': True, 'Unknown': nan}
    assert not set(column) - {nan, *mapping.keys()}
    return column.replace(mapping)


clinical_data.PrevTB = yes_no_unknown_to_binary(clinical_data.PrevTB)

counts = clinical_data.PrevTB.fillna('Unknown').value_counts()
counts.to_frame()

counts.plot(kind='bar');

patients_with_history_of_tb = clinical_data[clinical_data.PrevTB == True].PatientID

# ### 4. Form of previous TB

glossary['PrevTBForm']

clinical_data['PrevTBForm'].value_counts().to_frame()

# This column does not satisfy the first normal form - I will add binary columns (as suggested in the description):

clinical_data['was_previous_tb_pulmonary'] = clinical_data.apply(lambda patient:
    # ignore patients with no history of TB (set the value to NaN):
    nan if not (patient.PrevTB == True) else
    'Pulmonary' in str(patient.PrevTBForm)
, axis=1)

# Dr Rachel has confirmed that pleural does count as extrapulmonary:

clinical_data['was_previous_tb_extrapulmonary'] = clinical_data.apply(lambda patient:
    nan if not (patient.PrevTB == True) else
    any(
        term in str(patient.PrevTBForm)
        for term in ['TBM', 'Meningeal', 'Pleural']
    )
, axis=1)

# I will not be creating a separate column for *meningeal* now (it would be strongly correlated with *extrapulmonary*); I do not think that the difference of the single additional observation (one *pleural* patient more the *extrapulmonary* cohort) would make a big impact.

clinical_data[['PrevTB', 'PrevTBForm', 'was_previous_tb_pulmonary', 'was_previous_tb_extrapulmonary']].head(n=8)

# #### Sanity check

# If patient has a form of previous TB given, they should be marked as having history of TB:

patients_with_previous_tb_form_data = clinical_data[~clinical_data.PrevTBForm.isnull()].PatientID
assert not set(patients_with_previous_tb_form_data) - set(patients_with_history_of_tb)

# And it is indeed the case. How many patients have known history but no subtype given?

compare_sets(patients_with_history_of_tb, patients_with_previous_tb_form_data)

# Just one - that's very good!

# ### 5. Treatment for the previous TB

glossary['PrevTBTreat']

clinical_data['PrevTBTreat'].value_counts()

# I will create an additional, binary column for the previous treatment, with a different encoding for missing values:

clinical_data['previous_tb_treatment'] = clinical_data.apply(lambda patient:
    nan if not (patient.PrevTB == True) else
    (
        nan if patient.PrevTBTreat == 'Unknown' else
        patient.PrevTBTreat == 'Y'
    )
, axis=1)

subset = clinical_data[['PatientID', 'PrevTB', 'PrevTBTreat', 'previous_tb_treatment']]
subset.head(n=2)

subset.tail(40).head(8)

# An example of the advantage of the new column is that `148.HC` (58th row) has `NaN` assigned because this patient has no history of TB. The original, `PrevTBTreat` column has 'N' value for this patient. This may be slightly suspicious, given that some other patients (e.g. `149.TMD`/59th row) who also have no history of TB have `NaN` in the original column.

# #### Sanity check

patients_with_previous_tb_treatment = clinical_data[clinical_data.PrevTBTreat == 'Y'].PatientID
assert not set(patients_with_previous_tb_treatment) - set(patients_with_history_of_tb)

# ### 6. Grade of TBM

glossary['TBMGrade']

clinical_data['TBMGrade'].value_counts()

# Already in numeric format, that's great.
#
# **Note:** only two patients with grade three!

# ### 7. HIV Result

glossary['HIVResult']

clinical_data['HIVResult'].value_counts()

# ### 8. CD4

glossary['CD4']

# Distribution is skewed towards zero:

clinical_data['CD4'].hist();

clinical_data['CD4'].describe()

# And there are data for 78/95 patients.

# It seems that that CD4 count is related to HIV/AIDS. A [US governmental site](https://www.hiv.va.gov/patient/diagnosis/labs-CD4-count.asp) states that 200 is a threshold for AIDS diagnosis but they do not give units; I am going to find a more accurate and reliable source and also read [PMC3729334](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3729334/) to understand the context a better at some time.

# ### 9. Is patient on antiretroviral drug(s)?

glossary['ARV']

clinical_data['ARV'].value_counts()

# No need to reformat to another binary format as there are no "Unknown" values.

# Data 79 patients, similar to the 78 above.
# Is there a pattern in the missing data - e.g. are these only missing in healthy controls?

show_list(clinical_data[clinical_data['ARV'].isnull()].PatientID)

# No, it is not the case.

# ### 10. History of treatment with antiretroviral drugs?

glossary['UnARV']

clinical_data['UnARV'].value_counts()

# Interpretation: answering "does the patient has a history of treatment with antiretroviral drugs"?
#  - "ARV-naive" - no, no history of treatment with antiretroviral drugs
#  - "PreviousHistory" - yes

# Very little data - for only 41 patients. Which ones?

patients_with_avr_data = clinical_data[~clinical_data['UnARV'].isnull()]

patients_with_avr_data.condition.value_counts()

patients_with_avr_data.HIVResult.value_counts()

# HIV positive ones.

# **Does lack of data mean "No" or "Unknown"?**

# ### 11. Headache

clinical_data['Headache'].value_counts()

# Reformatting to binary format with NaN for unknown:

clinical_data['Headache'] = yes_no_unknown_to_binary(clinical_data['Headache'])

# ### 12. Duration of headache

glossary['HeadacheD']

clinical_data['HeadacheD'].hist();

clinical_data['HeadacheD'].describe()

# Is 600 days a real figure or a typo? I mean, if some patient reported having a headache for two years it would likely by closer to 730. And pressing an additional 0 when entering "60" is not so unlikely..

# Dr Rachel agrees that it might be a typo, but it is to be confirmed. Creating a column with corrected data:

clinical_data['headache_days_corrected'] = clinical_data['HeadacheD'].replace({600: 60})
clinical_data['headache_days_corrected'].max()

# Othwerise it looks ok, numeric format, no negative values.

# Are headache durations only present for patients with headache?

headache_patients = clinical_data[clinical_data.Headache == True].PatientID

compare_sets(headache_patients, clinical_data[clinical_data.HeadacheD > 0].PatientID)

# That is ok - apparently there are no data for this patient when it comes to headache duration. 

# ### 13. Lethargy

glossary['Lethargy']

clinical_data['Lethargy'].value_counts()

clinical_data['Lethargy'] = yes_no_unknown_to_binary(clinical_data['Lethargy'])

# ### 14. Duration of lethargy

glossary['LethargyD']

clinical_data['LethargyD'].hist();

# ### 15. Vomiting

glossary['Vomiting']

clinical_data['Vomiting'].value_counts()

clinical_data['Vomiting'] = yes_no_unknown_to_binary(clinical_data['Vomiting'])

# ### 16. Duration of vomiting

glossary['VomitingD']

clinical_data['VomitingD'].hist();

compare_sets(
    clinical_data[clinical_data.Vomiting == True].PatientID,
    clinical_data[clinical_data.VomitingD > 0].PatientID
)

# Ok!

# ### 17. Conscious

glossary['Conscious']

# **Note:** the column name may be misleading!

clinical_data['Conscious'].value_counts()

clinical_data['Conscious'] = yes_no_unknown_to_binary(clinical_data['Conscious'])

# ### 18. Conscious

glossary['ConsciousD']

clinical_data['ConsciousD'].hist();

compare_sets(
    clinical_data[clinical_data.Conscious == True].PatientID,
    clinical_data[clinical_data.ConsciousD > 0].PatientID
)

# Ok!

# ### 19. Treatment for TB

glossary['OnTBTreat']

clinical_data['OnTBTreat'].value_counts()

clinical_data['OnTBTreat'] = yes_no_unknown_to_binary(clinical_data['OnTBTreat'])

# Only patients with TB should be on the treatment:

clinical_data[clinical_data['OnTBTreat'] == True].condition.value_counts()

# This is surprising... Maybe these are patients who had a history of TB?

clinical_data[clinical_data['OnTBTreat'] == True].PrevTB.value_counts()

# Yes! This makes sense (the were on treatment at the time of admission).
#
# Are all the noes for patients with history of TB as well?

clinical_data[clinical_data['OnTBTreat'] == False].PrevTB.value_counts()

# Yes. What about current conditions?

clinical_data[clinical_data['OnTBTreat'] == False].condition.value_counts()

# **Note**: this is interesting that all the patients who are still treated for TB are assigned to the CM group. But such patients form the biggest group among all patients with history of TB, thus it might have been expected.

# ### 20. Start date of previous TB treatment

glossary['DateTBTreat']

clinical_data['DateTBTreat'].hist();

clinical_data['DateTBTreat'].describe()

# #### Interpeting the glossary

compare_sets(
    clinical_data[clinical_data.OnTBTreat == True].PatientID,
    clinical_data[~clinical_data.DateTBTreat.isnull()].PatientID
)

# Data inconsistency?

clinical_data[clinical_data.PatientID == '037.HC'][['condition', 'PrevTB', 'PrevTBTreat', 'OnTBTreat', 'DateTBTreat']]

# Maybe. But it may be a matter of description (glossary) interpretation - if the glossary entry was originally placed after the "PrevTBTreat", it would read:
#
# - 'Did the patient receive TB treatment in the past?'
# - 'If yes, when did it start?'
#
# which would make sense.

# ### 21. End date of previous TB treatment

glossary['DateTBTreatStop']

clinical_data['DateTBTreatStop'].hist();

# Because the earliest date is older than of the start date:
#  - we may have data censoring on both ends - such a wonderful data to experiment with survival models,
#  - can we be sure that the column names are not swapped?

# Observation - no overlap with start date (expected given the glossary description):

clinical_data[
    (~clinical_data.DateTBTreat.isnull())
    &
    (~clinical_data.DateTBTreatStop.isnull())
].empty

# #### Minor inconsistencies with treatment status?

compare_sets(
    clinical_data[clinical_data.OnTBTreat == False].PatientID,
    clinical_data[~clinical_data.DateTBTreatStop.isnull()].PatientID
)

# The one odd example:

# Ok!

columns_of_interest = [
    'PrevTB', 'PrevTBTreat', 'OnTBTreat', 'DateTBTreat', 'DateTBTreatStop',
    'AdmissionDate'
]
clinical_data[clinical_data.PatientID == '224.CM'][columns_of_interest]

# If the treatment stopped in 2014-09, why was the patient still on treatment when admitted?
# One explanation would be that the admission was before this date. But it was not.

# Dr Rachel notes that there is an explanation of such discrepancies: the patient may have stopped treatment early; then at the admission the doctor may have placed them back (immediately) on the treatment; this is thought not to be uncommon.

# Because it is just a single patient, I think that overall this column is ok.

# ### 22. Date of lumbar puncture

glossary['DateLP']

clinical_data['DateLP'].hist(bins=20, figsize=(10, 4));

# I expect that there are dates for all patients as all of them had the procedure:

assert not clinical_data['DateLP'].isnull().any()

(clinical_data['DateLP'] >= clinical_data['AdmissionDate']).value_counts()

time_to_puncture = clinical_data['DateLP'] - clinical_data['AdmissionDate']
time_to_puncture.value_counts()

outlier_row_id = time_to_puncture.idxmax()
clinical_data.loc[outlier_row_id][['PatientID', 'DateLP','AdmissionDate', 'DateCD4', 'DateAntiTB']]

# This is suspicious!

# **Update:** this record can be corrected - the electronic record (input date) hints us that the admission date for this patient had a typo.

new_admission_date = clinical_data.loc[outlier_row_id, 'AdmissionDate'].replace(year=2015)
new_admission_date

clinical_data.loc[outlier_row_id, 'AdmissionDate'] = new_admission_date

# ### 23. Appearance of CSF #1

glossary['App1']

clinical_data['App1'].value_counts()

# Perfect!

# ### 24. Appearance of CSF #2

glossary['App2']

clinical_data['App2'].value_counts()

# There are different terms than in the first one...

clinical_data[clinical_data.App1 == 'Mild cloudy/ opalescent'][['App1', 'App2']]

# And the values do not correlate. Interesting, though I do not understand this fully:
#  - were the two observations performed at different time points, or
#  - were the two looking for different things?

# **Update:** based on a comment from Dr Rachel, the observations were collected at the same time, but looking at different characteristics: clear/cloudy and color/colorless respectively.

# ### 25. Red blood cells count

clinical_data['RCC'].hist();

clinical_data['RCC'].describe()

# How was RCC measured? In CSF? That would seem likely given a lot of zeros!
#
# If this is the case, there should be a few healthy controls with RCC > 0:

clinical_data[clinical_data.RCC > 0].condition.value_counts()

# ### 26. White blood cells count

clinical_data['WCC'].hist();

clinical_data[clinical_data.WCC > 0].condition.value_counts()

# ### 27. Percent of neutrophils among WBC

clinical_data['%Neutro'].hist();

# The range is 0-100 - ok.

# The spike at zero could be due to many samples having zero WCC. Is it?

zero_percent = clinical_data[clinical_data['%Neutro'] == 0]
zero_percent[zero_percent.WCC == 0].empty

# No, it is not.

# Sanity check: %Neutro > 0 only if WCC > 0?

(clinical_data[clinical_data['%Neutro'] > 0].WCC > 0).all()

# ### 28. Percent of lymphocytes among WBC

clinical_data['%Lympho'].hist();

# Range is ok. %Lympho > 0 only if WCC > 0?

(clinical_data[clinical_data['%Lympho'] > 0].WCC > 0).all()

# Great!

# ### 29. Protein in CSF

glossary['Protein']

clinical_data['Protein'].hist();

# Protein >= 0 → ok.
# No other sanity checks if could think of.

# Only one missing value:

clinical_data['Protein'].isnull().sum()

# ### 30. Glucose in CSF

glossary['CSFGlucose']

clinical_data['CSFGlucose'].hist();

clinical_data['CSFGlucose'].isnull().sum()

# ### 31. Was anti-TB treatment started (this admission)?

clinical_data['AntiTB'].value_counts()

clinical_data[clinical_data['AntiTB'] == 'Y'].condition.value_counts()

# Given the Monday discussion I understand that the treatment may be decided upon before the final diagnosis is made, thus I just assume that the classes are ok.

clinical_data['AntiTB'] = yes_no_unknown_to_binary(clinical_data['AntiTB'])

# ### 32. The start date of anti-TB treatment (this admission)

clinical_data['DateAntiTB'].hist(bins=20, figsize=(10, 4));

# All who have a start date should also be marked as having the treatment:

assert clinical_data[~clinical_data.DateAntiTB.isnull()].AntiTB.all()

# ### 33. Was steroid started (this admission)?

glossary['SteroidsStarted']

clinical_data['SteroidsStarted'].value_counts()

clinical_data['SteroidsStarted'] = yes_no_unknown_to_binary(clinical_data['SteroidsStarted'])

# ### 34. The date of starting steroid (this admission)

clinical_data['SteroidDate'].hist(bins=20, figsize=(10, 4));

# All who have a start date should also be marked as having the treatment:

assert clinical_data[~clinical_data.SteroidDate.isnull()].SteroidsStarted.all()

# ### 35. Type of steroid

glossary['SteroidType']

clinical_data['SteroidType'].value_counts()

# Data for all the 34 - ok!

# ### 36. Death

glossary['Death']

clinical_data['Death'].value_counts()

clinical_data['Death'] = clinical_data['Death'].replace({'Alive': False, 'Death': True})

# It might be good to know when was the end of the study to know the right censoring horizon.

# ### 37. Date of death

clinical_data['DateDeath'].hist(bins=20, figsize=(10, 4));

assert clinical_data[~clinical_data.DateDeath.isnull()].Death.all()

# One could envision many sanity checks for this:

for date_column in ['SteroidDate', 'DateAntiTB', 'DateLP', 'AdmissionDate']: 
    assert not (clinical_data[date_column] > clinical_data['DateDeath']).any()

# ### 38. CD4 test date

glossary['DateCD4']

clinical_data.DateCD4.hist(bins=20, figsize=(10, 4));

(clinical_data.DateCD4 >= clinical_data.AdmissionDate).sum()

(clinical_data.DateCD4 < clinical_data.AdmissionDate).sum()

# A quick note: about a third of the CD4 tests was performed before the current admission; 
# it is plausible as the test might have been done immediately after the HIV diagnosis or around the time when someone was suspected to have HIV.
#
# The glossary entry reaffirms that this is ok.

compare_sets(
    clinical_data[~clinical_data.CD4.isnull()].PatientID,
    clinical_data[~clinical_data.DateCD4.isnull()].PatientID
)

# ### 39. Gender

clinical_data['Sex'].value_counts()

# Perfect!

# ## Variable exploration summary

# There are many variables which are by definition dependent one on another. These are:

# this dictionary was organised in shuch a way that the column with
# the most data is the key (on the left side, preferring the re-formatted ones)
relations = {
    'PrevTB': [
        'PrevTBForm', 'PrevTBTreat',
        # to be confirmed
        'OnTBTreat', 'DateTBTreat', 'DateTBTreatStop',
        # addded columns
        'was_previous_tb_pulmonary', 'was_previous_tb_extrapulmonary',
        'previous_tb_treatment'
    ],
    'HIVResult': ['CD4', 'ARV', 'UnARV'],
    'Headache': ['HeadacheD', 'headache_days_corrected'],
    'Lethargy': ['LethargyD'],
    'Vomiting': ['VomitingD'],
    'Conscious': ['ConsciousD'],
    # OnTBTreat not only defines the presence of the dates, but also which date was recorded!
    'OnTBTreat': ['DateTBTreat', 'DateTBTreatStop'],
    # These dates are mutually exclusive:
    'DateTBTreat': ['DateTBTreatStop'],
    'WCC': ['%Neutro', '%Lympho'],
    'AntiTB': ['DateAntiTB'],
    'SteroidsStarted': ['SteroidDate', 'SteroidType'],
    'Death': ['DateDeath'],
    'CD4': ['DateCD4']
}

# ## Data completeness

# ### Patient-wise

# Are there any patients for whom there is too little data in general?

# Exclude the dependent variables so that we don't have unequal weights of missing values just because there were some variables measured only in a subset of patients:

key_variables = set(clinical_data.columns) - {e for l in relations.values() for e in l}
len(key_variables)

subset = clinical_data[key_variables].set_index('PatientID')
missing_data = (100 * subset.isnull().sum(axis=1) / len(subset.columns))
missing_data = missing_data[missing_data > 0]

missing_data.to_frame()

# Up to 20% of missing "key" variables per patient (key variable defined as not being dependent - due to the data collection procedure - on other variables).
#
# It should not be a big problem.
#
# We may want to exclude the patient with 20% of missing data (also, see the note on PLS tolerance to missing values below). I will keep all the patients for now, because I may be able to fill the missing values for some of the variables and it may lower the percentage to more acceptable level.

# ### Variable-wise

missing_percent = (100 * subset.isnull().sum(axis=0) / len(subset))
missing = missing_percent[missing_percent != 0]
missing.sort_index().sort_values().to_frame()

# Up to 3.2% of missing data on the selected variables - good enough!

# ### Overall

ratio_missing = subset.isnull().sum(axis=1).sum() / (len(subset) * len(subset.columns))
f'{100 * ratio_missing:.2f}%'

# (<NAME> et. al, 2001) states that PLS can handle up to 10-20% of missing data for dataset with 20 observations / 20 variables if the values are not missing in systematic pattern (thus the exclusion of non-key columns) and that the tolerance to missing values increases with the increase in dataset size.
#
# Thus there should be no problem with 1% of missing data.
#
# Why is 20/20 anyhow relevant to this dataset?
# - there are 22 clinical variables; even though these will be changed (admission date or date of birth does not really count, but I will get age/ survival out of these), there will be ~20 major variables in the end
# - the subgroups (by condition) if analyzed separately have about 20 observations

clinical_data.set_index('PatientID').to_csv(clean_clinical_path)

StarcoderdataJupyterScriptsDedupFiltered