<jupyter_start><jupyter_text>Week3-InteractiveVisualizationThis notebook loads a multi-site dataset of subcortical nucleus (thalamus, globus pallidus, and striatum) volumes, performs site harmonization on these volumes using ComBat, then displays a number of interactive figures to explore the harmonized data.<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import neuroCombat as nc<jupyter_output><empty_output><jupyter_text>Load and Prepare Data<jupyter_code>data = pd.read_csv("data_trimmed.csv", index_col=0)
# Map site IDs to integers, as required by neuroCombat()
sites = data['Site'].unique()
site_dict = dict(zip(sites, range(1, len(sites)+1)))
data['Site_no'] = data['Site'].map(site_dict)
#data['Site'].map(site_dict)<jupyter_output><empty_output><jupyter_text>ComBat Site Harmonization<jupyter_code># The list of feature names
features = ['L_str_vol', 'L_GP_vol', 'L_thal_vol', 'R_str_vol', 'R_GP_vol', 'R_thal_vol']

# Perform harmonization with ComBat
harmonized_features = nc.neuroCombat(data[features].transpose(), covars=data[['Site_no', 'Age', 'DX']], batch_col='Site_no', categorical_cols=[], continuous_cols=[]).transpose()
harmonized_features = pd.DataFrame(harmonized_features)
harmonized_features.columns = features

# Add back site data and measures that were not harmonized.
harmonized_features['Site'] = data['Site']
harmonized_features['Age'] = data['Age']
harmonized_features['TBV'] = data['TBV']
harmonized_features['DX'] = data['DX']       # I sure hope that ComBat doesn't change the order of the rows!!!<jupyter_output>[neuroCombat] Creating design matrix
[neuroCombat] Standardizing data across features
[neuroCombat] Fitting L/S model and finding priors
[neuroCombat] Finding parametric adjustments
[neuroCombat] Final adjustment of data<jupyter_text>Compute linear regression models on the harmonized data.Regress nucleus volume against diagnosis and age.<jupyter_code># Linear regression on ComBat harmonized data
from statsmodels.formula.api import ols

models = pd.Series(dtype='object')
model_params = pd.Series(dtype = 'object')
for nucleus in harmonized_features.columns[:-4]:        # Don't include last (covariate) columns
    models[nucleus] = ols(formula = nucleus + ' ~ DX + Age', data = harmonized_features).fit()
    print(models[nucleus].summary())
    model_params[nucleus] = models[nucleus].params
    print(model_params[nucleus])<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:              L_str_vol   R-squared:                       0.001
Model:                            OLS   Adj. R-squared:                 -0.005
Method:                 Least Squares   F-statistic:                    0.1438
Date:                Thu, 11 Jun 2020   Prob (F-statistic):              0.866
Time:                        12:43:49   Log-Likelihood:                -3024.7
No. Observations:                 359   AIC:                             6055.
Df Residuals:                     356   BIC:                             6067.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
=================================================================================
                    coef    std err          t   [...]<jupyter_text>Compute Linear Regression on Unharmonized Data<jupyter_code># Linear regression on unharmonized data
models_unh = pd.Series(dtype='object')
model_unh_params = pd.Series(dtype = 'object')
for nucleus in features:
    models_unh[nucleus] = ols(formula = nucleus + ' ~ DX + Age', data = data).fit()
    print(models_unh[nucleus].summary())
    model_unh_params[nucleus] = models_unh[nucleus].params
    print(model_unh_params[nucleus])<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:              L_str_vol   R-squared:                       0.003
Model:                            OLS   Adj. R-squared:                 -0.003
Method:                 Least Squares   F-statistic:                    0.4976
Date:                Thu, 11 Jun 2020   Prob (F-statistic):              0.608
Time:                        12:43:49   Log-Likelihood:                -3034.0
No. Observations:                 359   AIC:                             6074.
Df Residuals:                     356   BIC:                             6086.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
=================================================================================
                    coef    std err          t   [...]<jupyter_text>Interactive Visualizations<jupyter_code># Set up ipywidgets for interactive chart.
