diff --git a/docs/outputs.rst b/docs/outputs.rst index 94a4a149a..0ebafd1dc 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -133,19 +133,19 @@ TEDICA codes Code Classification Description ===== =============== ======================================================== I001 rejected Manual exclusion -I002 rejected Rho greater than Kappa or more significant voxels - in S0 model than R2 model -I003 rejected S0 Dice is higher than R2 Dice and high variance +I002 rejected Rho greater than Kappa +I003 rejected More significant voxels in S0 model than R2 model +I004 rejected S0 Dice is higher than R2 Dice and high variance explained -I004 rejected Noise F-value is higher than signal F-value and high +I005 rejected Noise F-value is higher than signal F-value and high variance explained -I005 ignored No good components found -I006 rejected Mid-Kappa component -I007 ignored Low variance explained -I008 rejected Artifact candidate type A -I009 rejected Artifact candidate type B -I010 ignored ign_add0 -I011 ignored ign_add1 +I006 ignored No good components found +I007 rejected Mid-Kappa component +I008 ignored Low variance explained +I009 rejected Mid-Kappa artifact type A +I010 rejected Mid-Kappa artifact type B +I011 ignored ign_add0 +I012 ignored ign_add1 ===== =============== ======================================================== Visual reports diff --git a/tedana/decomposition/eigendecomp.py b/tedana/decomposition/eigendecomp.py index e8b2180c9..68ba1c9ec 100644 --- a/tedana/decomposition/eigendecomp.py +++ b/tedana/decomposition/eigendecomp.py @@ -60,13 +60,13 @@ def run_mlepca(data): return u, s, v -def kundu_tedpca(ct_df, n_echos, kdaw, rdaw, stabilize=False): +def kundu_tedpca(comptable, n_echos, kdaw, rdaw, stabilize=False): """ Select PCA components using Kundu's decision tree approach. Parameters ---------- - ct_df : :obj:`pandas.DataFrame` + comptable : :obj:`pandas.DataFrame` Component table with relevant metrics: kappa, rho, and normalized variance explained. n_echos : :obj:`int` @@ -83,85 +83,85 @@ def kundu_tedpca(ct_df, n_echos, kdaw, rdaw, stabilize=False): Returns ------- - ct_df : :obj:`pandas.DataFrame` + comptable : :obj:`pandas.DataFrame` Component table with components classified as 'accepted', 'rejected', or 'ignored'. """ - eigenvalue_elbow = getelbow(ct_df['normalized variance explained'], + eigenvalue_elbow = getelbow(comptable['normalized variance explained'], return_val=True) - diff_varex_norm = np.abs(np.diff(ct_df['normalized variance explained'])) + diff_varex_norm = np.abs(np.diff(comptable['normalized variance explained'])) lower_diff_varex_norm = diff_varex_norm[(len(diff_varex_norm) // 2):] varex_norm_thr = np.mean([lower_diff_varex_norm.max(), diff_varex_norm.min()]) - varex_norm_min = ct_df['normalized variance explained'][ + varex_norm_min = comptable['normalized variance explained'][ (len(diff_varex_norm) // 2) + np.arange(len(lower_diff_varex_norm))[lower_diff_varex_norm >= varex_norm_thr][0] + 1] - varex_norm_cum = np.cumsum(ct_df['normalized variance explained']) + varex_norm_cum = np.cumsum(comptable['normalized variance explained']) fmin, fmid, fmax = utils.getfbounds(n_echos) if int(kdaw) == -1: - lim_idx = utils.andb([ct_df['kappa'] < fmid, - ct_df['kappa'] > fmin]) == 2 - kappa_lim = ct_df.loc[lim_idx, 'kappa'].values + lim_idx = utils.andb([comptable['kappa'] < fmid, + comptable['kappa'] > fmin]) == 2 + kappa_lim = comptable.loc[lim_idx, 'kappa'].values kappa_thr = kappa_lim[getelbow(kappa_lim)] - lim_idx = utils.andb([ct_df['rho'] < fmid, ct_df['rho'] > fmin]) == 2 - rho_lim = ct_df.loc[lim_idx, 'rho'].values + lim_idx = utils.andb([comptable['rho'] < fmid, comptable['rho'] > fmin]) == 2 + rho_lim = comptable.loc[lim_idx, 'rho'].values rho_thr = rho_lim[getelbow(rho_lim)] stabilize = True LGR.info('kdaw set to -1. Switching TEDPCA method to ' 'kundu-stabilize') elif int(rdaw) == -1: - lim_idx = utils.andb([ct_df['rho'] < fmid, ct_df['rho'] > fmin]) == 2 - rho_lim = ct_df.loc[lim_idx, 'rho'].values + lim_idx = utils.andb([comptable['rho'] < fmid, comptable['rho'] > fmin]) == 2 + rho_lim = comptable.loc[lim_idx, 'rho'].values rho_thr = rho_lim[getelbow(rho_lim)] else: kappa_thr = np.average( - sorted([fmin, (getelbow(ct_df['kappa'], return_val=True) / 2), fmid]), + sorted([fmin, (getelbow(comptable['kappa'], return_val=True) / 2), fmid]), weights=[kdaw, 1, 1]) rho_thr = np.average( - sorted([fmin, (getelbow_cons(ct_df['rho'], return_val=True) / 2), fmid]), + sorted([fmin, (getelbow_cons(comptable['rho'], return_val=True) / 2), fmid]), weights=[rdaw, 1, 1]) # Reject if low Kappa, Rho, and variance explained - is_lowk = ct_df['kappa'] <= kappa_thr - is_lowr = ct_df['rho'] <= rho_thr - is_lowe = ct_df['normalized variance explained'] <= eigenvalue_elbow + is_lowk = comptable['kappa'] <= kappa_thr + is_lowr = comptable['rho'] <= rho_thr + is_lowe = comptable['normalized variance explained'] <= eigenvalue_elbow is_lowkre = is_lowk & is_lowr & is_lowe - ct_df.loc[is_lowkre, 'classification'] = 'rejected' - ct_df.loc[is_lowkre, 'rationale'] += 'P001;' + comptable.loc[is_lowkre, 'classification'] = 'rejected' + comptable.loc[is_lowkre, 'rationale'] += 'P001;' # Reject if low variance explained - is_lows = ct_df['normalized variance explained'] <= varex_norm_min - ct_df.loc[is_lows, 'classification'] = 'rejected' - ct_df.loc[is_lows, 'rationale'] += 'P002;' + is_lows = comptable['normalized variance explained'] <= varex_norm_min + comptable.loc[is_lows, 'classification'] = 'rejected' + comptable.loc[is_lows, 'rationale'] += 'P002;' # Reject if Kappa over limit - is_fmax1 = ct_df['kappa'] == F_MAX - ct_df.loc[is_fmax1, 'classification'] = 'rejected' - ct_df.loc[is_fmax1, 'rationale'] += 'P003;' + is_fmax1 = comptable['kappa'] == F_MAX + comptable.loc[is_fmax1, 'classification'] = 'rejected' + comptable.loc[is_fmax1, 'rationale'] += 'P003;' # Reject if Rho over limit - is_fmax2 = ct_df['rho'] == F_MAX - ct_df.loc[is_fmax2, 'classification'] = 'rejected' - ct_df.loc[is_fmax2, 'rationale'] += 'P004;' + is_fmax2 = comptable['rho'] == F_MAX + comptable.loc[is_fmax2, 'classification'] = 'rejected' + comptable.loc[is_fmax2, 'rationale'] += 'P004;' if stabilize: temp7 = varex_norm_cum >= 0.95 - ct_df.loc[temp7, 'classification'] = 'rejected' - ct_df.loc[temp7, 'rationale'] += 'P005;' - under_fmin1 = ct_df['kappa'] <= fmin - ct_df.loc[under_fmin1, 'classification'] = 'rejected' - ct_df.loc[under_fmin1, 'rationale'] += 'P006;' - under_fmin2 = ct_df['rho'] <= fmin - ct_df.loc[under_fmin2, 'classification'] = 'rejected' - ct_df.loc[under_fmin2, 'rationale'] += 'P007;' - - n_components = ct_df.loc[ct_df['classification'] == 'accepted'].shape[0] + comptable.loc[temp7, 'classification'] = 'rejected' + comptable.loc[temp7, 'rationale'] += 'P005;' + under_fmin1 = comptable['kappa'] <= fmin + comptable.loc[under_fmin1, 'classification'] = 'rejected' + comptable.loc[under_fmin1, 'rationale'] += 'P006;' + under_fmin2 = comptable['rho'] <= fmin + comptable.loc[under_fmin2, 'classification'] = 'rejected' + comptable.loc[under_fmin2, 'rationale'] += 'P007;' + + n_components = comptable.loc[comptable['classification'] == 'accepted'].shape[0] LGR.info('Selected {0} components with Kappa threshold: {1:.02f}, Rho ' 'threshold: {2:.02f}'.format(n_components, kappa_thr, rho_thr)) - return ct_df + return comptable def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, @@ -319,18 +319,18 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, # Normalize each component's time series vTmixN = stats.zscore(comp_ts, axis=1).T LGR.info('Making initial component selection guess from PCA results') - _, ct_df, betasv, v_T = model.fitmodels_direct( + _, comptable, betasv, v_T = model.fitmodels_direct( catd, comp_ts.T, eimum, t2s, t2sG, tes, combmode, ref_img, mmixN=vTmixN, full_sel=False, label='mepca_', verbose=verbose) # varex_norm overrides normalized varex computed by fitmodels_direct - ct_df['normalized variance explained'] = varex_norm + comptable['normalized variance explained'] = varex_norm pcastate = {'method': method, 'voxel_comp_weights': voxel_comp_weights, 'varex': varex, 'comp_ts': comp_ts, - 'comptable': ct_df} + 'comptable': comptable} # Save state LGR.info('Saving PCA results to: {}'.format(fname)) @@ -344,7 +344,7 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, voxel_comp_weights = pcastate['voxel_comp_weights'] varex = pcastate['varex'] comp_ts = pcastate['comp_ts'] - ct_df = pcastate['comptable'] + comptable = pcastate['comptable'] np.savetxt('mepca_mix.1D', comp_ts.T) @@ -357,22 +357,23 @@ def tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, io.filewrite(comp_maps, 'mepca_OC_components.nii', ref_img) # Add new columns to comptable for classification - ct_df['classification'] = 'accepted' - ct_df['rationale'] = '' + comptable['classification'] = 'accepted' + comptable['rationale'] = '' # Select components using decision tree if method == 'kundu': - ct_df = kundu_tedpca(ct_df, n_echos, kdaw, rdaw, stabilize=False) + comptable = kundu_tedpca(comptable, n_echos, kdaw, rdaw, stabilize=False) elif method == 'kundu-stabilize': - ct_df = kundu_tedpca(ct_df, n_echos, kdaw, rdaw, stabilize=True) + comptable = kundu_tedpca(comptable, n_echos, kdaw, rdaw, stabilize=True) elif method == 'mle': LGR.info('Selected {0} components with MLE dimensionality ' - 'detection'.format(ct_df.shape[0])) + 'detection'.format(comptable.shape[0])) - ct_df.to_csv('comp_table_pca.txt', sep='\t', index=True, - index_label='component', float_format='%.6f') + comptable['rationale'] = comptable['rationale'].str.rstrip(';') + comptable.to_csv('comp_table_pca.txt', sep='\t', index=True, + index_label='component', float_format='%.6f') - sel_idx = ct_df['classification'] == 'accepted' + sel_idx = comptable['classification'] == 'accepted' n_components = np.sum(sel_idx) voxel_kept_comp_weighted = (voxel_comp_weights[:, sel_idx] * varex[None, sel_idx]) diff --git a/tedana/selection/select_comps.py b/tedana/selection/select_comps.py index 5db15f78c..727bf46cd 100644 --- a/tedana/selection/select_comps.py +++ b/tedana/selection/select_comps.py @@ -83,6 +83,7 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): midk = [] ign = [] all_comps = np.arange(comptable.shape[0]) + # acc remains a full list that is whittled down over criteria acc = np.arange(comptable.shape[0]) # If user has specified @@ -95,17 +96,22 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): # Move decision columns to end comptable = comptable[[c for c in comptable if c not in cols_at_end] + [c for c in cols_at_end if c in comptable]] + comptable['rationale'] = comptable['rationale'].str.rstrip(';') return comptable """ - Do some tallies for no. of significant voxels + Tally number of significant voxels for cluster-extent thresholded R2 and S0 + model F-statistic maps. """ - countnoise = np.zeros(n_comps) comptable['countsigFR2'] = F_R2_clmaps.sum(axis=0) comptable['countsigFS0'] = F_S0_clmaps.sum(axis=0) """ - Make table of dice values + Generate Dice values for R2 and S0 models + - dice_FR2: Dice value of cluster-extent thresholded maps of R2-model betas + and F-statistics. + - dice_FS0: Dice value of cluster-extent thresholded maps of S0-model betas + and F-statistics. """ comptable['dice_FR2'] = np.zeros(all_comps.shape[0]) comptable['dice_FS0'] = np.zeros(all_comps.shape[0]) @@ -119,16 +125,24 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): comptable.loc[np.isnan(comptable['dice_FS0']), 'dice_FS0'] = 0 """ - Make table of noise gain + Generate three metrics of component noise: + - countnoise: Number of "noise" voxels (voxels highly weighted for + component, but not from clusters) + - signal-noise_t: T-statistic for two-sample t-test of F-statistics from + "signal" voxels (voxels in clusters) against "noise" voxels (voxels not + in clusters) for R2 model. + - signal-noise_p: P-value from t-test. """ comptable['countnoise'] = 0 comptable['signal-noise_t'] = 0 comptable['signal-noise_p'] = 0 for i_comp in all_comps: + # index voxels significantly loading on component but not from clusters comp_noise_sel = ((np.abs(Z_maps[:, i_comp]) > 1.95) & (Z_clmaps[:, i_comp] == 0)) comptable.loc[i_comp, 'countnoise'] = np.array( comp_noise_sel, dtype=np.int).sum() + # NOTE: Why only compare distributions of *unique* F-statistics? noise_FR2_Z = np.log10(np.unique(F_R2_maps[comp_noise_sel, i_comp])) signal_FR2_Z = np.log10(np.unique( F_R2_maps[Z_clmaps[:, i_comp] == 1, i_comp])) @@ -140,13 +154,23 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): comptable.loc[np.isnan(comptable['signal-noise_p']), 'signal-noise_p'] = 0 """ - Assemble decision table + Assemble decision table with five metrics: + - Kappa values ranked from largest to smallest + - R2-model F-score map/beta map Dice scores ranked from largest to smallest + - Signal F > Noise F t-statistics ranked from largest to smallest + - Number of "noise" voxels (voxels highly weighted for component, but not + from clusters) ranked from smallest to largest + - Number of voxels with significant R2-model F-scores within clusters + ranked from largest to smallest + + Smaller values (i.e., higher ranks) across metrics indicate more BOLD + dependence and less noise. """ d_table_rank = np.vstack([ n_comps - stats.rankdata(comptable['kappa'], method='ordinal'), n_comps - stats.rankdata(comptable['dice_FR2'], method='ordinal'), n_comps - stats.rankdata(comptable['signal-noise_t'], method='ordinal'), - stats.rankdata(countnoise, method='ordinal'), + stats.rankdata(comptable['countnoise'], method='ordinal'), n_comps - stats.rankdata(comptable['countsigFR2'], method='ordinal')]).T n_decision_metrics = d_table_rank.shape[1] comptable['d_table_score'] = d_table_rank.sum(axis=1) @@ -155,24 +179,35 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): Step 1: Reject anything that's obviously an artifact a. Estimate a null variance """ - temp_rej0 = all_comps[(comptable['rho'] > comptable['kappa']) | - ((comptable['countsigFS0'] > comptable['countsigFR2']) & - (comptable['countsigFR2'] > 0))] - comptable.loc[temp_rej0, 'classification'] = 'rejected' - comptable.loc[temp_rej0, 'rationale'] += 'I002;' - + # Rho is higher than Kappa + temp_rej0a = all_comps[(comptable['rho'] > comptable['kappa'])] + comptable.loc[temp_rej0a, 'classification'] = 'rejected' + comptable.loc[temp_rej0a, 'rationale'] += 'I002;' + + # Number of significant voxels for S0 model is higher than number for R2 + # model *and* number for R2 model is greater than zero. + temp_rej0b = all_comps[((comptable['countsigFS0'] > comptable['countsigFR2']) & + (comptable['countsigFR2'] > 0))] + comptable.loc[temp_rej0b, 'classification'] = 'rejected' + comptable.loc[temp_rej0b, 'rationale'] += 'I003;' + rej = np.union1d(temp_rej0a, temp_rej0b) + + # Dice score for S0 maps is higher than Dice score for R2 maps and variance + # explained is higher than the median across components. temp_rej1 = all_comps[(comptable['dice_FS0'] > comptable['dice_FR2']) & (comptable['variance explained'] > np.median(comptable['variance explained']))] comptable.loc[temp_rej1, 'classification'] = 'rejected' - comptable.loc[temp_rej1, 'rationale'] += 'I003;' - rej = np.union1d(temp_rej0, temp_rej1) + comptable.loc[temp_rej1, 'rationale'] += 'I004;' + rej = np.union1d(temp_rej1, rej) + # T-value is less than zero (noise has higher F-statistics than signal in + # map) and variance explained is higher than the median across components. temp_rej2 = acc[(comptable.loc[acc, 'signal-noise_t'] < 0) & (comptable.loc[acc, 'variance explained'] > np.median(comptable['variance explained']))] comptable.loc[temp_rej2, 'classification'] = 'rejected' - comptable.loc[temp_rej2, 'rationale'] += 'I004;' + comptable.loc[temp_rej2, 'rationale'] += 'I005;' rej = np.union1d(temp_rej2, rej) acc = np.setdiff1d(acc, rej) @@ -188,6 +223,8 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): f. Estimate a low and high variance """ # Step 2a + # Upper limit for variance explained is median across components with high + # Kappa values. High Kappa is defined as Kappa above Kappa elbow. varex_upper_p = np.median( comptable.loc[comptable['kappa'] > getelbow(comptable['kappa'], return_val=True), 'variance explained']) @@ -195,80 +232,90 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): # NOTE: We're not sure why this is done, nor why it's specifically done # three times. Need to look into this deeper, esp. to make sure the 3 # isn't a hard-coded reference to the number of echoes. + # Reduce components to investigate as "good" to ones in which change in + # variance explained is less than the limit defined above.... What? for nn in range(3): ncls = comptable.loc[ncls].loc[ comptable.loc[ ncls, 'variance explained'].diff() < varex_upper_p].index.values - # Compute elbows - kappas_lim = comptable.loc[comptable['kappa'] < utils.getfbounds(n_echos)[-1], 'kappa'] - kappa_elbow = np.min((getelbow(kappas_lim, return_val=True), + # Compute elbows from other elbows + kappas_under_f01 = (comptable.loc[comptable['kappa'] < + utils.getfbounds(n_echos)[-1], 'kappa']) + # NOTE: Would an elbow from all Kappa values *ever* be lower than one from + # a subset of lower values? + kappa_elbow = np.min((getelbow(kappas_under_f01, return_val=True), getelbow(comptable['kappa'], return_val=True))) rho_elbow = np.mean((getelbow(comptable.loc[ncls, 'rho'], return_val=True), getelbow(comptable['rho'], return_val=True), utils.getfbounds(n_echos)[0])) - # Initial guess of good components based on Kappa and Rho elbows - good_guess = ncls[(comptable.loc[ncls, 'kappa'] >= kappa_elbow) & - (comptable.loc[ncls, 'rho'] < rho_elbow)] + # Provisionally accept components based on Kappa and Rho elbows + acc_prov = ncls[(comptable.loc[ncls, 'kappa'] >= kappa_elbow) & + (comptable.loc[ncls, 'rho'] < rho_elbow)] - if len(good_guess) == 0: + if len(acc_prov) == 0: LGR.warning('No BOLD-like components detected') ign = sorted(np.setdiff1d(all_comps, rej)) comptable.loc[ign, 'classification'] = 'ignored' - comptable.loc[ign, 'rationale'] += 'I005;' + comptable.loc[ign, 'rationale'] += 'I006;' # Move decision columns to end comptable = comptable[[c for c in comptable if c not in cols_at_end] + [c for c in cols_at_end if c in comptable]] + comptable['rationale'] = comptable['rationale'].str.rstrip(';') return comptable - kappa_rate = ((np.max(comptable.loc[good_guess, 'kappa']) - - np.min(comptable.loc[good_guess, 'kappa'])) / - (np.max(comptable.loc[good_guess, 'variance explained']) - - np.min(comptable.loc[good_guess, 'variance explained']))) + # Calculate "rate" for kappa: kappa range divided by variance explained + # range, for potentially accepted components + # NOTE: What is the logic behind this? + kappa_rate = ((np.max(comptable.loc[acc_prov, 'kappa']) - + np.min(comptable.loc[acc_prov, 'kappa'])) / + (np.max(comptable.loc[acc_prov, 'variance explained']) - + np.min(comptable.loc[acc_prov, 'variance explained']))) kappa_ratios = kappa_rate * comptable['variance explained'] / comptable['kappa'] varex_lower = stats.scoreatpercentile( - comptable.loc[good_guess, 'variance explained'], LOW_PERC) + comptable.loc[acc_prov, 'variance explained'], LOW_PERC) varex_upper = stats.scoreatpercentile( - comptable.loc[good_guess, 'variance explained'], HIGH_PERC) + comptable.loc[acc_prov, 'variance explained'], HIGH_PERC) """ Step 3: Get rid of midk components; i.e., those with higher than max decision score and high variance """ - max_good_d_score = EXTEND_FACTOR * len(good_guess) * n_decision_metrics + max_good_d_score = EXTEND_FACTOR * len(acc_prov) * n_decision_metrics midk = acc[(comptable.loc[acc, 'd_table_score'] > max_good_d_score) & (comptable.loc[acc, 'variance explained'] > EXTEND_FACTOR * varex_upper)] comptable.loc[midk, 'classification'] = 'rejected' - comptable.loc[midk, 'rationale'] += 'I006;' + comptable.loc[midk, 'rationale'] += 'I007;' acc = np.setdiff1d(acc, midk) + acc_prov = np.setdiff1d(acc_prov, midk) """ Step 4: Find components to ignore """ - good_guess = np.setdiff1d(good_guess, midk) - loaded = np.union1d(good_guess, acc[comptable.loc[acc, 'variance explained'] > varex_lower]) - ign = np.setdiff1d(acc, loaded) + high_varex = np.union1d(acc_prov, acc[comptable.loc[acc, 'variance explained'] > varex_lower]) + ign = np.setdiff1d(acc, high_varex) # ignore low variance components ign = np.setdiff1d( ign, ign[comptable.loc[ign, 'd_table_score'] < max_good_d_score]) ign = np.setdiff1d(ign, ign[comptable.loc[ign, 'kappa'] > kappa_elbow]) comptable.loc[ign, 'classification'] = 'ignored' - comptable.loc[ign, 'rationale'] += 'I007;' + comptable.loc[ign, 'rationale'] += 'I008;' acc = np.setdiff1d(acc, ign) """ - Step 5: Scrub the set + Step 5: Scrub the set if there are components that haven't been rejected or + ignored, but are still not listed in the possible accepted group. """ - if len(acc) > len(good_guess): + if len(acc) > len(acc_prov): + comptable['d_table_score_scrub'] = np.nan # Recompute the midk steps on the limited set to clean up the tail d_table_rank = np.vstack([ len(acc) - stats.rankdata(comptable.loc[acc, 'kappa'], method='ordinal'), len(acc) - stats.rankdata(comptable.loc[acc, 'dice_FR2'], method='ordinal'), len(acc) - stats.rankdata(comptable.loc[acc, 'signal-noise_t'], method='ordinal'), - stats.rankdata(countnoise[acc], method='ordinal'), + stats.rankdata(comptable.loc[acc, 'countnoise'], method='ordinal'), len(acc) - stats.rankdata(comptable.loc[acc, 'countsigFR2'], method='ordinal')]).T - comptable['d_table_score_scrub'] = np.nan comptable.loc[acc, 'd_table_score_scrub'] = d_table_rank.sum(1) num_acc_guess = int(np.mean([ np.sum((comptable.loc[acc, 'kappa'] > kappa_elbow) & @@ -284,7 +331,7 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): candartA, candartA[comptable.loc[candartA, 'variance explained'] > varex_upper * EXTEND_FACTOR]) comptable.loc[candartA, 'classification'] = 'rejected' - comptable.loc[candartA, 'rationale'] += 'I008;' + comptable.loc[candartA, 'rationale'] += 'I009;' midk = np.union1d(midk, candartA) # Rejection candidate based on artifact type B: candartB @@ -294,11 +341,11 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): candartB = np.intersect1d( candartB, candartB[comptable.loc[candartB, 'variance explained'] > varex_lower * EXTEND_FACTOR]) - midk = np.union1d(midk, candartB) comptable.loc[candartB, 'classification'] = 'rejected' - comptable.loc[candartB, 'rationale'] += 'I009;' + comptable.loc[candartB, 'rationale'] += 'I010;' + midk = np.union1d(midk, candartB) - # Find comps to ignore + # Find components to ignore new_varex_lower = stats.scoreatpercentile( comptable.loc[acc[:num_acc_guess], 'variance explained'], LOW_PERC) @@ -310,7 +357,7 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): candart) ign_add0 = np.setdiff1d(ign_add0, midk) comptable.loc[ign_add0, 'classification'] = 'ignored' - comptable.loc[ign_add0, 'rationale'] += 'I010;' + comptable.loc[ign_add0, 'rationale'] += 'I011;' ign = np.union1d(ign, ign_add0) ign_add1 = np.intersect1d( @@ -318,11 +365,10 @@ def selcomps(seldict, comptable, mmix, manacc, n_echos): acc[comptable.loc[acc, 'variance explained'] > new_varex_lower]) ign_add1 = np.setdiff1d(ign_add1, midk) comptable.loc[ign_add1, 'classification'] = 'ignored' - comptable.loc[ign_add1, 'rationale'] += 'I011;' - ign = np.union1d(ign, ign_add1) - acc = np.setdiff1d(acc, np.union1d(midk, ign)) + comptable.loc[ign_add1, 'rationale'] += 'I012;' # Move decision columns to end comptable = comptable[[c for c in comptable if c not in cols_at_end] + [c for c in cols_at_end if c in comptable]] + comptable['rationale'] = comptable['rationale'].str.rstrip(';') return comptable