tutorial: Meta-Analysis and real life data using Python

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math

# Set path
# Ubuntu
work_dir = "/home/julien/Documents/Projects/datasets/meta_analysis"
# Windows
#work_dir = "C:\\Users\hernandez\Dropbox (Personal)\Python\meta_analysis"
os.chdir(work_dir)
# Verify the path
cwd = os.getcwd() 
# Print the current directory 
print("Current working directory is:", cwd) 

files = os.listdir(cwd)
df = pd.read_excel(io = files[files.index("Effect_size_pooled_c-fos_v2.xlsx")], 
                   sheet_name ='c-fos_only',
                   header     = 1,
                   index_col  = None)

# lower case to avoid later conflicts
df.columns = [c.lower() for c in df.columns]
# display columns
df.shape

# Create dictionnary of the species encoded
specie_dict ={0 : 'rat', 1 : 'mouse'}
# Map the binary coding onto the dictionnary
df['specie_str'] = df['specie '].map(specie_dict)

# Select data of interest
DOI = df[['study id','specie ', 
          'specie_str', 
          'structure / tissue', 
          'structure_renamed',              
          'standardized es',    
          'standardized es se', 
          'weight factor', 
          'n']]

DOI.head(10)

(DOI
.groupby(['specie_str', 'structure_renamed'])
.count()
.sort_values('n', ascending=False)
)

(DOI
.groupby(['specie_str'])
.count()['study id']
)

# This function will add a few computations to the dataset.
# Computations are self explanatory, and equations are in methods
def compute_meta_data(dataf):
    return (dataf
            .assign(vi_SE2           = lambda d:  1/d['weight factor']) # Equation 15
            .assign(standardized_ESU = lambda d:  d['standardized es'] 
                    * (1-(3/(4*(d['n'])-9)))) # Equation 9
            .assign(wES              = lambda d:  d['standardized es'] 
                    * d['weight factor'])
            .assign(wES2             = lambda d:  d['weight factor']   
                    * d['standardized es'] * d['standardized es']) 
            .assign(w2               = lambda d:  d['weight factor']   
                    * d['weight factor']))

# This function will compute, for each specie&structure, the random effect parameters
def compute_random_effects(dataf):
    return (dataf
            .groupby(['specie_str','structure_renamed'])
            .apply(lambda d: pd.Series({  
            "sum_wi":          d['weight factor'].sum(),
            "sum_wi2":         d['w2'].sum(),
            "k":               d['w2'].count(), #num studies
            "wxES":            d['wES'].sum(),
            "wxES2":           d['wES2'].sum()     
            })))

# Compute the Q and v0 parameters (see methods)
def compute_random_variance(dataf):
    return (dataf
            .assign(Q  = lambda d:  abs(d['wxES2']-(d['wxES']**2/d['sum_wi']))) 
            # For struture with only 1 effect size, the denominator of v0 will 
            # automatically be equal to 0. hence v0 will not be computed.
            .assign(v0 = lambda d:  (d['Q']-(d['k']-1))/(d['sum_wi']-
                                                        (d['sum_wi2']/d['sum_wi'])))
           )

# Turn v0 to 0 if negative (missing reference!)
def zero_if_negative(dataf):
        num = dataf['v0']
        num[num < 0] = 0
        dataf['v0'] = num
        return dataf

# Compute corrected parameters (see methods)
def apply_corrections(dataf):
    return (dataf
            .assign(v0_plus_vi       = lambda d:  d['vi_SE2']+d['v0'])  
            .assign(wi_corr          = lambda d:  1/d['v0_plus_vi'])  
            .assign(wxES_corr        = lambda d:  d['wi_corr']*d['standardized_ESU'])  
            .assign(wxESsq_corr      = lambda d:  d['wi_corr']*(d['standardized_ESU']**2))
            .assign(SE_corr          = lambda d:  (1/(d['wi_corr'])**0.5))  
            .fillna(0) # NaN are replaced by 0
            .replace(np.inf, 0)
)

