1 /*---------------------------------------------------------------------------*/
2 /*
3 This file contains routines for performing deconvolution analysis.
4
5 File: Deconvolve.c
6 Author: B. Douglas Ward
7 Date: 31 August 1998
8
9 Mod: Restructured matrix calculations to improve execution speed.
10 Date: 16 December 1998
11
12 Mod: Report mean square error from full model.
13 Date: 04 January 1999
14
15 Mod: Earlier termination if unable to invert X matrix.
16 (Avoids redundant error messages.)
17 Date: 06 January 1999
18
19 Mod: Modifications for matrix calculation of general linear tests.
20 Date: 02 July 1999
21
22 Mod: Additional statistical output (partial R^2 statistics).
23 Date: 07 September 1999
24
25 Mod: Allow reading of multiple input stimulus functions from a single
26 file by selection of individual columns.
27 Date: 09 November 1999
28
29 Mod: Added options for writing the fitted full model time series (-fitts)
30 and the residual error time series (-errts) to 3d+time datasets.
31 Date: 22 November 1999
32
33 Mod: Added -censor option to allow operator to eliminate individual
34 time points from the analysis.
35 Date: 21 June 2000
36
37 Mod: Added -censor option to allow operator to eliminate individual
38 time points from the analysis.
39 Date: 21 June 2000
40
41 Mod: Added screen display of p-values.
42 Date: 22 June 2000
43
44 Mod: Added -glt_label option for labeling the GLT matrix.
45 Date: 23 June 2000
46
47 Mod: Added -concat option for analysis of a concatenated 3d+time dataset.
48 Date: 26 June 2000
49
50 Mod: Increased size of screen output buffer.
51 Date: 27 July 2000
52
53 Mod: Additional output with -nodata option (norm.std.dev.'s for
54 GLT linear constraints).
55 Date: 11 August 2000
56
57 Mod: Added -stim_nptr option, to allow input stim functions to be sampled
58 at a multiple of the 1/TR rate.
59 Date: 02 January 2001
60
61 Mod: Changes to eliminate constraints on number of stimulus time series,
62 number of regressors, number of GLT's, and number of GLT linear
63 constraints.
64 Date: 10 May 2001
65
66 Mod: Made fstat_t2p() a static function to avoid conflicts on CYGWIN.
67 Date: 08 Jan 2002 - RWCox
68
69 Mod: Added static function student_t2p() for display of p-values
70 corresponding to the t-statistics.
71 Date: 28 January 2002
72
73 Mod: Modifications to glt_analysis and report_results for enhanced screen
74 output: Display of p-values for individual stim function regression
75 coefficients. Display of t-stats and p-values for individual linear
76 constraints within a GLT.
77 Date: 29 January 2002
78
79 Mod: Allow user to specify which input stimulus functions are part of
80 the baseline model.
81 Date: 02 May 2002
82
83 Mod: Increased size of screen output buffer (again).
84 Date: 02 December 2002
85
86 Mod: global variable legendre_polort replaces x^m with P_m(x)
87 Date 15 Jul 2004 - RWCox
88
89 */
90
91 /*---------------------------------------------------------------------------*/
92 /*
93 Include linear regression analysis software.
94 */
95
96 #include "RegAna.c"
97 #include "misc_math.h" /* 21 Jun 2010 [rickr] */
98
99
100 /*---------------------------------------------------------------------------*/
101 static int legendre_polort = 1 ; /* 15 Jul 2004 */
102 static int demean_base = 1 ; /* 12 Aug 2004 */
103
104 #define SPOL ((legendre_polort) ? "P_" : "t^") /* string label for polynomials */
105
106 /*---------------------------------------------------------------------------*/
107
108 #undef IBOT
109 #undef ITOP
110 #ifdef USE_BASIS
111 # define IBOT(ss) ((basis_stim[ss]!=NULL) ? 0 : min_lag[ss])
112 # define ITOP(ss) ((basis_stim[ss]!=NULL) ? basis_stim[ss]->nparm-1 : max_lag[ss])
113 #else
114 # define IBOT(ss) min_lag[ss]
115 # define ITOP(ss) max_lag[ss]
116 #endif
117
118 static matrix xfull ; /* X matrix for init_indep_var_matrix() below */
119 static int is_xfull=0 ;
120
121 static matrix xfull_plus ; /* X matrix with 0-1 augmentation [16 Aug 2019] */
122 static int is_xfull_plus=0 ;
123
124 /*---------------------------------------------------------------------------*/
125 /*
126 Initialize independent variable X matrix
127 */
128
init_indep_var_matrix(int p,int qp,int polort,int nt,int N,int * good_list,int * block_list,int num_blocks,int num_stimts,float ** stimulus,int * stim_length,int * min_lag,int * max_lag,int * nptr,int * stim_base,matrix * xgood)129 int init_indep_var_matrix
130 (
131 int p, /* number of parameters in the full model */
132 int qp, /* number of poly. trend baseline parameters */
133 int polort, /* degree of polynomial for baseline model */
134 int nt, /* number of images in input 3d+time dataset */
135 int N, /* total number of images used in the analysis */
136 int * good_list, /* list of usable time points */
137 int * block_list, /* list of block (run) starting points */
138 int num_blocks, /* number of blocks (runs) */
139 int num_stimts, /* number of stimulus time series */
140 float ** stimulus, /* stimulus time series arrays */
141 int * stim_length, /* length of stimulus time series arrays */
142 int * min_lag, /* minimum time delay for impulse response */
143 int * max_lag, /* maximum time delay for impulse response */
144 int * nptr, /* number of stim fn. time points per TR */
145 int * stim_base , /* flag for being in the baseline [12 Aug 2004] */
146 matrix * xgood /* independent variable matrix */
147 )
148
149 {
150 int m; /* X matrix column index */
151 int n; /* X matrix row index */
152 int is; /* input stimulus index */
153 int ilag; /* time lag index */
154 int ib; /* block (run) index */
155 int nfirst, nlast; /* time boundaries of a block (run) */
156 int mfirst=0, mlast=0; /* column boundaries of baseline parameters
157 for a block (run) */
158
159 float * stim_array; /* stimulus function time series */
160
161 int mold ; /* 12 Aug 2004 */
162 int ibot,itop ;
163
164 ENTRY("init_indep_var_matrix") ;
165
166
167 /*----- Initialize X matrix -----*/
168
169 if( is_xfull ){ matrix_destroy(&xfull); is_xfull=0; }
170
171 if( is_xfull_plus ){ matrix_destroy(&xfull_plus); is_xfull_plus=0; }
172
173 STATUS("create x matrix" ) ;
174 matrix_initialize (&xfull);
175 matrix_create (nt, p, &xfull);
176
177
178 /*----- Set up columns of X matrix corresponding to
179 the baseline (null hypothesis) signal model -----*/
180
181 STATUS("loop over blocks") ;
182
183 for (ib = 0; ib < num_blocks; ib++)
184 {
185 nfirst = block_list[ib]; /* start time index for this run */
186 if (ib+1 < num_blocks)
187 nlast = block_list[ib+1]; /* last+1 time index for this run */
188 else
189 nlast = nt;
190
191 if(PRINT_TRACING){
192 char str[256] ;
193 sprintf(str," block #%d = %d .. %d",ib,nfirst,nlast-1) ; STATUS(str) ;
194 }
195
196 for (n = nfirst; n < nlast; n++)
197 {
198 mfirst = ib * (polort+1); /* first column index */
199 mlast = (ib+1) * (polort+1); /* last+1 column index */
200
201 if( !legendre_polort ){ /* the old way: powers */
202 for (m = mfirst; m < mlast; m++)
203 xfull.elts[n][m] = pow ((double)(n-nfirst), (double)(m-mfirst));
204
205 } else { /* 15 Jul 2004: the new way: Legendre - RWCox */
206
207 double xx , aa=2.0/(nlast-nfirst-1.0) ; /* map nfirst..nlast-1 */
208 for( m=mfirst ; m < mlast ; m++ ){ /* to interval [-1,1] */
209 xx = aa*(n-nfirst) - 1.0 ;
210 xfull.elts[n][m] = legendre( xx , m-mfirst ) ;
211 }
212 }
213 }
214
215 if( mfirst+1 < mlast && demean_base ){ /* 12 Aug 2004: remove means? */
216 float sum ;
217 for( m=mfirst+1 ; m < mlast ; m++ ){
218 sum = 0.0f ;
219 for( n=nfirst ; n < nlast ; n++ ) sum += xfull.elts[n][m] ;
220 sum /= (nlast-nfirst) ;
221 for( n=nfirst ; n < nlast ; n++ ) xfull.elts[n][m] -= sum ;
222 }
223 }
224 }
225
226
227 /*----- Set up columns of X matrix corresponding to
228 time delayed versions of the input stimulus -----*/
229
230 STATUS("loop over stimulus time series") ;
231
232 m = qp;
233 for (is = 0; is < num_stimts; is++){
234 #ifdef USE_BASIS
235 if( basis_vect[is] != NULL ){ /* 16 Aug 2004 */
236 float *bv=MRI_FLOAT_PTR(basis_vect[is]) ;
237 int nf=basis_vect[is]->ny , jj ;
238 if( PRINT_TRACING ){
239 char str[256] ;
240 sprintf(str," stim #%d: expanding into %d basis vectors",is,nf) ;
241 STATUS(str) ;
242 }
243 for( jj=0 ; jj < nf ; jj++ ){
244 for( n=0 ; n < nt ; n++ ) xfull.elts[n][m] = bv[n+jj*nt] ;
245 m++ ;
246 }
247 }
248 else {
249 #endif
250 if (stim_length[is] < nt*nptr[is])
251 {
252 DC_error ("Input stimulus time series is too short");
253 RETURN (0);
254 }
255 stim_array = stimulus[is]; mold = m ; /* mold = col index we start at */
256 ibot = IBOT(is) ; itop = ITOP(is) ;
257 if( PRINT_TRACING ){
258 char str[256] ;
259 sprintf(str," stim #%d: ibot=%d itop=%d",is,ibot,itop) ;
260 STATUS(str) ;
261 }
262 for( ilag=ibot ; ilag <= itop ; ilag++ )
263 {
264 for (n = 0; n < nt; n++)
265 {
266 if (n*nptr[is] < ilag)
267 xfull.elts[n][m] = 0.0;
268 else
269 xfull.elts[n][m] = stim_array[n*nptr[is]-ilag];
270 }
271 m++;
272 }
273
274 /* 12 Aug 2004: remove mean of baseline columns? */
275 /* 07 Feb 2005: Oops -- used [m] instead of [mm] in the for(n) loops! */
276
277 if( stim_base != NULL && stim_base[is] && demean_base ){
278 int mm ; float sum ;
279 STATUS(" remove baseline mean") ;
280 for( mm=mold ; mm < m ; mm++ ){
281 sum = 0.0f ;
282 for( n=0 ; n < nt ; n++ ) sum += xfull.elts[n][mm] ;
283 sum /= nt ;
284 for( n=0 ; n < nt ; n++ ) xfull.elts[n][mm] -= sum ;
285 }
286 }
287 #ifdef USE_BASIS
288 }
289 #endif
290 }
291
292
293 /*----- Keep only those rows of the X matrix which correspond to
294 usable time points (i.e., apply censoring here) -----*/
295
296 STATUS("extract xgood matrix") ;
297
298 matrix_extract_rows (xfull, N, good_list, xgood);
299
300 is_xfull = 1 ; /* original (uncensored) X matrix saved in xfull */
301
302 /* Augment full X matrix with 0-1 columns: censor via regression [16 Aug 2019] */
303
304 { int ii , nbad , *glist , *blist ;
305 glist = (int *)calloc(sizeof(int),nt) ; /* 0-1 list of good points */
306 blist = (int *)calloc(sizeof(int),nt) ; /* index list of bad points */
307 for( ii=0 ; ii < N ; ii++ ){ glist[good_list[ii]] = 1 ; }
308 for( nbad=ii=0 ; ii < nt ; ii++ ){ if( glist[ii] == 0 ) blist[nbad++] = ii ; }
309 matrix_augment_01_columns( xfull, nbad , blist , &xfull_plus ) ;
310 is_xfull_plus = 1 ;
311 free(blist) ; free(glist) ;
312 }
313
314 RETURN (1);
315 }
316
317
318 /*---------------------------------------------------------------------------*/
319 /*
320 Initialization for the regression analysis.
321 */
322
init_regression_analysis(int p,int qp,int num_stimts,int * baseline,int * min_lag,int * max_lag,matrix xdata,matrix * x_full,matrix * xtxinv_full,matrix * xtxinvxt_full,matrix * x_base,matrix * xtxinvxt_base,matrix * x_rdcd,matrix * xtxinvxt_rdcd)323 int init_regression_analysis
324 (
325 int p, /* number of parameters in full model */
326 int qp, /* number of poly. trend baseline parameters */
327 int num_stimts, /* number of stimulus time series */
328 int * baseline, /* flag for stim function in baseline model */
329 int * min_lag, /* minimum time delay for impulse response */
330 int * max_lag, /* maximum time delay for impulse response */
331 matrix xdata, /* independent variable matrix */
332 matrix * x_full, /* extracted X matrix for full model */
333 matrix * xtxinv_full, /* matrix: 1/(X'X) for full model */
334 matrix * xtxinvxt_full, /* matrix: (1/(X'X))X' for full model */
335 matrix * x_base, /* extracted X matrix for baseline model */
336 matrix * xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */
337 matrix * x_rdcd, /* extracted X matrices for reduced models */
338 matrix * xtxinvxt_rdcd /* matrix: (1/(X'X))X' for reduced models */
339 )
340
341 {
342 int * plist = NULL; /* list of model parameters */
343 int ip, it; /* parameter indices */
344 int is, js; /* stimulus indices */
345 int im, jm; /* lag index */
346 int ok; /* flag for successful matrix calculation */
347 matrix xtxinv_temp; /* intermediate results */
348 int ibot,itop ;
349
350
351 ENTRY("init_regression_analysis") ;
352
353 /*----- Initialize matrix -----*/
354 matrix_initialize (&xtxinv_temp);
355
356
357 /*----- Initialize matrices for the baseline model -----*/
358 plist = (int *) malloc (sizeof(int) * p); MTEST(plist);
359 for (ip = 0; ip < qp; ip++)
360 plist[ip] = ip;
361 it = ip = qp;
362 for (is = 0; is < num_stimts; is++)
363 {
364 ibot = IBOT(is) ; itop = ITOP(is) ;
365 for (im = ibot; im <= itop; im++)
366 {
367 if (baseline[is])
368 {
369 plist[ip] = it;
370 ip++;
371 }
372 it++;
373 }
374 }
375 ok = calc_matrices (xdata, ip, plist, x_base, &xtxinv_temp, xtxinvxt_base);
376 if (!ok) { matrix_destroy (&xtxinv_temp); RETURN (0); };
377
378
379 /*----- Initialize matrices for stimulus functions -----*/
380 for (is = 0; is < num_stimts; is++)
381 {
382 for (ip = 0; ip < qp; ip++)
383 {
384 plist[ip] = ip;
385 }
386
387 it = ip = qp;
388
389 for (js = 0; js < num_stimts; js++)
390 {
391 ibot = IBOT(js) ; itop = ITOP(js) ;
392 for (jm = ibot; jm <= itop; jm++)
393 {
394 if (is != js){ plist[ip] = it; ip++; }
395 it++;
396 }
397 }
398
399 ibot = IBOT(is) ; itop = ITOP(is) ;
400 ok = calc_matrices (xdata, p-(itop-ibot+1),
401 plist, &(x_rdcd[is]), &xtxinv_temp, &(xtxinvxt_rdcd[is]));
402 if (!ok) { matrix_destroy (&xtxinv_temp); RETURN (0); };
403 }
404
405
406 /*----- Initialize matrices for full model -----*/
407 for (ip = 0; ip < p; ip++)
408 plist[ip] = ip;
409 ok = calc_matrices (xdata, p, plist, x_full, xtxinv_full, xtxinvxt_full);
410 if (!ok) { matrix_destroy (&xtxinv_temp); RETURN (0); };
411
412
413 /*----- Destroy matrix -----*/
414 matrix_destroy (&xtxinv_temp);
415
416 if (plist != NULL) { free(plist); plist = NULL; }
417
418 RETURN (1);
419 }
420
421
422 /*---------------------------------------------------------------------------*/
423 /*
424 Initialization for the general linear test analysis.
425 */
426
init_glt_analysis(matrix xtxinv,int glt_num,matrix * glt_cmat,matrix * glt_amat,matrix * cxtxinvct)427 int init_glt_analysis
428 (
429 matrix xtxinv, /* matrix: 1/(X'X) for full model */
430 int glt_num, /* number of general linear tests */
431 matrix * glt_cmat, /* general linear test matrices */
432 matrix * glt_amat, /* constant GLT matrices for later use */
433 matrix * cxtxinvct /* matrices: C(1/(X'X))C' for GLT */
434
435 )
436
437 {
438 int iglt; /* index for general linear test */
439 int ok; /* flag for successful matrix inversion */
440
441
442 ENTRY("init_glt_analysis") ;
443
444 for (iglt = 0; iglt < glt_num; iglt++)
445 {
446 ok = calc_glt_matrix (xtxinv, glt_cmat[iglt], &(glt_amat[iglt]),
447 &(cxtxinvct[iglt]));
448 if (! ok) RETURN (0);
449 }
450
451 RETURN (1);
452 }
453
454
455 /*---------------------------------------------------------------------------*/
456 /*
457 Calculate results for this voxel.
458 */
459
regression_analysis(int N,int p,int q,int num_stimts,int * min_lag,int * max_lag,matrix x_full,matrix xtxinv_full,matrix xtxinvxt_full,matrix x_base,matrix xtxinvxt_base,matrix * x_rdcd,matrix * xtxinvxt_rdcd,vector y,float rms_min,float * mse,vector * coef_full,vector * scoef_full,vector * tcoef_full,float * fpart,float * rpart,float * ffull,float * rfull,int * novar,float * fitts,float * errts)460 void regression_analysis
461 (
462 int N, /* number of usable data points */
463 int p, /* number of parameters in full model */
464 int q, /* number of parameters in baseline model */
465 int num_stimts, /* number of stimulus time series */
466 int * min_lag, /* minimum time delay for impulse response */
467 int * max_lag, /* maximum time delay for impulse response */
468 matrix x_full, /* extracted X matrix for full model */
469 matrix xtxinv_full, /* matrix: 1/(X'X) for full model */
470 matrix xtxinvxt_full, /* matrix: (1/(X'X))X' for full model */
471 matrix x_base, /* extracted X matrix for baseline model */
472 matrix xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */
473 matrix * x_rdcd, /* extracted X matrices for reduced models */
474 matrix * xtxinvxt_rdcd, /* matrix: (1/(X'X))X' for reduced models */
475 vector y, /* vector of measured data */
476 float rms_min, /* minimum variation in data to fit full model */
477 float * mse, /* mean square error from full model */
478 vector * coef_full, /* regression parameters */
479 vector * scoef_full, /* std. devs. for regression parameters */
480 vector * tcoef_full, /* t-statistics for regression parameters */
481 float * fpart, /* partial F-statistics for the stimuli */
482 float * rpart, /* partial R^2 stats. for the stimuli */
483 float * ffull, /* full model F-statistics */
484 float * rfull, /* full model R^2 stats. */
485 int * novar, /* flag for insufficient variation in data */
486 float * fitts, /* full model fitted time series */
487 float * errts /* full model residual error time series */
488 )
489
490 {
491 int is; /* input stimulus index */
492 float sse_base; /* error sum of squares, baseline model */
493 float sse_rdcd; /* error sum of squares, reduced model */
494 float sse_full; /* error sum of squares, full model */
495 vector coef_temp; /* intermediate results */
496 int ibot,itop ;
497
498
499 ENTRY("regression_analysis") ;
500
501 /*----- Initialization -----*/
502 vector_initialize (&coef_temp);
503
504
505 /*----- Calculate regression coefficients for baseline model -----*/
506 calc_coef (xtxinvxt_base, y, &coef_temp);
507
508
509 /*----- Calculate the error sum of squares for the baseline model -----*/
510 sse_base = calc_sse (x_base, coef_temp, y);
511
512
513 /*----- Stop here if variation about baseline is sufficiently low -----*/
514 if (sqrt(sse_base/N) < rms_min)
515 {
516 int it;
517
518 *novar = 1;
519 vector_create (p, coef_full);
520 vector_create (p, scoef_full);
521 vector_create (p, tcoef_full);
522 for (is = 0; is < num_stimts; is++)
523 {
524 fpart[is] = 0.0;
525 rpart[is] = 0.0;
526 }
527 for (it = 0; it < N; it++)
528 {
529 fitts[it] = 0.0;
530 errts[it] = 0.0;
531 }
532 *mse = 0.0;
533 *rfull = 0.0;
534 *ffull = 0.0;
535 vector_destroy (&coef_temp);
536 EXRETURN;
537 }
538 else
539 *novar = 0;
540
541
542 /*----- Calculate regression coefficients for the full model -----*/
543 calc_coef (xtxinvxt_full, y, coef_full);
544
545
546 /*----- Calculate the error sum of squares for the full model -----*/
547 sse_full = calc_sse_fit (x_full, *coef_full, y, fitts, errts);
548 *mse = sse_full / (N-p);
549
550
551 /*----- Calculate t-statistics for the regression coefficients -----*/
552 calc_tcoef (N, p, sse_full, xtxinv_full,
553 *coef_full, scoef_full, tcoef_full);
554
555
556 /*----- Determine significance of the individual stimuli -----*/
557 for (is = 0; is < num_stimts; is++)
558 {
559
560 /*----- Calculate regression coefficients for reduced model -----*/
561 calc_coef (xtxinvxt_rdcd[is], y, &coef_temp);
562
563
564 /*----- Calculate the error sum of squares for the reduced model -----*/
565 sse_rdcd = calc_sse (x_rdcd[is], coef_temp, y);
566
567
568 /*----- Calculate partial F-stat for significance of the stimulus -----*/
569 ibot = IBOT(is) ; itop = ITOP(is) ;
570 fpart[is] = calc_freg (N, p, p-(itop-ibot+1), sse_full, sse_rdcd);
571
572
573 /*----- Calculate partial R^2 for this stimulus -----*/
574 rpart[is] = calc_rsqr (sse_full, sse_rdcd);
575
576 }
577
578
579 /*----- Calculate coefficient of multiple determination R^2 -----*/
580 *rfull = calc_rsqr (sse_full, sse_base);
581
582
583 /*----- Calculate the total regression F-statistic -----*/
584 *ffull = calc_freg (N, p, q, sse_full, sse_base);
585
586
587 /*----- Dispose of vector -----*/
588 vector_destroy (&coef_temp);
589
590 EXRETURN ;
591 }
592
593
594 /*---------------------------------------------------------------------------*/
595 /*
596 Perform the general linear test analysis for this voxel.
597 */
598
glt_analysis(int N,int p,matrix x,vector y,float ssef,vector coef,int novar,matrix * cxtxinvct,int glt_num,int * glt_rows,matrix * glt_cmat,matrix * glt_amat,vector * glt_coef,vector * glt_tcoef,float * fglt,float * rglt)599 void glt_analysis
600 (
601 int N, /* number of data points */
602 int p, /* number of parameters in the full model */
603 matrix x, /* X matrix for full model */
604 vector y, /* vector of measured data */
605 float ssef, /* error sum of squares from full model */
606 vector coef, /* regression parameters for full model */
607 int novar, /* flag for insufficient variation in data */
608 matrix * cxtxinvct, /* matrices: C(1/(X'X))C' for GLT */
609 int glt_num, /* number of general linear tests */
610 int * glt_rows, /* number of linear constraints in glt */
611 matrix * glt_cmat, /* general linear test matrices */
612 matrix * glt_amat, /* constant matrices */
613 vector * glt_coef, /* linear combinations from GLT matrices */
614 vector * glt_tcoef, /* t-statistics for GLT linear combinations */
615 float * fglt, /* F-statistics for the general linear tests */
616 float * rglt /* R^2 statistics for the general linear tests */
617 )
618
619 {
620 int iglt; /* index for general linear test */
621 int q; /* number of parameters in the rdcd model */
622 float sser; /* error sum of squares, reduced model */
623 vector rcoef; /* regression parameters for reduced model */
624 vector scoef; /* std. devs. for regression parameters */
625
626 ENTRY("glt_analysis") ;
627
628 /*----- Initialization -----*/
629 vector_initialize (&rcoef);
630 vector_initialize (&scoef);
631
632
633 /*----- Loop over multiple general linear tests -----*/
634 for (iglt = 0; iglt < glt_num; iglt++)
635 {
636 /*----- Test for insufficient variation in data -----*/
637 if (novar)
638 {
639 vector_create (glt_rows[iglt], &glt_coef[iglt]);
640 vector_create (glt_rows[iglt], &glt_tcoef[iglt]);
641 fglt[iglt] = 0.0;
642 rglt[iglt] = 0.0;
643 }
644 else
645 {
646 /*----- Calculate the GLT linear combinations -----*/
647 calc_lcoef (glt_cmat[iglt], coef, &glt_coef[iglt]);
648
649 /*----- Calculate t-statistics for GLT linear combinations -----*/
650 calc_tcoef (N, p, ssef, cxtxinvct[iglt],
651 glt_coef[iglt], &scoef, &glt_tcoef[iglt]);
652
653 /*----- Calculate regression parameters for the reduced model -----*/
654 /*----- (that is, the model in the column space of X but ) -----*/
655 /*----- (orthogonal to the restricted column space of XC') -----*/
656 calc_rcoef (glt_amat[iglt], coef, &rcoef);
657
658 /*----- Calculate error sum of squares for the reduced model -----*/
659 sser = calc_sse (x, rcoef, y);
660
661 /*----- Calculate the F-statistic for this GLT -----*/
662 q = p - glt_rows[iglt];
663 fglt[iglt] = calc_freg (N, p, q, ssef, sser);
664
665 /*----- Calculate the R^2 statistic for this GLT -----*/
666 rglt[iglt] = calc_rsqr (ssef, sser);
667
668 }
669 }
670
671
672 /*----- Dispose of vectors -----*/
673 vector_destroy (&rcoef);
674 vector_destroy (&scoef);
675
676 EXRETURN ;
677 }
678
679
680 /*---------------------------------------------------------------------------*/
681 /*
682 Convert t-values and F-values to p-value.
683 These routines were copied and modified from: mri_stats.c
684 16 May 2005: Change names from '_t2p' to '_t2pp' to avoid library
685 name conflict on Mac OS X 10.4 (stupid system).
686 */
687
688
student_t2pp(double tt,double dof)689 static double student_t2pp( double tt , double dof )
690 {
691 double bb , xx , pp ;
692
693 tt = fabs(tt);
694
695 if( dof < 1.0 ) return 1.0 ;
696
697 if (tt >= 1000.0) return 0.0;
698
699 bb = lnbeta( 0.5*dof , 0.5 ) ;
700 xx = dof/(dof + tt*tt) ;
701 pp = incbeta( xx , 0.5*dof , 0.5 , bb ) ;
702 return pp ;
703 }
704
705
fstat_t2pp(double ff,double dofnum,double dofden)706 static double fstat_t2pp( double ff , double dofnum , double dofden )
707 {
708 int which , status ;
709 double p , q , f , dfn , dfd , bound ;
710
711 if (ff >= 1000.0) return 0.0;
712
713 which = 1 ;
714 p = 0.0 ;
715 q = 0.0 ;
716 f = ff ;
717 dfn = dofnum ;
718 dfd = dofden ;
719
720 cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
721
722 if( status == 0 ) return q ;
723 else return 1.0 ;
724 }
725
726 /*---------------------------------------------------------------------------*/
727
728 static char lbuf[65536]; /* character string containing statistical summary */
729 static char sbuf[512];
730
731
report_results(int N,int qp,int q,int p,int polort,int * block_list,int num_blocks,int num_stimts,char ** stim_label,int * baseline,int * min_lag,int * max_lag,vector coef,vector tcoef,float * fpart,float * rpart,float ffull,float rfull,float mse,int glt_num,char ** glt_label,int * glt_rows,vector * glt_coef,vector * glt_tcoef,float * fglt,float * rglt,char ** label)732 void report_results
733 (
734 int N, /* number of usable data points */
735 int qp, /* number of poly. trend baseline parameters */
736 int q, /* number of baseline model parameters */
737 int p, /* number of full model parameters */
738 int polort, /* degree of polynomial for baseline model */
739 int * block_list, /* list of block (run) starting points */
740 int num_blocks, /* number of blocks (runs) */
741 int num_stimts, /* number of stimulus time series */
742 char ** stim_label, /* label for each stimulus */
743 int * baseline, /* flag for stim function in baseline model */
744 int * min_lag, /* minimum time delay for impulse response */
745 int * max_lag, /* maximum time delay for impulse response */
746 vector coef, /* regression parameters */
747 vector tcoef, /* t-statistics for regression parameters */
748 float * fpart, /* partial F-statistics for the stimuli */
749 float * rpart, /* partial R^2 stats. for the stimuli */
750 float ffull, /* full model F-statistic */
751 float rfull, /* full model R^2 stat. */
752 float mse, /* mean square error from full model */
753 int glt_num, /* number of general linear tests */
754 char ** glt_label, /* label for general linear tests */
755 int * glt_rows, /* number of linear constraints in glt */
756 vector * glt_coef, /* linear combinations from GLT matrices */
757 vector * glt_tcoef, /* t-statistics for GLT linear combinations */
758 float * fglt, /* F-statistics for the general linear tests */
759 float * rglt, /* R^2 statistics for the general linear tests */
760 char ** label /* statistical summary for ouput display */
761 )
762
763 {
764 const int MAXBUF = 65000; /* maximum buffer string length */
765 int m; /* coefficient index */
766 int is; /* stimulus index */
767 int ilag; /* time lag index */
768
769 int iglt; /* general linear test index */
770 int ilc; /* linear combination index */
771
772 double pvalue; /* p-value corresponding to F-value */
773 int r; /* number of parameters in the reduced model */
774
775 int ib; /* block (run) index */
776 int mfirst, mlast; /* column boundaries of baseline parameters
777 for a block (run) */
778 int ibot,itop ;
779
780
781 lbuf[0] = '\0' ; /* make this a 0 length string to start */
782
783 /** for each reference, make a string into sbuf **/
784
785
786 /*----- Statistical results for baseline fit -----*/
787 if (num_blocks == 1)
788 {
789 sprintf (sbuf, "\nBaseline: \n");
790 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf); else goto finisher ;
791 for (m=0; m < qp; m++)
792 {
793 sprintf (sbuf, "%s%d coef = %10.4f ", SPOL,m, coef.elts[m]);
794 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
795 sprintf (sbuf, "%s%d t-st = %10.4f ", SPOL,m, tcoef.elts[m]);
796 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
797 pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
798 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
799 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
800 }
801 }
802 else
803 {
804 for (ib = 0; ib < num_blocks; ib++)
805 {
806 sprintf (sbuf, "\nBaseline for Run #%d: \n", ib+1);
807 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf); else goto finisher ;
808
809 mfirst = ib * (polort+1);
810 mlast = (ib+1) * (polort+1);
811 for (m = mfirst; m < mlast; m++)
812 {
813 sprintf (sbuf, "%s%d coef = %10.4f ",
814 SPOL,m - mfirst, coef.elts[m]);
815 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
816 sprintf (sbuf, "%s%d t-st = %10.4f ",
817 SPOL,m - mfirst, tcoef.elts[m]);
818 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
819 pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
820 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
821 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
822 }
823 }
824 }
825
826
827 /*----- Statistical results for stimulus response -----*/
828 m = qp;
829 for (is = 0; is < num_stimts; is++)
830 {
831 if (baseline[is])
832 sprintf (sbuf, "\nBaseline: %s \n", stim_label[is]);
833 else
834 sprintf (sbuf, "\nStimulus: %s \n", stim_label[is]);
835 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf); else goto finisher ;
836 ibot = IBOT(is) ; itop = ITOP(is) ;
837 for (ilag = ibot; ilag <= itop; ilag++)
838 {
839 sprintf (sbuf,"h[%2d] coef = %10.4f ", ilag, coef.elts[m]);
840 if (strlen(lbuf) < MAXBUF) strcat(lbuf,sbuf) ; else goto finisher ;
841 sprintf (sbuf,"h[%2d] t-st = %10.4f ", ilag, tcoef.elts[m]);
842 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
843 pvalue = student_t2pp ((double)tcoef.elts[m], (double)(N-p));
844 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
845 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
846 m++;
847 }
848
849 sprintf (sbuf, " R^2 = %10.4f ", rpart[is]);
850 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
851 r = p - (itop-ibot+1);
852 sprintf (sbuf, "F[%2d,%3d] = %10.4f ", p-r, N-p, fpart[is]);
853 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
854 pvalue = fstat_t2pp ((double)fpart[is], (double)(p-r), (double)(N-p));
855 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
856 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
857 }
858
859
860 /*----- Statistical results for full model -----*/
861 sprintf (sbuf, "\nFull Model: \n");
862 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
863
864 sprintf (sbuf, " MSE = %10.4f \n", mse);
865 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
866
867 sprintf (sbuf, " R^2 = %10.4f ", rfull);
868 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
869
870 sprintf (sbuf, "F[%2d,%3d] = %10.4f ", p-q, N-p, ffull);
871 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
872 pvalue = fstat_t2pp ((double)ffull, (double)(p-q), (double)(N-p));
873 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
874 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
875
876
877 /*----- Statistical results for general linear test -----*/
878 if (glt_num > 0)
879 {
880 for (iglt = 0; iglt < glt_num; iglt++)
881 {
882 sprintf (sbuf, "\nGeneral Linear Test: %s \n", glt_label[iglt]);
883 if (strlen(lbuf) < MAXBUF) strcat (lbuf,sbuf); else goto finisher ;
884 for (ilc = 0; ilc < glt_rows[iglt]; ilc++)
885 {
886 sprintf (sbuf, "LC[%d] coef = %10.4f ",
887 ilc, glt_coef[iglt].elts[ilc]);
888 if (strlen(lbuf) < MAXBUF) strcat (lbuf,sbuf); else goto finisher ;
889 sprintf (sbuf, "LC[%d] t-st = %10.4f ",
890 ilc, glt_tcoef[iglt].elts[ilc]);
891 if (strlen(lbuf) < MAXBUF) strcat (lbuf,sbuf); else goto finisher ;
892 pvalue = student_t2pp ((double)glt_tcoef[iglt].elts[ilc],
893 (double)(N-p));
894 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
895 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
896 }
897
898 sprintf (sbuf, " R^2 = %10.4f ", rglt[iglt]);
899 if (strlen(lbuf) < MAXBUF) strcat (lbuf,sbuf); else goto finisher ;
900
901 r = p - glt_rows[iglt];
902 sprintf (sbuf, "F[%2d,%3d] = %10.4f ", p-r, N-p, fglt[iglt]);
903 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
904 pvalue = fstat_t2pp ((double)fglt[iglt],
905 (double)(p-r), (double)(N-p));
906 sprintf (sbuf, "p-value = %12.4e \n", pvalue);
907 if (strlen(lbuf) < MAXBUF) strcat (lbuf, sbuf); else goto finisher ;
908 }
909 }
910
911 finisher:
912 if (strlen(lbuf) >= MAXBUF)
913 strcat (lbuf, "\n\nWarning: Screen output buffer is full. \n");
914
915 *label = lbuf ; /* send address of lbuf back in what label points to */
916
917 }
918