<s> """<STR_LIT>""" <EOL> import logging <EOL> import numpy as np <EOL> from scipy . stats import norm <EOL> from neurosynth . base import imageutils <EOL> from neurosynth . analysis import stats <EOL> from os . path import join , exists <EOL> from os import makedirs <EOL> logger = logging . getLogger ( '<STR_LIT>' ) <EOL> def analyze_features ( dataset , features = None , image_type = '<STR_LIT>' , <EOL> threshold = <NUM_LIT> , q = <NUM_LIT> , output_dir = None , <EOL> prefix = None ) : <EOL> """<STR_LIT>""" <EOL> if features is None : <EOL> features = dataset . get_feature_names ( ) <EOL> if output_dir is None : <EOL> result = np . zeros ( ( dataset . masker . n_vox_in_mask , len ( features ) ) ) <EOL> for i , f in enumerate ( features ) : <EOL> ids = dataset . get_studies ( features = f , frequency_threshold = threshold ) <EOL> ma = MetaAnalysis ( dataset , ids , q = q ) <EOL> if output_dir is None : <EOL> result [ : , i ] = ma . images [ image_type ] <EOL> else : <EOL> pfx = f if prefix is None else prefix + '<STR_LIT:_>' + f <EOL> ma . save_results ( output_dir = output_dir , prefix = pfx ) <EOL> if output_dir is None : <EOL> return result <EOL> class MetaAnalysis ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , dataset , ids , ids2 = None , q = <NUM_LIT> , prior = <NUM_LIT:0.5> , <EOL> min_studies = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> self . dataset = dataset <EOL> mt = dataset . image_table <EOL> self . selected_ids = list ( set ( mt . ids ) & set ( ids ) ) <EOL> self . selected_id_indices = np . in1d ( mt . ids , ids ) <EOL> unselected_id_indices = ~ self . selected_id_indices if ids2 == None else np . in1d ( mt . ids , ids2 ) <EOL> logger . debug ( "<STR_LIT>" ) <EOL> n_selected = len ( self . selected_ids ) <EOL> n_unselected = np . sum ( unselected_id_indices ) <EOL> n_mappables = n_selected + n_unselected <EOL> n_selected_active_voxels = mt . data . dot ( self . selected_id_indices ) <EOL> n_unselected_active_voxels = mt . data . dot ( unselected_id_indices ) <EOL> pF = ( n_selected * <NUM_LIT:1.0> ) / n_mappables <EOL> pA = np . array ( ( mt . data . sum ( axis = <NUM_LIT:1> ) * <NUM_LIT:1.0> ) / n_mappables ) . squeeze ( ) <EOL> logger . debug ( "<STR_LIT>" ) <EOL> pAgF = n_selected_active_voxels * <NUM_LIT:1.0> / n_selected <EOL> pAgU = n_unselected_active_voxels * <NUM_LIT:1.0> / n_unselected <EOL> pFgA = pAgF * pF / pA <EOL> logger . debug ( "<STR_LIT>" ) <EOL> pAgF_prior = prior * pAgF + ( <NUM_LIT:1> - prior ) * pAgU <EOL> pFgA_prior = pAgF * prior / pAgF_prior <EOL> def p_to_z ( p , sign ) : <EOL> p = p / <NUM_LIT:2> <EOL> p [ p < <NUM_LIT> ] = <NUM_LIT> <EOL> z = np . abs ( norm . ppf ( p ) ) * sign <EOL> z [ np . isinf ( z ) ] = norm . ppf ( <NUM_LIT> ) * - <NUM_LIT:1> <EOL> return z <EOL> p_vals = stats . one_way ( <EOL> np . squeeze ( n_selected_active_voxels ) , n_selected ) <EOL> p_vals [ p_vals < <NUM_LIT> ] = <NUM_LIT> <EOL> z_sign = np . sign ( <EOL> n_selected_active_voxels - np . mean ( <EOL> n_selected_active_voxels ) ) . ravel ( ) <EOL> pAgF_z = p_to_z ( p_vals , z_sign ) <EOL> fdr_thresh = stats . fdr ( p_vals , q ) <EOL> pAgF_z_FDR = imageutils . threshold_img ( <EOL> pAgF_z , fdr_thresh , p_vals , mask_out = '<STR_LIT>' ) <EOL> cells = np . squeeze ( <EOL> np . array ( [ [ n_selected_active_voxels , n_unselected_active_voxels ] , <EOL> [ n_selected - n_selected_active_voxels , n_unselected - <EOL> n_unselected_active_voxels ] ] ) . T ) <EOL> p_vals = stats . two_way ( cells ) <EOL> p_vals [ p_vals < <NUM_LIT> ] = <NUM_LIT> <EOL> z_sign = np . sign ( pAgF - pAgU ) . ravel ( ) <EOL> pFgA_z = p_to_z ( p_vals , z_sign ) <EOL> fdr_thresh = stats . fdr ( p_vals , q ) <EOL> pFgA_z_FDR = imageutils . threshold_img ( <EOL> pFgA_z , fdr_thresh , p_vals , mask_out = '<STR_LIT>' ) <EOL> self . images = { <EOL> '<STR_LIT>' : pA , <EOL> '<STR_LIT>' : pAgF , <EOL> '<STR_LIT>' : pFgA , <EOL> ( '<STR_LIT>' % prior ) : pAgF_prior , <EOL> ( '<STR_LIT>' % prior ) : pFgA_prior , <EOL> '<STR_LIT>' : pAgF_z , <EOL> '<STR_LIT>' : pFgA_z , <EOL> ( '<STR_LIT>' % q ) : pAgF_z_FDR , <EOL> ( '<STR_LIT>' % q ) : pFgA_z_FDR <EOL> } <EOL> if min_studies > <NUM_LIT:0> : <EOL> if isinstance ( min_studies , int ) : <EOL> min_studies = float ( <EOL> min_studies ) / n_mappables <EOL> vox_to_exclude = np . where ( pA < min_studies ) [ <NUM_LIT:0> ] <EOL> for k in self . images : <EOL> self . images [ k ] [ vox_to_exclude ] = <NUM_LIT:0> <EOL> def save_results ( self , output_dir = '<STR_LIT:.>' , prefix = '<STR_LIT>' , prefix_sep = '<STR_LIT:_>' , <EOL> image_list = None ) : <EOL> """<STR_LIT>""" <EOL> if prefix == '<STR_LIT>' : <EOL> prefix_sep = '<STR_LIT>' <EOL> if not exists ( output_dir ) : <EOL> makedirs ( output_dir ) <EOL> logger . debug ( "<STR_LIT>" ) <EOL> if image_list is None : <EOL> image_list = self . images . keys ( ) <EOL> for suffix , img in self . images . items ( ) : <EOL> if suffix in image_list : <EOL> filename = prefix + prefix_sep + suffix + '<STR_LIT>' <EOL> outpath = join ( output_dir , filename ) <EOL> imageutils . save_img ( img , outpath , self . dataset . masker ) </s>
