|
| 1 | +import numpy as np |
| 2 | + |
| 3 | +def vcorrcoef(X, y): # return a correlation between each row of X with y |
| 4 | + Xm = np.reshape(np.mean(X, axis=1), (X.shape[0], 1)) |
| 5 | + ym = np.mean(y) |
| 6 | + r_num = np.sum((X-Xm)*(y-ym), axis=1) |
| 7 | + r_den = np.sqrt(np.sum((X-Xm)**2, axis=1)*np.sum((y-ym)**2)) |
| 8 | + r = r_num/r_den |
| 9 | + return r |
| 10 | + |
| 11 | + |
| 12 | +def elementwise_corrcoef(X, Y): |
| 13 | + # X and Y are each of shape num_observations X num_element |
| 14 | + # computes the correlation between each element of X and Y |
| 15 | + Xm = X.mean(axis=0) |
| 16 | + Ym = Y.mean(axis=0) |
| 17 | + r_num = np.sum((X-Xm)*(Y-Ym), axis=0) |
| 18 | + |
| 19 | + r_den = np.sqrt(np.sum((X-Xm)**2, axis=0)*np.sum((Y-Ym)**2, axis=0)) |
| 20 | + r = r_num/r_den |
| 21 | + return r |
| 22 | + |
| 23 | + |
| 24 | +def elementwise_spearman(X,Y): |
| 25 | + order = X.argsort(axis=0) |
| 26 | + X_ranks = order.argsort(axis=0) |
| 27 | + order = Y.argsort(axis=0) |
| 28 | + Y_ranks = order.argsort(axis=0) |
| 29 | + return elementwise_corrcoef(X_ranks, Y_ranks) |
| 30 | + |
| 31 | + |
| 32 | +def dice_coefficient(mask1,mask2): |
| 33 | + dice = np.sum(mask1*mask2)*2.0 / (np.sum(mask1) + np.sum(mask2)) |
| 34 | + return dice |
| 35 | + |
| 36 | + |
| 37 | +''' |
| 38 | +LINEAR REGRESSION --- CLOSED-FORM SOLUTION |
| 39 | +''' |
| 40 | + |
| 41 | + |
| 42 | +def closed_form(X, Y, intercept=False): # functions that computes the Least Squares Estimates |
| 43 | + if intercept: |
| 44 | + X = np.concatenate((X, np.ones([X.shape[0], 1])), axis=1) |
| 45 | + return np.linalg.inv(X.transpose().dot(X)).dot(X.transpose()).dot(Y) |
| 46 | + |
| 47 | + |
| 48 | +def mse(X, Y, w): # function that computes the Mean Square Error (MSE) |
| 49 | + return np.mean((Y-np.matmul(X, w))**2) |
| 50 | + |
| 51 | + |
| 52 | + |
| 53 | +''' |
| 54 | +
|
| 55 | +def closed_form_3d(X,Y): |
| 56 | + return np.matmul(np.matmul(np.linalg.inv(np.matmul(X.transpose(0,2,1),X)),X.transpose(0,2,1)),Y) |
| 57 | +
|
| 58 | +def lme_stats_3d(X,Y): |
| 59 | + #add an intercept |
| 60 | + X=np.concatenate((X,np.ones((X.shape[0],X.shape[1],1))),axis=2) |
| 61 | + [num_comparisons,num_observations,num_predictors] = X.shape |
| 62 | + [num_comparisons,num_observations,num_features] = Y.shape |
| 63 | +
|
| 64 | + w=closed_form_3d(X,Y) |
| 65 | +
|
| 66 | + residuals = Y-np.matmul(X, w) |
| 67 | + MSE = (((residuals)**2).sum(axis=1)/(num_observations-num_predictors)) |
| 68 | +
|
| 69 | +
|
| 70 | + var_b = np.expand_dims(MSE, axis=1)*np.expand_dims(np.linalg.inv(np.matmul(X.transpose(0,2,1),X)).diagonal(axis1=1,axis2=2), axis=2) |
| 71 | + sd_b = np.sqrt(var_b) # standard error on the Betas |
| 72 | + ts_b = w/sd_b # calculate t-values for the Betas |
| 73 | + p_values =[2*(1-stats.t.cdf(np.abs(ts_b[:,i,:]),(num_observations-num_predictors))) for i in range(ts_b.shape[1])] # calculate a p-value map for each predictor |
| 74 | +
|
| 75 | + return ts_b,p_values,w,residuals |
| 76 | +
|
| 77 | +non_nan_idx = (np.isnan(voxelwise_array).sum(axis=(0,1))==0) |
| 78 | +
|
| 79 | +# take out the voxels which have null values |
| 80 | +X=voxelwise_array[:,:6,non_nan_idx].transpose(2,0,1) |
| 81 | +Y=voxelwise_array[:,6:,non_nan_idx].transpose(2,0,1) |
| 82 | +
|
| 83 | +
|
| 84 | +ts_b,p_values,w,residuals = lme_stats_3d(X,Y) |
| 85 | +
|
| 86 | +x_name=['Somatomotor','Dorsal Comp','DMN', 'Prior Modeling 1', 'Prior Modeling 2', 'Prior Modeling 3'] |
| 87 | +y_name=['group','temporal_std','VE_spatial','GS_corr','DVARS_corr','FD_corr'] |
| 88 | +
|
| 89 | +fig,axes = plt.subplots(nrows=len(x_name), ncols=len(y_name),figsize=(12*len(y_name),3*len(x_name))) |
| 90 | +
|
| 91 | +
|
| 92 | +for i,x_label in zip(range(len(x_name)),x_name): |
| 93 | + for j,y_label in zip(range(len(y_name)),y_name): |
| 94 | + ax=axes[i,j] |
| 95 | +
|
| 96 | + stat_map=np.zeros(voxelwise_array.shape[2]) |
| 97 | + stat_map[non_nan_idx]=ts_b[:,j,i] |
| 98 | +
|
| 99 | + ax.set_title('T-value of {} on {}'.format(y_label,x_label), fontsize=15) |
| 100 | + plot_stat_map(analysis_functions.recover_3D(mask_file,stat_map),bg_img='DSURQE.nii.gz', axes=ax, cut_coords=(0,1,2,3,4,5), display_mode='y') |
| 101 | +''' |
0 commit comments