Home > marsbar > @mardo_99 > mars_spm_graph.m

mars_spm_graph

PURPOSE ^

Graphical display of adjusted data

SYNOPSIS ^

function [r_st,marsD,changef] = mars_spm_graph(marsD,rno,Ic)

DESCRIPTION ^

 Graphical display of adjusted data
 FORMAT [r_st,marsD,changef] = mars_spm_graph(marsD,rno,Ic)

 marsD    - SPM design object
        required fields in des_struct are:
        xX     - Design Matrix structure
        - (see spm_spm.m for structure)
        betas  - the betas!
        ResidualMS  - the residual mean square
        xCon   - contrast definitions
          - required fields are:
          .c  - contrast vector/matrix
          (see spm_FcUtil.m for details of contrast structure... )
        marsY  - the data itself (as a data object)

 rno    - region number (index for marsD.marsY)
 Ic     - contrast number (optional)
 
 Returns
 r_st   - return structure, with fields
          Y      - fitted   data for the selected voxel
          y      - adjusted data for the selected voxel
          beta   - parameter estimates
          SE     - standard error of parameter estimates
          cbeta  = betas multiplied by contrast
 marsD   - design structure, with possibly added contrasts
 changef - set to 1 if design has changed

 see spm99 version of spm_graph for details
_______________________________________________________________________
 @(#)spm_graph.m    2.29 Karl Friston 00/09/14
 Modified for MarsBar - Matthew Brett 1/11/01 NRN
 Added return of cbetas at some point before 28/1/03

 $Id$

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [r_st,marsD,changef] = mars_spm_graph(marsD,rno,Ic)
0002 % Graphical display of adjusted data
0003 % FORMAT [r_st,marsD,changef] = mars_spm_graph(marsD,rno,Ic)
0004 %
0005 % marsD    - SPM design object
0006 %        required fields in des_struct are:
0007 %        xX     - Design Matrix structure
0008 %        - (see spm_spm.m for structure)
0009 %        betas  - the betas!
0010 %        ResidualMS  - the residual mean square
0011 %        xCon   - contrast definitions
0012 %          - required fields are:
0013 %          .c  - contrast vector/matrix
0014 %          (see spm_FcUtil.m for details of contrast structure... )
0015 %        marsY  - the data itself (as a data object)
0016 %
0017 % rno    - region number (index for marsD.marsY)
0018 % Ic     - contrast number (optional)
0019 %
0020 % Returns
0021 % r_st   - return structure, with fields
0022 %          Y      - fitted   data for the selected voxel
0023 %          y      - adjusted data for the selected voxel
0024 %          beta   - parameter estimates
0025 %          SE     - standard error of parameter estimates
0026 %          cbeta  = betas multiplied by contrast
0027 % marsD   - design structure, with possibly added contrasts
0028 % changef - set to 1 if design has changed
0029 %
0030 % see spm99 version of spm_graph for details
0031 %_______________________________________________________________________
0032 % @(#)spm_graph.m    2.29 Karl Friston 00/09/14
0033 % Modified for MarsBar - Matthew Brett 1/11/01 NRN
0034 % Added return of cbetas at some point before 28/1/03
0035 %
0036 % $Id$
0037 
0038 if ~is_mars_estimated(marsD)
0039   error('Need estimated design for plot');
0040 end
0041 if nargin < 2
0042   rno = [];
0043 end
0044 if nargin < 3
0045   Ic = [];
0046 end
0047 changef = 0;
0048 
0049 % make values ready for return
0050 def_r_st = struct(...
0051     'Y', [],...
0052     'y', [],...
0053     'beta', [],...
0054     'SE', [],...
0055     'cbeta',[]);
0056 cbeta = [];
0057 
0058 % get stuff from object
0059 mRes = des_struct(marsD);
0060 xCon = mRes.xCon;
0061 
0062 % Check if we want to, and can, assume region no is 1
0063 if isempty(rno) 
0064   if n_regions(mRes.marsY) > 1
0065     error('Need to specify region number');
0066   end
0067   rno = 1;
0068 end
0069 
0070 %-Get Graphics figure handle
0071 %-----------------------------------------------------------------------
0072 Fgraph = spm_figure('GetWin','Graphics');
0073 
0074 %-Delete previous axis and their pagination controls (if any)
0075 %-----------------------------------------------------------------------
0076 spm_results_ui('Clear',Fgraph,2);
0077 
0078 % Get required data
0079 sY   = summary_data(mRes.marsY); 
0080 y    = apply_filter(marsD, sY(:, rno));
0081 
0082 % Get design matrix for convenience
0083 xX = mRes.xX;
0084 
0085 % Label for region
0086 XYZstr = region_name(mRes.marsY, rno);
0087 XYZstr = XYZstr{1};
0088 
0089 %-Get parameter estimates, ResidualMS, (compute) fitted data & residuals
0090 %=======================================================================
0091 
0092 %-Parameter estimates: beta = xX.pKX*xX.K*y;
0093 %-----------------------------------------------------------------------
0094 beta  = mRes.betas(:, rno);
0095 
0096 %-Compute residuals
0097 %-----------------------------------------------------------------------
0098 R   = spm_sp('r',xX.xKXs,pr_spm_filter('apply',xX.K,y));
0099 
0100 %-Residual mean square: ResidualMS = sum(R.^2)/xX.trRV;
0101 %-----------------------------------------------------------------------
0102 ResidualMS = mRes.ResidualMS(rno);
0103 Bcov  = xX.Bcov;
0104 SE    = sqrt(ResidualMS*diag(Bcov));
0105 COL   = ['r','b','g','c','y','m','r','b','g','c','y','m'];
0106 
0107 
0108 %-Plot
0109 %=======================================================================
0110 
0111 % find out what to plot
0112 %-----------------------------------------------------------------------
0113 Cplot = {    'Contrast of parameter estimates',...
0114          'Fitted and adjusted responses',...
0115          'Event/epoch-related responses',...
0116          'Plots of parametric responses',...
0117         'Volterra Kernels'};
0118 
0119 
0120 % ensure options are appropriate
0121 %-----------------------------------------------------------------------
0122 if ~isfield(mRes,'Sess')
0123 
0124     Cplot = Cplot(1:2);
0125 else
0126     Sess  = mRes.Sess;
0127 end
0128 Cp     = spm_input('Plot',-1,'m',Cplot);
0129 Cplot  = Cplot{Cp};
0130 
0131 switch Cplot
0132 
0133 % select contrast to plot and compute fitted and adjusted data
0134 %----------------------------------------------------------------------
0135 case {'Contrast of parameter estimates','Fitted and adjusted responses'}
0136 
0137     % determine current contrasts
0138     %---------------------------------------------------------------
0139     if isempty(Ic)
0140       % determine which contrast
0141       %---------------------------------------------------------------
0142       [Ic marsD changef] = ui_get_contrasts(...
0143           marsD, 'T|F',1, 'Select contrast...', ' for plot', 1);
0144       if changef, xCon = get_contrasts(marsD); end
0145     end
0146     TITLE = {Cplot xCon(Ic).name};
0147 
0148     % fitted (corrected) data (Y = X1o*beta)
0149     %---------------------------------------------------------------
0150     Y     = spm_FcUtil('Yc',xCon(Ic),xX.xKXs,beta);
0151 
0152     % adjusted data
0153     %---------------------------------------------------------------
0154     y     = Y + R;
0155 
0156 end
0157 
0158 switch Cplot
0159 
0160 % plot parameter estimates
0161 %----------------------------------------------------------------------
0162 case 'Contrast of parameter estimates'
0163 
0164     % comute contrast of parameter estimates and standard error
0165     %--------------------------------------------------------------
0166     c     = xCon(Ic).c;
0167     cbeta = c'*beta;
0168     SE    = sqrt(ResidualMS*diag(c'*Bcov*c));
0169 
0170     % bar chart
0171     %--------------------------------------------------------------
0172     figure(Fgraph)
0173     subplot(2,1,2)
0174     h     = bar(cbeta);
0175     set(h,'FaceColor',[1 1 1]*.8)
0176     for j = 1:length(cbeta)
0177         line([j j],([SE(j) 0 - SE(j)] + cbeta(j)),...
0178                 'LineWidth',3,'Color','r')
0179     end
0180     set(gca,'XLim',[0.4 (length(cbeta) + 0.6)])
0181     XLAB  = 'effect';
0182     YLAB  = ['size of effect',XYZstr];
0183 
0184 
0185 % all fitted effects or selected effects
0186 %-----------------------------------------------------------------------
0187 case 'Fitted and adjusted responses'
0188 
0189     % get ordinates
0190     %---------------------------------------------------------------
0191     Xplot = {    'an explanatory variable',...
0192             'scan or time',...
0193             'a user specified ordinate'};
0194 
0195     Cx    = spm_input('plot against','!+1','m',Xplot);
0196 
0197     if     Cx == 1
0198 
0199         str  = 'Which column or effect?';
0200         i    = spm_input(str,'!+1','m',xX.Xnames);
0201         x    = xX.xKXs.X(:,i);
0202         XLAB = xX.Xnames{i};
0203 
0204     elseif Cx == 2
0205 
0206         if isfield(xX,'RT') & ~isempty(xX.RT)
0207             x    = xX.RT*[1:size(Y,1)]';
0208             XLAB = 'time {seconds}';
0209         else
0210             x    = [1:size(Y,1)]';
0211             XLAB = 'scan number';
0212         end
0213 
0214     elseif Cx == 3
0215 
0216         x    = spm_input('enter ordinate','!+1','e','',size(Y,1));
0217         XLAB = 'ordinate';
0218 
0219     end
0220 
0221     % plot
0222     %---------------------------------------------------------------
0223     figure(Fgraph)
0224     subplot(2,1,2)
0225     [p q] = sort(x);
0226     if all(diff(x(q)))
0227         plot(x(q),y(q),':b'); hold on
0228         plot(x(q),y(q),'.b','MarkerSize',8); hold on
0229         plot(x(q),Y(q),'r' ); hold off
0230 
0231     else
0232         plot(x(q),y(q),'.b','MarkerSize',8); hold on
0233         plot(x(q),Y(q),'.r','MarkerSize',16); hold off
0234         xlim = get(gca,'XLim');
0235         xlim = [-1 1]*diff(xlim)/4 + xlim;
0236         set(gca,'XLim',xlim)
0237 
0238     end
0239     YLAB  = ['response',XYZstr];
0240 
0241 
0242 % modeling evoked responses based on Sess
0243 %----------------------------------------------------------------------
0244 case 'Event/epoch-related responses'
0245 
0246 
0247     % average over sessions?
0248     %--------------------------------------------------------------
0249     ss    = length(Sess);
0250     if ss > 1
0251 
0252         % determine if the same basis functions have been used
0253         %------------------------------------------------------
0254         for s = 1:ss
0255             rep     = length(Sess{s}.bf) == length(Sess{1}.bf);
0256             if ~rep, break, end
0257             for t   = 1:length(Sess{s}.bf)
0258             rep = all(size(Sess{s}.bf{t}) == size(Sess{1}.bf{t}));
0259             if ~rep, break, end
0260             rep = all(all( Sess{s}.bf{t}  == Sess{1}.bf{t} ));
0261             if ~rep, break, end
0262             end
0263             if ~rep, break, end
0264         end
0265 
0266         % average over sessions?
0267         %------------------------------------------------------
0268         if rep
0269             str   = 'average over sessions?';
0270             rep   = spm_input(str,'+1','y/n',[1 0]);
0271         end
0272 
0273         % selected sessions
0274         %------------------------------------------------------
0275         if rep
0276             ss    = 1:ss;
0277         else
0278             str   = sprintf('which session (1 to %d)',ss);
0279             ss    = spm_input(str,'+1','n','1',1,ss);
0280         end
0281     end
0282 
0283     % get plot type
0284     %--------------------------------------------------------------
0285     Rplot = {    'fitted response',...
0286             'fitted response and PSTH',...
0287             'fitted response +/- standard error',...
0288             'fitted response and adjusted data'};
0289 
0290     if isempty(y), Rplot = Rplot([1 3]); end
0291     Cr      = spm_input('plot in terms of','+1','m',Rplot);
0292     TITLE   = Rplot{Cr};
0293     YLAB    = ['response',XYZstr];
0294     XLAB{1} = 'peri-stimulus time {secs}';
0295 
0296 
0297     % get selected trials
0298     %--------------------------------------------------------------
0299     tr    = length(Sess{ss(1)}.pst);
0300     if tr > 1
0301         str   = sprintf('which trials or conditions (1 to %d)',tr);
0302         tr    = spm_input(str,'+1','n');
0303     end
0304 
0305 
0306     % plot
0307     %--------------------------------------------------------------
0308     switch TITLE
0309     case 'fitted response and PSTH'
0310             str = 'bin size for PSTH {secs}';
0311             BIN = spm_input(str,'+1','r','2',1);
0312 
0313     otherwise,    BIN = 2;   end
0314 
0315     % reconstruct response with filtering
0316     %--------------------------------------------------------------
0317     figure(Fgraph)
0318     subplot(2,1,2)
0319     hold on
0320     dx    = xX.dt;
0321     XLim  = 0;
0322     u     = 1;
0323     for t = tr
0324 
0325         for s = ss
0326 
0327             % basis functions, filter and parameters
0328             %----------------------------------------------
0329             j    = 1:size(Sess{s}.sf{t},2):length(Sess{s}.ind{t});
0330             j    = Sess{s}.col(Sess{s}.ind{t}(j));
0331             B    = beta(j);
0332             X    = Sess{s}.bf{t};
0333             q    = 1:size(X,1);
0334             x    = (q - 1)*dx;
0335             K{1} = struct(  'HChoice',    'none',...
0336                      'HParam',    [],...
0337                     'LChoice',    xX.K{s}.LChoice,...
0338                     'LParam',    xX.K{s}.LParam,...
0339                     'row',        q,...
0340                     'RT',        dx);
0341 
0342             % fitted responses with standard error
0343             %----------------------------------------------
0344             KX       = pr_spm_filter('apply',K,X);
0345             Y(q,s)   = KX*B;
0346             se(:,s)  = sqrt(diag(X*Bcov(j,j)*X')*ResidualMS);
0347         end
0348 
0349         % average over sessions
0350         %------------------------------------------------------
0351         Y     = sum(Y,2)/length(ss);
0352 
0353         % peristimulus times and adjusted data (Y + R)
0354         %------------------------------------------------------
0355         pst   = [];
0356         y     = [];
0357         for s = ss
0358             p       = Sess{s}.pst{t}(:);
0359             bin     = round(p/dx);
0360             q       = find((bin >= 0) & (bin < size(X,1)));
0361             pst     = [pst; p];
0362             p       = R(Sess{s}.row(:));
0363             p(q)    = p(q) + Y(bin(q) + 1);
0364             y       = [y; p];
0365         end
0366 
0367 
0368         % PSTH
0369         %------------------------------------------------------
0370         INT    = -BIN:BIN:max(pst);
0371         PSTH   = [];
0372         SEM    = [];
0373         PST    = [];
0374         for k  = 1:(length(INT) - 1)
0375             q = find(pst > INT(k) & pst <= INT(k + 1));
0376             n = length(q);
0377             if n
0378                 PSTH = [PSTH mean(y(q))];
0379                 SEM  = [SEM std(y(q))/sqrt(n)];
0380                 PST  = [PST mean(pst(q))];
0381             end
0382         end
0383 
0384         % plot
0385         %------------------------------------------------------
0386         switch TITLE
0387 
0388             case 'fitted response'
0389             %----------------------------------------------
0390             plot(x,Y,COL(u))
0391 
0392             case 'fitted response and PSTH'
0393             %----------------------------------------------
0394             errorbar(PST,PSTH,SEM,[':' COL(u)])
0395             plot(PST,PSTH,['.' COL(u)],'MarkerSize',16)
0396             plot(PST,PSTH,COL(u),'LineWidth',2)
0397             plot(x,Y,['-.' COL(u)])
0398 
0399             case 'fitted response +/- standard error'
0400             %----------------------------------------------
0401             plot(x,Y,COL(u))
0402             plot(x,Y + se,['-.' COL(u)],x,Y - se,['-.' COL(u)])
0403 
0404             case 'fitted response and adjusted data'
0405             %----------------------------------------------
0406             plot(x,Y,COL(u),pst,y,['.' COL(u)],'MarkerSize',8)
0407 
0408         end
0409 
0410         % xlabel
0411         %------------------------------------------------------
0412         XLAB{end + 1} = [Sess{s}.name{t} ' - ' COL(u)];
0413         u    = u + 1;
0414         XLim = max([XLim max(x)]);
0415 
0416     end
0417 
0418     hold off; axis on
0419     set(gca,'XLim',[-4 XLim])
0420 
0421 % modeling evoked responses based on Sess
0422 %----------------------------------------------------------------------
0423 case 'Plots of parametric responses'
0424 
0425     % Get session
0426     %--------------------------------------------------------------
0427     s     = length(Sess);
0428     if  s > 1
0429         s = spm_input('which session','+1','n1',[],s);
0430     end
0431 
0432     % Get [parametric] trial
0433     %--------------------------------------------------------------
0434     Vname = {};
0435     j     = [];
0436     for i = 1:length(Sess{s}.Pv)
0437         if length(Sess{s}.Pv{i})
0438             Vname{end + 1} = Sess{s}.name{i};
0439             j              = [j i];
0440         end
0441     end
0442     if isempty(Vname)
0443       warning('No parametric responses specified');
0444     end
0445     t     = j(spm_input('which effect','+1','m',Vname));
0446 
0447     % parameter estimates and fitted response
0448     %-------------------------------------------------------------
0449     B     = beta(Sess{s}.col(Sess{s}.ind{t}));
0450     Q     = Sess{s}.Pv{t};
0451     SF    = Sess{s}.sf{t};
0452     SF    = SF(find(SF(:,1)),:);
0453     X     = Sess{s}.bf{t};
0454     q     = 1:size(X,1);
0455     x     = q*xX.dt;
0456     K{1}  = struct(        'HChoice',    'none',...
0457                 'HParam',    [],...
0458                 'LChoice',    xX.K{s}.LChoice,...
0459                 'LParam',    xX.K{s}.LParam,...
0460                 'row',        q,...
0461                 'RT',        xX.dt);
0462 
0463     KX    = pr_spm_filter('apply',K,X);
0464     p     = size(SF,2);
0465     b     = [];
0466     for i = 1:size(KX,2)
0467         b = [b SF*B([1:p] + (i - 1)*p)];
0468     end
0469     Y     = KX*b';
0470 
0471     % plot
0472     %-------------------------------------------------------------
0473     figure(Fgraph)
0474     subplot(2,1,2)
0475     surf(x',Q',Y')
0476     shading flat
0477     TITLE = Sess{s}.name{t};
0478     XLAB  = 'perstimulus time (secs)';
0479     YLAB  = Sess{s}.Pname{t};
0480     zlabel(['responses',XYZstr]);
0481 
0482 
0483 % modeling evoked responses based on Sess
0484 %----------------------------------------------------------------------
0485 case 'Volterra Kernels'
0486 
0487 
0488     % Get session
0489     %--------------------------------------------------------------
0490     s     = length(Sess);
0491     if  s > 1
0492         s = spm_input('which session','+1','n1',[],s);
0493     end
0494 
0495     % Get [non-parametric] trial
0496     %--------------------------------------------------------------
0497     Vname = {};
0498     j     = [];
0499     for i = 1:length(Sess{s}.name)
0500         Vname{end + 1} = Sess{s}.name{i};
0501     end
0502     t     = spm_input('which effect','+1','m',Vname);
0503 
0504     % Parameter estimates
0505     %--------------------------------------------------------------
0506     B     = beta(Sess{s}.col(Sess{s}.ind{t}));
0507 
0508     % plot
0509     %--------------------------------------------------------------
0510     figure(Fgraph)
0511     subplot(2,1,2)
0512 
0513     % second order kernel
0514     %--------------------------------------------------------------
0515     if iscell(Sess{s}.bf{t})
0516 
0517         Y     = 0;
0518         for i = 1:length(Sess{s}.bf{t})
0519             Y = Y + B(i)*Sess{s}.bf{t}{i};
0520         end
0521         p     = ([1:size(Y,2)] - 1)*xX.dt;
0522         q     = ([1:size(Y,1)] - 1)*xX.dt;
0523         imagesc(p,q,Y)
0524         axis xy
0525 
0526         TITLE = {'Second order Volterra Kernel' Sess{s}.name{t}};
0527         XLAB  = 'perstimulus time (secs)';
0528         YLAB  = 'perstimulus time (secs)';
0529 
0530     % first  order kernel
0531     %--------------------------------------------------------------
0532     else
0533         j     = 1:size(Sess{s}.sf{t},2):length(Sess{s}.ind{t});
0534         Y     = Sess{s}.bf{t}*B(j);
0535         p     = ([1:length(Y)] - 1)*xX.dt;
0536         plot(p,Y)
0537         grid on
0538 
0539         TITLE = {'First order Volterra Kernel' Sess{s}.name{t}};
0540         XLAB  = 'perstimulus time (secs)';
0541         YLAB  = ['responses',XYZstr];
0542 
0543     end
0544 end
0545 
0546 %-Label and call Plot UI
0547 %----------------------------------------------------------------------
0548 axis square
0549 if strcmp(get(get(gca,'Children'),'type'),'image')
0550     axis image
0551 end
0552 xlabel(XLAB,'FontSize',10)
0553 ylabel(YLAB,'FontSize',10)
0554 title(TITLE,'FontSize',16)
0555 
0556 spm_results_ui('PlotUi',gca)
0557 
0558 % Complete return values
0559 r_st = mars_struct('fillafromb', def_r_st, struct(...
0560     'Y', Y, 'y', y, 'beta', beta, 'SE', SE, 'cbeta', cbeta));
Generated on Wed 11-May-2022 16:26:09 by m2html © 2003-2019