def calculate_constants(dataf):
    return (dataf
            .groupby(['specie_str', 'structure_renamed'])
            .apply(lambda d: pd.Series({  
                        "sum_wxES_corr":        d['wxES_corr'].sum(),
                        "sum_wxESsq_corr":      d['wxESsq_corr'].sum(), #num studies                  
                        "sum_wi_corr":          abs(d['wi_corr']).sum(), 
             })) 
             .fillna(0) # NaN are replaced by 0
             .replace(np.inf, 0)
           )

def I2(Q,dfg):
#     if pd.isnull((Q < dfg)):
#         I2 = 0
#         return I2
#     else:
        I2=((Q-dfg)/Q)*100
        return I2

def calculate_mean_se(dataf):
    return (dataf
            .assign(ES_mean       = lambda d: d['sum_wxES_corr']/d['sum_wi_corr']) 
            .assign(SE_mean       = lambda d: (1/d['sum_wi_corr'])**0.5)
            .assign(z             = lambda d: d['ES_mean']/d['SE_mean'])
            .assign(high_CI       = lambda d: d['ES_mean']+(1.96*d['SE_mean']))
            .assign(low_CI        = lambda d: d['ES_mean']-(1.96*d['SE_mean']))                                              
           # .assign(k_val         = lambda d: d['w2'].count())#num studies  
            .assign(I_val         = lambda d: I2(d['Q'],d['k']-1))
           )

def calculate_number_es(dataf):
    return (dataf
            .groupby(['specie_str', 'structure_renamed'])
            .apply(lambda d: pd.Series({  
                        "k_val": d['w2'].count(),
             })) 
           )

# Compute & add meta variables to dataset
cfos_df = DOI.pipe(compute_meta_data)
# Compute & add random effect parameters
rand_df = (cfos_df
.pipe(compute_random_effects)
.pipe(compute_random_variance)
.pipe(zero_if_negative)
)
# Merge datasets
cfos_ana = pd.merge(cfos_df, rand_df, on=['structure_renamed','specie_str'])
cfos_ana = apply_corrections(cfos_ana)
# Compute constants for each structure
constants = calculate_constants(cfos_ana)
cfos_ana = pd.merge(cfos_ana, constants, on=['structure_renamed','specie_str'])
# Compute mean & se for each structure
cfos_ana = (cfos_ana
.pipe(calculate_mean_se))
# Calculate number of ES per structure
nES = (cfos_ana
.pipe(calculate_number_es)           
)
cfos_ana = pd.merge(cfos_ana, nES, on=['structure_renamed','specie_str'])

# Sample dataframe
cfos_ana.head(3)

# By convention, Q = 0 if Q < k-1, so that the precision of a random effects summary estimate will not exceed the precision of a fixed effect summary estimate
cfos_ana.loc[cfos_ana['Q'] < (cfos_ana['k_val']-1), 'I_val'] = 0 

# Prepare data for plotting
to_plot = (cfos_ana
#.groupby(['specie_str'])
.drop_duplicates(subset =["specie_str", "structure_renamed"], inplace = False)
.dropna()
#.groupby(['specie_str', 'structure_renamed'])
#.head(20)
)
to_plot = to_plot[['specie_str', 
                   'structure_renamed',
                   'ES_mean', 
                   'low_CI',
                   'high_CI', 
                   'z', 
                   'I_val', 
                   'k_val']]
# Visualize higher effect sizes per species
(to_plot
.groupby(["specie_str"])
.apply(lambda x: x.sort_values(["ES_mean"], ascending = False))
.reset_index(drop=True)
)

#import plotnine as p9
import plotnine as p9
from plotnine import ggplot, geom_point, aes, theme, element_text, save_as_pdf_pages
p9.options.figure_size = (6.4, 4.8)

# Plot the effect sizes for each structure, separately for each specie
(p9.ggplot() +
p9.geom_point(data = cfos_ana, 
              mapping = p9.aes('structure_renamed', 
                               'k_val', 
                               color='ES_mean', 
                               size='ES_mean'
                              )
             ) +
p9.facet_wrap("specie_str") + 
theme(axis_text_x=element_text(rotation=90, hjust=1))
)