!jupyter nbextension enable --py widgetsnbextension
import ipywidgets as widgets
from ipywidgets import interact, fixed<jupyter_output>Enabling notebook extension jupyter-js-widgets/extension...
      - Validating: [32mOK[0m<jupyter_text>Scatterplots of nucleus volume against Age and Total Brain Volume Split violin plot showing original and harmonized distributions for each site<jupyter_code># Function to plot original and harmonized distributions together, by site.
def plot_dual_distributions(x, y, data1, data2, scale = 'count', ax=None):
    # Plots two distibutions (e.g. harmonized and unharmonized data) as the two halves of stacked violin plots.
    plt.figure()
    vp_data = pd.DataFrame({x: data1[x], y: data1[y], 'Legend': 'Unharmonized'})
    vp_data = vp_data.append(pd.DataFrame({x: data2[x], y: data2[y], 'Legend': 'Harmonized'}))
    ax_h_uh = sns.violinplot(y = y, x = x, data = vp_data, hue = 'Legend', split = True, scale = scale, ax = ax)
    #ax_h_uh.legend().remove()
    plt.show()

# Plot a panel of figures useful for quick data exploration and verification.
def plot_panel(nucleus, data1, data2):
    fig, axes = plt.subplots(ncols=2)
    plt.subplots_adjust(wspace = 0.5)
    sns.scatterplot('Age', nucleus, data = data, hue = "DX", ax=axes[0])
    sns.scatterplot('TBV', nucleus, data = data, hue = "DX", ax=axes[1])
    plot_dual_distributions(x=nucleus, y = 'Site', data1 = data1, data2 = data2, scale = 'count')

interact(plot_panel, nucleus = features, data1 = fixed(data), data2=fixed(harmonized_features))<jupyter_output><empty_output><jupyter_text>Plot fitted values against residuals for the ComBat-harmonized data<jupyter_code># Function to plot fitted values against residuals
def plot_fitted_vs_residual(models, ex_data, x_fitted_col, y_resid_col, display_col, hue_col, title_add = None):
    plt.figure()
    fitted_y = models[display_col].fittedvalues
    resid = models[display_col].resid
    # The following is not ideal. Need to revisit this solution.
    fit_resid_df = pd.DataFrame({'Residual': resid, 'Fitted_Value': fitted_y, 'Site': ex_data['Site'], 'DX': ex_data['DX']})
    ax = sns.scatterplot(fitted_y, resid, hue = hue_col, data=fit_resid_df)
    ax.set_title('Fitted vs. Residual Plot for ' + display_col + title_add)
    ax.set(xlabel = 'Fitted', ylabel = 'Residual')
    ax.legend().remove()
    plt.show()    
# Plot fitted values vs. residuals of ComBat harmonized data

interact(plot_fitted_vs_residual, 
         models = fixed(models), 
         ex_data = fixed(harmonized_features),
         x_fitted_col = fixed('Fitted_Value'), 
         y_resid_col = fixed('Residual'), 
         display_col = features, 
         hue_col = ['DX', 'Site'],
         title_add = fixed(' ComBat Harmonized'))
# Plot fitted values vs. residuals of unharmonized data

interact(plot_fitted_vs_residual, 
         models = fixed(models_unh), 
         ex_data = fixed(data),
         x_fitted_col = fixed('Fitted_Value'), 
         y_resid_col = fixed('Residual'), 
         display_col = features, 
         hue_col = ['DX', 'Site'],
         title_add = fixed(' Unharmonized'))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1378715008015785, 'fraction_numerical': 0.015414971019854482, 'mean_word_length': 2.9599609375, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '144374', 'license': 'CC0-1.0', 'n_tokens_mistral': 2619, 'n_tokens_neox': 2375, 'n_words': 616, 'path': 'Week_3_Deliverable_Data_Visualization/Testing.ipynb', 'repo_name': 'brainhack-school2020/dnmacdon_multisite_smri', 'size': '8109'}