99from scipy .stats import ttest_ind
1010from bct import degrees_und
1111
12+
1213def communicability_bin (adjacency , normalize = False ):
1314 """
1415 Computes the communicability of pairs of nodes in `adjacency`
@@ -101,79 +102,86 @@ def communicability_wei(adjacency):
101102
102103 return cmc
103104
104- def rich_feeder_peripheral (x , sc , stat = 'median' ):
105+
106+ def rich_feeder_peripheral (x , sc , stat = 'median' ):
105107 """
106- Calculates median or mean "connectivity" in rich (hub to hub), feeder (hub to non-hub),
107- and peripheral (non-hub to non-hub) links.
108-
109- This function takes a square symmetric correlation/connectivity matrix `x`
110- and computes the median or mean value within rich, feeder, and peripheral links.
111- Hubs are defined againt a threshold (from 0 to the maximum degree `k`,
112- defined on the structural connectivity matrix, `sc`).
113-
108+ Calculates median or mean "connectivity" in rich (hub to hub), feeder
109+ (hub to non-hub), and peripheral (non-hub to non-hub) links.
110+
111+ This function takes a square symmetric correlation/connectivity matrix `x`
112+ and computes the median or mean value within rich, feeder, and peripheral
113+ links. Hubs are defined againt a threshold (from 0 to the maximum degree
114+ `k`, defined on the structural connectivity matrix, `sc`).
114115
115116 Parameters
116117 ----------
117- x : (N, N) symmetric numpy.ndarray.
118- sc : (N, N) symmetric numpy.ndarray
119- binary structural connectivity matrix
120- stat : 'median' (default) or 'mean'.
121- statistic to use over rich/feeder/peripheral links.
122-
118+ x : (N, N) numpy.ndarray
119+ Symmetric correlation or connectivity matrix
120+ sc : (N, N) numpy.ndarray
121+ binary structural connectivity matrix
122+ stat : str
123+ statistic to use over rich/feeder/peripheral links
124+ 'mean' or 'median' (default)
125+
123126 Returns
124127 -------
125- rfp : (3, k) numpy.ndarray of median rich (0), feeder (1), and peripheral (2)
126- values, defined by `x`. `k` is the maximum degree defined on `sc`.
127- pvals : (3, k) numpy.ndarray p-value for each link, computed using Welch's t-test.
128- Rich links are compared against non-rich links. Feeder links are
129- compared against peripheral links. Peripheral links are compared
130- against feeder links. T-test is two-sided.
131-
128+ rfp : (3, k) numpy.ndarray
129+ Array of median rich (0), feeder (1), and peripheral (2)
130+ values, defined by `x`. `k` is the maximum degree defined on `sc`.
131+ pvals : (3, k) numpy.ndarray
132+ p-value for each link, computed using Welch's t-test.
133+ Rich links are compared against non-rich links. Feeder links are
134+ compared against peripheral links. Peripheral links are compared
135+ against feeder links. T-test is two-sided.
136+
132137 Author
133138 ------
134139 This code was written by Justine Hansen who promises to fix and even
135140 optimize the code should any issues arise, provided you let her know.
136-
137141 """
142+
138143 nnodes = len (sc )
139- mask = np .triu (np .ones (nnodes ),1 ) > 0
144+ mask = np .triu (np .ones (nnodes ), 1 ) > 0
140145 node_degree = degrees_und (sc )
141146 k = np .max (node_degree ).astype (np .int64 )
142- rfp_label = np .zeros ((len (sc [mask ]),k ))
143-
147+ rfp_label = np .zeros ((len (sc [mask ]), k ))
148+
144149 for i in range (k ): # for each degree threshold
145150 hub_idx = np .where (node_degree >= i ) # find the hubs
146- hub = np .zeros ([nnodes ,1 ])
147- hub [hub_idx ,:] = 1
148-
149- rfp = np .zeros ([nnodes ,nnodes ]) # for each link, define rfp
151+ hub = np .zeros ([nnodes , 1 ])
152+ hub [hub_idx , :] = 1
153+
154+ rfp = np .zeros ([nnodes , nnodes ]) # for each link, define rfp
150155 for ii in range (nnodes ):
151156 for iii in range (nnodes ):
152157 if hub [ii ] + hub [iii ] == 2 :
153- rfp [ii ,iii ] = 1 # rich
158+ rfp [ii , iii ] = 1 # rich
154159 if hub [ii ] + hub [iii ] == 1 :
155- rfp [ii ,iii ] = 2 # feeder
160+ rfp [ii , iii ] = 2 # feeder
156161 if hub [ii ] + hub [iii ] == 0 :
157- rfp [ii ,iii ] = 3 # peripheral
158- rfp_label [:,i ] = rfp [mask ]
159-
160- rfp = np .zeros ([3 ,k ])
161- pvals = np .zeros ([3 ,k ])
162+ rfp [ii , iii ] = 3 # peripheral
163+ rfp_label [:, i ] = rfp [mask ]
164+
165+ rfp = np .zeros ([3 , k ])
166+ pvals = np .zeros ([3 , k ])
162167 for i in range (k ):
163168
164169 if stat == 'median' :
165- rfp [0 ,i ] = np .median (x [rfp_label [:,i ] == 1 ]) # rich
166- rfp [1 ,i ] = np .median (x [rfp_label [:,i ] == 2 ]) # feeder
167- rfp [2 ,i ] = np .median (x [rfp_label [:,i ] == 3 ]) # peripheral
168-
170+ rfp [0 , i ] = np .median (x [rfp_label [:, i ] == 1 ]) # rich
171+ rfp [1 , i ] = np .median (x [rfp_label [:, i ] == 2 ]) # feeder
172+ rfp [2 , i ] = np .median (x [rfp_label [:, i ] == 3 ]) # peripheral
173+
169174 if stat == 'mean' :
170- rfp [0 ,i ] = np .mean (x [rfp_label [:,i ] == 1 ]) # rich
171- rfp [1 ,i ] = np .mean (x [rfp_label [:,i ] == 2 ]) # feeder
172- rfp [2 ,i ] = np .mean (x [rfp_label [:,i ] == 3 ]) # peripheral
173-
174- # p-value
175- _ , pvals [0 ,i ] = ttest_ind (x [rfp_label [:,i ] == 1 ], x [rfp_label [:,i ] != 1 ], equal_var = False ) # Welch's t-test
176- _ , pvals [1 ,i ] = ttest_ind (x [rfp_label [:,i ] == 2 ], x [rfp_label [:,i ] == 3 ], equal_var = False )
177- _ , pvals [2 ,i ] = ttest_ind (x [rfp_label [:,i ] == 3 ], x [rfp_label [:,i ] == 2 ], equal_var = False )
178-
179- return rfp , pvals
175+ rfp [0 , i ] = np .mean (x [rfp_label [:, i ] == 1 ]) # rich
176+ rfp [1 , i ] = np .mean (x [rfp_label [:, i ] == 2 ]) # feeder
177+ rfp [2 , i ] = np .mean (x [rfp_label [:, i ] == 3 ]) # peripheral
178+
179+ # p-value (Welch's t-test)
180+ _ , pvals [0 ,i ] = ttest_ind (x [rfp_label [:, i ] == 1 ], \
181+ x [rfp_label [:, i ] != 1 ], equal_var = False )
182+ _ , pvals [1 ,i ] = ttest_ind (x [rfp_label [:, i ] == 2 ], \
183+ x [rfp_label [:, i ] == 3 ], equal_var = False )
184+ _ , pvals [2 ,i ] = ttest_ind (x [rfp_label [:, i ] == 3 ], \
185+ x [rfp_label [:, i ] == 2 ], equal_var = False )
186+
187+ return rfp , pvals
0 commit comments