1 #include "mrilib.h"
2
3 #undef PMAX
4 #define PMAX 0.9999f /*** don't process p-values >= PMAX ***/
5
6 #undef QTOZ /* convert q-value to z-score */
7 #undef ZTOQ /* vice-versa */
8 #define QTOZ(x) qginv(0.5*x)
9 #define ZTOQ(x) (2.0*qg(x))
10
11 #undef ZTOP
12 #undef QBOT
13 #define ZTOP 9.0 /* largest value of z(q) to output */
14 #define QBOT 2.25718e-19 /* smallest value of q to return */
15 #define PBOT 1.e-15 /* smallest value of p to use */
16
17 /*--------------------------------------------------------------------------*/
18 #undef INMASK
19 #define INMASK(vv) ( FDRmask == NULL || FDRmask[vv] != 0 )
20 static byte *FDRmask = NULL ;
21
mri_fdr_setmask(byte * mmm)22 void mri_fdr_setmask( byte *mmm ){ FDRmask = mmm ; } /* 27 Mar 2009 */
23
24 /*--------------------------------------------------------------------------*/
25 #undef QSTHRESH
26 #define QSTHRESH 0.15 /* q threshold for computing FDR_mdfv */
27
28 static floatvec *FDR_mdfv = NULL ; /* missed detection fraction vs. log10(p) */
29
mri_fdr_getmdf(void)30 floatvec *mri_fdr_getmdf(void){ return FDR_mdfv; } /* 22 Oct 2008 */
31
32 /*--------------------------------------------------------------------------*/
33 /* Estimate m1 = number of true positives.
34 Actually, estimates m0 = number of true negatives, then m1 = nq-m0.
35 Method:
36 * Build a histogram of large-ish p-values [0.15..0.95], which should be
37 approximately uniformly distributed.
38 * Then find m0 by estimating the average-ish level of this histogram.
39 * Then multiply by the number of such bins (20) that would cover the
40 entire p-range of 0..1 to get the m0 estimate.
41 * If something bad happens, return value is -1.
42 * Note that this is the NUMBER of true positives, NOT the fraction
43 of true positives (which is how I use the symbol m_1 in presentations)
44 ----------------------------------------------------------------------------*/
45
estimate_m1(int nq,float * qq)46 static int estimate_m1( int nq , float *qq )
47 {
48 int jj , kk , nh=0 , hist[16] , ma,mb, mone ;
49
50 if( nq < 233 || qq == NULL ) return -1 ; /* not enuf data to bother with */
51
52 /* build histogram with 16 bins of length 0.05,
53 to count p-values ranging from 0.15 to 0.95 */
54
55 for( kk=0 ; kk < 16 ; kk++ ) hist[kk] = 0.0f ;
56 for( jj=0 ; jj < nq ; jj++ ){ /* jj-th point is in the kk-th bin */
57 kk = (int)( (qq[jj]-0.15f)*20.0f ) ; if( kk < 0 || kk > 15 ) continue ;
58 hist[kk]++ ; nh++ ;
59 }
60 if( nh < 160 ) return -1 ; /* too few p-values in [0.15..0.95] range!! */
61 qsort_int( 16 , hist ) ; /* sort; use central values to get 'average' */
62 /* level of the uniform part of the histogram */
63
64 /* estimate m1 two different ways, take the smaller value */
65
66 /* m0 from median 4 bins averaged */
67 ma = nq - 20.0f * ( hist[6] + 2*hist[7] + 2*hist[8] + hist[9] ) / 6.0f ;
68
69 /* m0 from median 6 bins averaged */
70 mb = nq - 20.0f * ( hist[5] + 2*hist[6] + 2*hist[7]
71 + 2*hist[8] + 2*hist[9] + hist[10] ) / 10.0f ;
72
73 mone = MIN(ma,mb) ; return mone ;
74 }
75
76 /*--------------------------------------------------------------------------*/
77 /*! Take an image of statistics and convert to FDR-ized z-scores (in place):
78 - im must be in float format
79 - if statcode > 0, the data is a statistic to be converted to
80 a p-value first
81 (input data values == 0 are masked out as having p=1)
82 - otherwise, the data is already p-value-ized
83 (input data values < 0 or >= 1 will be masked out)
84 - if flags&1==1, then the function tries to be compatible with the
85 3dFDR program in '-old' mode:
86 - in 3dFDR -old, processed input values that give p==1 are
87 still counted in the number of thresholdings performed,
88 which will make the output z-scores smaller
89 - if flags&1==0, then this function will NOT count such
90 p==1 voxels at all, which will make the z-scores larger
91 - since this function actually sorts the p-values while 3dFDR
92 just bins them, small differences will be present anyhoo
93 - also, in '-old' mode, the m0 correction will not be made
94 - if flags&2==1, then the q-values are corrected for arbitrary
95 correlation structure -- this is not usually necessary for FMRI
96 task activation data
97 - if flags&4==1, then the output is q-values, not z-values
98 - to mask, set input values to a statistic that will give p==1
99 (e.g., 0.0 for t, F, or rho; 1.0 for p) -- and set flags=0;
100 masked out voxels will be set to 0 (or 1 if flags&4 is set).
101 *//*------------------------------------------------------------------------*/
102
mri_fdrize(MRI_IMAGE * im,int statcode,float * stataux,int flags)103 int mri_fdrize( MRI_IMAGE *im, int statcode, float *stataux, int flags )
104 {
105 float *far ;
106 int ii,jj, nvox, doz, qsmal=0 ;
107 float *qq, nthr, fbad, mone=0.0f,qfac=1.0f ; int *iq, nq ; double qval,qmin ;
108
109 ENTRY("mri_fdrize") ;
110
111 KILL_floatvec(FDR_mdfv) ; FDR_mdfv = NULL ; /* erase the past */
112
113 if( im == NULL || im->kind != MRI_float ) RETURN(0) ;
114 far = MRI_FLOAT_PTR(im); if( far == NULL ) RETURN(0) ;
115 nvox = im->nvox ;
116 doz = (flags&4) == 0 ;
117
118 /*----- convert to p-value first? -----*/
119
120 if( FUNC_IS_STAT(statcode) ){ /* conversion to p-value */
121 STATUS("convert to p-value") ;
122 for( ii=0 ; ii < nvox ; ii++ ){
123 if( far[ii] != 0.0f && INMASK(ii) )
124 far[ii] = THD_stat_to_pval( fabsf(far[ii]), statcode,stataux ) ;
125 else
126 far[ii] = 1.0f ; /* will be ignored */
127 }
128 } else { /* already supposed to be p-value */
129 STATUS("input is p-value") ;
130 for( ii=0 ; ii < nvox ; ii++ ) /* scan for bad values */
131 if( far[ii] < 0.0f || far[ii] > 1.0f || !INMASK(ii) ) far[ii] = 1.0f ;
132 }
133
134 qq = (float *)malloc(sizeof(float)*nvox) ; /* array of p-values */
135 iq = (int *)malloc(sizeof(int )*nvox) ; /* voxel indexes */
136
137 if( qq == NULL || iq == NULL ){
138 ERROR_message("mri_fdrize: out of memory!") ;
139 if( qq != NULL ) free(qq) ;
140 if( iq != NULL ) free(iq) ;
141 RETURN(0) ;
142 }
143
144 /*---- build array of p-values into qq and source voxel indexes into iq ----*/
145
146 STATUS("find reasonable p-values") ;
147 fbad = (doz) ? 0.0f : 1.0f ; /* output value for masked voxels */
148 for( nq=ii=0 ; ii < nvox ; ii++ ){
149 if( far[ii] >= 0.0f && far[ii] < PMAX ){ /* reasonable p-value */
150 qq[nq] = (far[ii] > PBOT) ? far[ii] : PBOT ;
151 iq[nq] = ii ; nq++ ;
152 } else {
153 far[ii] = fbad ; /* clear out such criminal voxels in the result array */
154 }
155 }
156
157 /*----- process p-values to get q-values -----*/
158
159 if( nq > 19 ){ /* something to process! */
160
161 if( nvox-nq > 3333 ){ /* free up space, if significant */
162 qq = (float *)realloc( qq , sizeof(float)*nq ) ;
163 iq = (int *)realloc( iq , sizeof(int )*nq ) ;
164 }
165
166 STATUS("sorting p-values") ;
167 qsort_floatint( nq , qq , iq ) ; /* sort into increasing order; */
168 /* iq[] tracks where it's from */
169
170 /* downward scan from large p's to calculate q's */
171
172 qmin = 1.0 ;
173 nthr = (flags&1) ? nvox : nq ;
174 if( (flags&2) && nthr > 1 ){
175 float nold = nthr ;
176 nthr *= (logf(nthr)+0.5772157f+0.5f/nthr) ;
177 if( PRINT_TRACING ){
178 char str[256] ;
179 sprintf(str, "expand count %g to %g to allow for dependence",nold,nthr) ;
180 STATUS(str) ;
181 }
182 }
183 STATUS("convert p to q") ;
184 for( jj=nq-1 ; jj >= 0 ; jj-- ){ /* convert to q, then z */
185 qval = (nthr * qq[jj]) / (jj+1.0) ;
186 if( qval > qmin ) qval = qmin; else qmin = qval;
187 if( qsmal == 0 && qval <= QSTHRESH ) qsmal = jj ;
188 far[iq[jj]] = (float)qval ; /* store q into result array */
189 }
190
191 /** estimate NUMBER (not fraction) of true positives into mone (m1) **/
192
193 if( qsmal && qq[0] > 0.0f && (flags&1)==0 ) mone = (float)estimate_m1(nq,qq) ;
194
195 /* 26 Mar 2009: scale q down using m1 estimate above */
196
197 if( mone > 0.0f && !AFNI_yesenv("AFNI_DONT_ADJUST_FDR") ){
198 qfac = (nq-mone)/(float)nq; if( qfac < 0.5f ) qfac = 0.25f+qfac*qfac ;
199 /* ININFO_message("FDR: using mone = %g qfac = %g",mone,qfac) ; */
200 if( PRINT_TRACING ){
201 char str[256] ; sprintf(str,"Adjust q by %.3f",qfac) ; STATUS(str) ;
202 }
203 for( jj=0 ; jj < nq ; jj++ ) far[iq[jj]] *= qfac ;
204 }
205
206 /* convert to z-score, if ordered to do so */
207
208 if( doz ){
209 STATUS("convert to z") ;
210 for( jj=0 ; jj < nq ; jj++ ){
211 qval = (double)far[iq[jj]] ;
212 if( qval < QBOT ) qval = ZTOP ; /* honking big z-score */
213 else if( qval >= 1.0 ) qval = 0.0 ; /* very non-significant */
214 else qval = QTOZ(qval) ; /* meaningful z(q) */
215 far[iq[jj]] = (float)qval ;
216 }
217 }
218
219 free(iq) ; iq = NULL ; /* done with the original indexes */
220
221 /* compute missed detection fraction (MDF) vs. log10(p) */
222
223 if( mone > 8.0f ){
224 float ms , dpl ; int jtop,jbot , npp , kk ;
225 float *mdf=(float *)malloc(sizeof(float)*nq) ;
226 floatvec *fv ; float p1,p2,m1,m2,pf,mf , pl,pv ;
227 if( mdf == NULL ){
228 ERROR_message("mri_fdrize: out of memory!") ; goto finished ;
229 }
230
231 STATUS("computing mdf") ;
232 qmin = 1.0 ;
233 for( jj=nq-1 ; jj >=0 ; jj-- ){ /* scan down again to get q */
234 qval = (nthr * qq[jj]) / (jj+1.0) ; /* NOT adjusted by qfac! */
235 if( qval > qmin ) qval = qmin; else qmin = qval;
236 if( qval < QBOT ) qval = QBOT;
237
238 /* The reasoning involved to get MDF as a function of threshold:
239 * Number of values more signifcant than this threshold is jj+1;
240 * Approximately qval*(jj+1) of these are false positive detections
241 (that's what FDR means, dude or dudette);
242 * So about (1-qval)*(jj+1) are true positive detections;
243 * So about (1-qval)*(jj+1)/mone is the ratio
244 of true detections to the number of true positives;
245 * So about 1-(1-qval)*(jj+1)/mone is about the
246 fraction of missed true detections (MDF);
247 * The whole thing depends on the accuracy of the mone estimate
248 of the number of true positives (m1=mone) in the data.
249 * If you believe this, I've got a bridge to Brooklyn for sale! */
250
251 mdf[jj] = 1.0 - (1.0-qval)*(jj+1) / mone ;
252 if( mdf[jj] < 0.0f ) mdf[jj] = 0.0f ; /* make sure mdf */
253 else if( mdf[jj] > 1.0f ) mdf[jj] = 1.0f ; /* is reasonable */
254 }
255
256 /* However, MDF should only decrease as p-value increases */
257
258 for( jj=1 ; jj < nq ; jj++ ){
259 if( mdf[jj] > mdf[jj-1] ) mdf[jj] = mdf[jj-1] ;
260 }
261
262 /* cheapo trick: adjust MDF downwards to make sure it -> 0 as p -> 1 */
263
264 ms = mdf[nq-1] ; /* last MDF, nearest to p=1 */
265 if( ms > 0.0f ){
266 float alp=1.0f/(1.0f-ms) , bet=alp*ms ;
267 for( jj=0 ; jj < nq ; jj++ ) mdf[jj] = alp*mdf[jj]-bet ;
268 mdf[nq-1] = 0.0f ;
269 }
270
271 /* now find last nonzero mdf */
272
273 for( jj=nq-2 ; jj > 0 && mdf[jj] == 0.0f ; jj-- ) ; /*nada*/
274 if( jj <= 1 || qq[jj+1] <= qq[0] ){
275 if( jj <= 1 ) STATUS("all mdf zero? ==> no mdf") ;
276 if( qq[jj+1] <= qq[0] ) STATUS("qq=const? ==> no mdf") ;
277 free(mdf) ; goto finished ;
278 }
279 jtop = jj+1 ; /* mdf[jtop] = 0 */
280
281 /* now find first mdf below 99.9% */
282
283 for( jj=1 ; jj < jtop && mdf[jj] >= 0.999f ; jj++ ) ; /*nada*/
284 jbot = (jj < jtop-9) ? jj-1 : 0 ;
285
286 /* build a table of mdf vs log10(p) */
287
288 ms = log10( qq[jtop] / qq[jbot] ) ; /* will be positive */
289 npp = (int)( 0.99f + 5.0f * ms ) ; /* number of grid points */
290
291 if( npp < 3 ){
292 if( PRINT_TRACING ){
293 char str[256] ;
294 sprintf(str,"nq=%d npp=%d jbot=%d jtop=%d qq[jtop]=%g qq[jbot]=%g ==> no mdf",
295 nq , npp , jbot,jtop , qq[jtop] , qq[jbot] ) ;
296 STATUS(str) ;
297 }
298 free(mdf); goto finished;
299 }
300
301 dpl = ms / (npp-1) ; /* grid spacing in log10(p) */
302 MAKE_floatvec(fv,npp) ; fv->x0 = log10(qq[jbot]) ; fv->dx = dpl ;
303 fv->ar[0] = mdf[jbot] ;
304 for( jj=jbot,kk=1 ; kk < npp-1 ; kk++ ){
305 pl = fv->x0 + kk*dpl ; /* kk-th grid point in log10(p) */
306 pv = powf(10.0f,pl) ; /* kk-th grid point in p */
307 for( ; jj < jtop && qq[jj] < pv ; jj++ ) ; /*nada*/
308 /* linearly interpolate mdf in log10(p) */
309 p1 = log10(qq[jj-1]) ; p2 = log10(qq[jj]) ; pf = (pl-p1)/(p2-p1) ;
310 m1 = mdf[jj-1] ; m2 = mdf[jj] ; mf = pf*m2 + (1.0f-pf)*m1;
311 fv->ar[kk] = mf ;
312 }
313 fv->ar[npp-1] = 0.0f ;
314 FDR_mdfv = fv ; /* record the results for posterity */
315 if( PRINT_TRACING ){
316 char str[256]; sprintf(str,"MDF: npp=%d x0=%g dx=%g",npp,fv->x0,fv->dx);
317 STATUS(str) ;
318 }
319 free(mdf) ; /* end of producing mdf values */
320
321 } else {
322 if( !qsmal ) STATUS("no qsmal ==> no mdf") ;
323 if( nq < 200 ) STATUS("small nq ==> no mdf") ;
324 if( qq[0] <= 0.0f ) STATUS("bad qq[0] ==> no mdf") ;
325 }
326
327 } else {
328 /* warn user of impending doom 28 Jan, 2015 [rickr] */
329 WARNING_message("mri_fdrize: will not process only %d values (min=20)",
330 nq);
331 RETURN(0);
332 } /* end of producing q-values */
333
334 finished:
335 STATUS("finished") ;
336 if( iq != NULL ) free(iq);
337 if( qq != NULL ) free(qq);
338 RETURN(nq);
339 }
340
341 /*--------------------------------------------------------------------------*/
342
343 #undef NCURV
344 #define NCURV 101
345
346 /*! Create a curve that gives the FDR z(q) value vs. the statistical
347 threshold. Stored as a floatvec struct, with the statistical
348 value given as x0+i*dx and the corresponding z(q) value in ar[i]. */
349
mri_fdr_curve(MRI_IMAGE * im,int statcode,float * stataux)350 floatvec * mri_fdr_curve( MRI_IMAGE *im, int statcode, float *stataux )
351 {
352 MRI_IMAGE *cim ;
353 float *car , *far , *zar , *tar ;
354 int nvox , ii , nq , *iq ;
355 floatvec *fv=NULL ;
356 float tbot,ttop , zbot,ztop , dt,tt,zz , t1,t2,z1,z2 , tf,zf ;
357 int kk,klast,jj , flags=0 ;
358
359 ENTRY("mri_fdr_curve") ;
360
361 /* check for legal inputs */
362
363 KILL_floatvec(FDR_mdfv) ; /* erase the past */
364
365 if( !FUNC_IS_STAT(statcode) ) RETURN(NULL) ;
366 if( im == NULL || im->kind != MRI_float ) RETURN(NULL) ;
367 far = MRI_FLOAT_PTR(im); if( far == NULL ) RETURN(NULL) ;
368
369 /* make a copy of the statistics and convert them to z(q) scores */
370
371 STATUS("copy statistics image") ;
372 cim = mri_to_float(im) ; car = MRI_FLOAT_PTR(cim) ;
373
374 if( AFNI_yesenv("AFNI_NON_INDEPENDENT_FDR") ) flags = 2 ; /* 21 May 2012 */
375
376 nq = mri_fdrize( cim , statcode , stataux , flags ) ;
377 if( nq < 9 ){ mri_free(cim); RETURN(NULL); } /* bad FDR-izing */
378
379 /* create iq[] = list of voxels that were validly FDR-ized */
380
381 STATUS("make list of valid q's") ;
382 nvox = im->nvox ;
383 iq = (int *)malloc(sizeof(int)*nvox) ;
384 if( iq == NULL ){
385 ERROR_message("mri_fdr_curve: out of memory!") ;
386 mri_free(cim); RETURN(NULL);
387 }
388
389 for( nq=ii=0 ; ii < nvox ; ii++ ){ /* make list of voxels with */
390 if( car[ii] > 0.0f ) iq[nq++] = ii ; /* meaningful z(q) values */
391 }
392 if( nq < 9 ){ free(iq); mri_free(cim); RETURN(NULL); } /* bad */
393
394 /* create list of z(q) scores and corresponding statistics */
395
396 zar = (float *)malloc(sizeof(float)*nq) ; /* z(q) values */
397 tar = (float *)malloc(sizeof(float)*nq) ; /* statistics */
398 if( zar == NULL || tar == NULL ){
399 ERROR_message("mri_fdr_curve: out of memory!") ;
400 if( zar != NULL ) free(zar) ;
401 if( tar != NULL ) free(tar) ;
402 free(iq); mri_free(cim); RETURN(NULL);
403 }
404
405 STATUS("make list of z(q) vs. statistic") ;
406 for( ii=0 ; ii < nq ; ii++ ){
407 zar[ii] = car[iq[ii]] ; tar[ii] = fabsf(far[iq[ii]]) ;
408 }
409 free(iq) ; mri_free(cim) ; /* toss the trash */
410
411 STATUS("sort the list") ;
412 qsort_floatfloat( nq , zar , tar ) ; /* sort into increasing z */
413
414 /* find the largest z(q) that isn't beyond the top value we like */
415
416 for( klast=nq-1 ; klast > 0 && zar[klast] >= ZTOP ; klast-- ) ; /*nada */
417 if( klast == 0 ){ free(tar); free(zar); RETURN(NULL); }
418 if( klast < nq-1 ) klast++ ;
419
420 /* make the floatvec, evenly spaced in the statistic (tar) values */
421
422 tbot = tar[0] ; ttop = tar[klast] ; dt = (ttop-tbot)/(NCURV-1) ;
423 zbot = zar[0] ; ztop = zar[klast] ;
424
425 STATUS("make the floatvec") ;
426 MAKE_floatvec(fv,NCURV) ; fv->dx = dt ; fv->x0 = tbot ;
427 fv->ar[0] = zbot ;
428 for( jj=ii=1 ; ii < NCURV-1 ; ii++ ){
429 tt = tbot + ii*dt ; /* the statistic for this point on the curve */
430 for( ; jj < nq && tar[jj] < tt ; jj++ ) ; /*nada*/
431 t1 = tar[jj-1] ; t2 = tar[jj] ; tf = (tt-t1)/(t2-t1) ; /* linearly */
432 z1 = zar[jj-1] ; z2 = zar[jj] ; zf = tf*z2 + (1.0f-tf)*z1 ; /* interp z */
433 fv->ar[ii] = zf ; /* to this tt */
434 }
435 fv->ar[NCURV-1] = ztop ;
436
437 STATUS("finished") ;
438 free(tar) ; free(zar) ; RETURN(fv) ;
439 }
440