Home > marsbar > @mardo_2 > 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  - MarsBaR 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
          Bcov   - covariance of parameter estimates
          cbeta  = betas multiplied by contrast
 marsD   - design structure, with possibly added contrasts
 changef - set to 1 if design has changed

 see spm2 version of spm_graph for details
_______________________________________________________________________
 @(#)spm_graph.m    2.43 Karl Friston 03/12/19

 $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  - MarsBaR 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 %          Bcov   - covariance 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 spm2 version of spm_graph for details
0031 %_______________________________________________________________________
0032 % @(#)spm_graph.m    2.43 Karl Friston 03/12/19
0033 %
0034 % $Id$
0035 
0036 if ~is_mars_estimated(marsD)
0037   error('Need estimated design for plot');
0038 end
0039 if nargin < 2
0040   rno = [];
0041 end
0042 if nargin < 3
0043   Ic = [];
0044 end
0045 changef = 0;
0046 
0047 % make values ready for return
0048 def_r_st = struct(...
0049     'Y', [],...
0050     'y', [],...
0051     'beta', [],...
0052     'SE', [],...
0053     'cbeta',[], ...
0054     'PSTH',[]);
0055 cbeta = []; PSTH = [];
0056 
0057 % get stuff from object
0058 SPM = des_struct(marsD);
0059 xCon = SPM.xCon;
0060   
0061 % Check if we want to, and can, assume region no is 1
0062 if isempty(rno) 
0063   if n_regions(SPM.marsY) > 1
0064     error('Need to specify region number');
0065   end
0066   rno = 1;
0067 end
0068 
0069 % Get required data and filter
0070 sY    = summary_data(SPM.marsY);
0071 y     = apply_filter(marsD, sY(:, rno));
0072 
0073 % Get design matrix for convenience
0074 xX = SPM.xX;
0075 
0076 % Label for region
0077 XYZstr = region_name(SPM.marsY, rno);
0078 XYZstr = XYZstr{1};
0079 
0080 %-Get parameter estimates, ResidualMS, (compute) fitted data & residuals
0081 %=======================================================================
0082 
0083 %-Parameter estimates: beta = xX.pKX*xX.K*y;
0084 %-----------------------------------------------------------------------
0085 beta  = SPM.betas(:, rno);
0086 
0087 %-Residual mean square: ResidualMS = sum(R.^2)/xX.trRV;
0088 %-----------------------------------------------------------------------
0089 ResidualMS = SPM.ResidualMS(rno);
0090 Bcov  = ResidualMS*SPM.xX.Bcov;
0091 
0092 %-Get Graphics figure handle
0093 %-----------------------------------------------------------------------
0094 Fgraph = spm_figure('GetWin','Graphics');
0095 
0096 
0097 %-Delete previous axis and their pagination controls (if any)
0098 %-----------------------------------------------------------------------
0099 spm_results_ui('Clear',Fgraph,2);
0100 
0101 %-Compute residuals
0102 %-----------------------------------------------------------------------
0103 if isempty(y)
0104 
0105     % make R = NaN so it will not be plotted
0106     %---------------------------------------------------------------
0107     R   = NaN*ones(size(SPM.xX.X,1),1);
0108 
0109 else
0110     % residuals (non-whitened)
0111     %---------------------------------------------------------------
0112     R   = spm_sp('r',SPM.xX.xKXs,y);
0113 
0114 end
0115 
0116 %-Get parameter and hyperparameter estimates
0117 %=======================================================================
0118 %-Parameter estimates:   beta = xX.pKX*xX.K*y;
0119 %-Residual mean square: ResidualMS = sum(R.^2)/xX.trRV
0120 %---------------------------------------------------------------
0121 
0122 P05_Z = 1.6449;                    % = spm_invNcdf(1 - 0.05);
0123 CI    = P05_Z;
0124 
0125 %-Colour specifications and index;
0126 %-----------------------------------------------------------------------
0127 Col   = [0 0 0; .8 .8 .8; 1 .5 .5];
0128 
0129 %-Plot
0130 %=======================================================================
0131 
0132 % find out what to plot
0133 %-----------------------------------------------------------------------
0134 Cplot = {    'Contrast estimates and 90% C.I.',...
0135          'Fitted responses',...
0136          'Event-related responses',...
0137          'Parametric responses',...
0138         'Volterra Kernels'};
0139 
0140 
0141 % ensure options are appropriate
0142 %-----------------------------------------------------------------------
0143 try
0144     Sess  = SPM.Sess;
0145 catch
0146     Cplot = Cplot(1:2);    
0147 end
0148 Cplot  = Cplot{spm_input('Plot',-1,'m',Cplot)};
0149 
0150 switch Cplot
0151 
0152 % select contrast if
0153 %----------------------------------------------------------------------
0154 case {'Contrast estimates and 90% C.I.','Fitted responses'}
0155 
0156         if isempty(Ic)
0157       % determine which contrast
0158       %---------------------------------------------------------------
0159       [Ic marsD changef] = ui_get_contrasts(...
0160           marsD, 'T|F',1, 'Select contrast...', ' for plot', 1);
0161       if changef, SPM.xCon = get_contrasts(marsD); end
0162     end
0163     TITLE = {Cplot SPM.xCon(Ic).name};
0164 
0165 % select session and trial if
0166 %----------------------------------------------------------------------
0167 case {'Event-related responses','Parametric responses','Volterra Kernels'}
0168 
0169     % get session
0170     %--------------------------------------------------------------
0171     s     = length(Sess);
0172     if  s > 1
0173         s = spm_input('which session','+1','n1',1,s);
0174     end
0175 
0176     % effect names
0177     %--------------------------------------------------------------
0178     switch Cplot
0179     case 'Volterra Kernels'
0180         u = length(Sess(s).Fc);
0181     otherwise
0182         u = length(Sess(s).U);
0183     end
0184     Uname = {};
0185     for i = 1:u
0186         Uname{i} = Sess(s).Fc(i).name;
0187     end
0188 
0189     % get effect
0190     %--------------------------------------------------------------
0191     str   = sprintf('which effect');
0192     u     = spm_input(str,'+1','m',Uname);
0193 
0194     % bin size
0195     %--------------------------------------------------------------
0196     dt    = SPM.xBF.dt;
0197 
0198 end
0199 
0200 switch Cplot
0201 
0202 % plot parameter estimates
0203 %----------------------------------------------------------------------
0204 case 'Contrast estimates and 90% C.I.'
0205 
0206     % compute contrast of parameter estimates and 90% C.I.
0207     %--------------------------------------------------------------
0208     cbeta = SPM.xCon(Ic).c'*beta;
0209     CI    = CI*sqrt(diag(SPM.xCon(Ic).c'*Bcov*SPM.xCon(Ic).c));
0210 
0211     % bar chart
0212     %--------------------------------------------------------------
0213     figure(Fgraph)
0214     subplot(2,1,2)
0215     cla
0216     hold on
0217 
0218     % estimates
0219     %--------------------------------------------------------------
0220     h     = bar(cbeta);
0221     set(h,'FaceColor',Col(2,:))
0222 
0223     % standard error
0224     %--------------------------------------------------------------
0225     for j = 1:length(cbeta)
0226         line([j j],([CI(j) 0 - CI(j)] + cbeta(j)),...
0227                 'LineWidth',6,'Color',Col(3,:))
0228     end
0229 
0230     title(TITLE,'FontSize',12)
0231     xlabel('contrast')
0232     ylabel(['contrast estimate',XYZstr])
0233     set(gca,'XLim',[0.4 (length(cbeta) + 0.6)])
0234     hold off
0235 
0236     % set Y to empty so outputs are assigned
0237     %-------------------------------------------------------------
0238     Y = [];
0239 
0240 % all fitted effects or selected effects
0241 %-----------------------------------------------------------------------
0242 case 'Fitted responses'
0243 
0244     % predicted or adjusted response
0245     %---------------------------------------------------------------
0246     str   = 'predicted or adjusted response?';
0247     if spm_input(str,'!+1','b',{'predicted','adjusted'},[1 0]);
0248 
0249         % fitted (predicted) data (Y = X1*beta)
0250         %--------------------------------------------------------
0251         Y = SPM.xX.X*SPM.xCon(Ic).c*pinv(SPM.xCon(Ic).c)*beta;
0252     else
0253 
0254         % fitted (corrected)  data (Y = X1o*beta)
0255         %-------------------------------------------------------
0256         Y = spm_FcUtil('Yc',SPM.xCon(Ic),SPM.xX.xKXs,beta);
0257 
0258     end
0259 
0260     % adjusted data
0261     %---------------------------------------------------------------
0262     y     = Y + R;
0263 
0264     % get ordinates
0265     %---------------------------------------------------------------
0266     Xplot = {    'an explanatory variable',...
0267             'scan or time',...
0268             'a user specified ordinate'};
0269     Cx    = spm_input('plot against','!+1','m',Xplot);
0270 
0271     % an explanatory variable
0272     %---------------------------------------------------------------
0273     if     Cx == 1
0274 
0275         str  = 'Which explanatory variable?';
0276         i    = spm_input(str,'!+1','m',SPM.xX.name);
0277         x    = SPM.xX.xKXs.X(:,i);
0278         XLAB = SPM.xX.name{i};
0279 
0280     % scan or time
0281     %---------------------------------------------------------------
0282     elseif Cx == 2
0283 
0284         if isfield(SPM.xY,'RT')
0285             x    = SPM.xY.RT*[1:size(Y,1)]';
0286             XLAB = 'time {seconds}';
0287         else
0288             x    = [1:size(Y,1)]';
0289             XLAB = 'scan number';
0290         end
0291 
0292     % user specified
0293     %---------------------------------------------------------------
0294     elseif Cx == 3
0295 
0296         x    = spm_input('enter ordinate','!+1','e','',size(Y,1));
0297         XLAB = 'ordinate';
0298 
0299     end
0300 
0301     % plot
0302     %---------------------------------------------------------------
0303     figure(Fgraph)
0304     subplot(2,1,2)
0305     cla
0306     hold on
0307     [p q] = sort(x);
0308     if all(diff(x(q)))
0309         plot(x(q),Y(q),'LineWidth',4,'Color',Col(2,:));
0310         plot(x(q),y(q),':','Color',Col(1,:));
0311         plot(x(q),y(q),'.','MarkerSize',8, 'Color',Col(3,:)); 
0312 
0313     else
0314         plot(x(q),Y(q),'.','MarkerSize',16,'Color',Col(1,:));
0315         plot(x(q),y(q),'.','MarkerSize',8, 'Color',Col(2,:));
0316         xlim = get(gca,'XLim');
0317         xlim = [-1 1]*diff(xlim)/4 + xlim;
0318         set(gca,'XLim',xlim)
0319 
0320     end
0321     title(TITLE,'FontSize',12)
0322     xlabel(XLAB)
0323     ylabel(['response',XYZstr])
0324     legend('fitted','plus error')
0325     hold off
0326 
0327 % modeling evoked responses based on Sess
0328 %----------------------------------------------------------------------
0329 case 'Event-related responses'
0330 
0331     % get plot type
0332     %--------------------------------------------------------------
0333     Rplot   = {    'fitted response and PSTH',...
0334             'fitted response and 90% C.I.',...
0335             'fitted response and adjusted data'};
0336 
0337     if isempty(y)
0338         TITLE = Rplot{2};
0339     else
0340         TITLE = Rplot{spm_input('plot in terms of','+1','m',Rplot)};
0341     end
0342 
0343     % plot
0344     %--------------------------------------------------------------
0345     switch TITLE
0346     case 'fitted response and PSTH'
0347 
0348 
0349         % build a simple FIR model subpartition (X); bin size = TR
0350         %------------------------------------------------------
0351         str         = 'bin size (secs)';
0352                 BIN         = sprintf('%0.2f',SPM.xY.RT);
0353         BIN         = spm_input(str,'!+1','r',BIN);
0354         xBF         = SPM.xBF;
0355         U           = Sess(s).U(u);
0356         U.u         = U.u(:,1);
0357         xBF.name    = 'Finite Impulse Response';
0358         xBF.order   = round(32/BIN);
0359         xBF.length  = xBF.order*BIN;
0360         xBF         = pr_spm_get_bf(xBF);
0361         BIN         = xBF.length/xBF.order;
0362         X           = pr_spm_volterra(U,xBF.bf,1);
0363         k           = SPM.nscan(s);
0364         X           = X([0:(k - 1)]*SPM.xBF.T + SPM.xBF.T0 + 32,:);
0365 
0366         % place X in SPM.xX.X
0367         %------------------------------------------------------
0368         jX          = Sess(s).row;
0369         iX          = Sess(s).col(Sess(s).Fc(u).i);
0370         iX0         = [1:size(SPM.xX.X,2)];
0371         iX0(iX)     = [];
0372         X           = [X SPM.xX.X(jX,iX0)];
0373         X           = SPM.xX.W(jX,jX)*X;
0374         X           = [X SPM.xX.K(s).X0];
0375 
0376         % Re-estimate to get PSTH and CI
0377         %------------------------------------------------------
0378         j           = xBF.order;
0379         xX          = spm_sp('Set',X);
0380         pX          = spm_sp('x-',xX);
0381         PSTH        = pX*y(jX);
0382         res         = spm_sp('r',xX,y(jX));
0383         df          = size(X,1) - size(X,2);
0384         bcov        = pX*pX'*sum(res.^2)/df;
0385         PSTH        = PSTH(1:j)/dt;
0386         PST         = [1:j]*BIN - BIN/2;
0387         PCI         = CI*sqrt(diag(bcov(1:j,(1:j))))/dt;
0388     end
0389 
0390     % basis functions and parameters
0391     %--------------------------------------------------------------
0392     X     = SPM.xBF.bf/dt;
0393     x     = ([1:size(X,1)] - 1)*dt;
0394     j     = Sess(s).col(Sess(s).Fc(u).i(1:size(X,2)));
0395     B     = beta(j);
0396     
0397     % fitted responses with standard error
0398     %--------------------------------------------------------------
0399     Y     = X*B;
0400     CI    = CI*sqrt(diag(X*Bcov(j,j)*X'));
0401 
0402     % peristimulus times and adjusted data (y = Y + R)
0403     %--------------------------------------------------------------
0404     pst   = Sess(s).U(u).pst;
0405     bin   = round(pst/dt);
0406     q     = find((bin >= 0) & (bin < size(X,1)));
0407     y     = R(Sess(s).row(:));
0408     pst   = pst(q);
0409     y     = y(q) + Y(bin(q) + 1);
0410 
0411     % plot
0412     %--------------------------------------------------------------
0413     figure(Fgraph)
0414     subplot(2,1,2)
0415     hold on
0416     switch TITLE
0417 
0418         case 'fitted response and PSTH'
0419         %------------------------------------------------------
0420         errorbar(PST,PSTH,PCI)
0421         plot(PST,PSTH,'LineWidth',4,'Color',Col(2,:))
0422         plot(x,Y,'-.','Color',Col(3,:))
0423 
0424         case 'fitted response and 90% C.I.'
0425         %------------------------------------------------------
0426         plot(x,Y,'Color',Col(2,:),'LineWidth',4)
0427         plot(x,Y + CI,'-.',x,Y - CI,'-.','Color',Col(1,:))
0428 
0429         case 'fitted response and adjusted data'
0430         %------------------------------------------------------
0431         plot(x,Y,'Color',Col(2,:),'LineWidth',4)
0432         plot(pst,y,'.','Color',Col(3,:))
0433 
0434     end
0435 
0436     % label
0437     %-------------------------------------------------------------
0438     [i j] = max(Y);
0439     text(ceil(1.1*x(j)),i,Sess(s).Fc(u).name,'FontSize',8);
0440     title(TITLE,'FontSize',12)
0441     xlabel('peristimulus time {secs}')
0442     ylabel(['response',XYZstr])
0443     hold off
0444 
0445 
0446 % modeling evoked responses based on Sess
0447 %----------------------------------------------------------------------
0448 case 'Parametric responses'
0449 
0450 
0451     % return gracefully if no parameters
0452     %--------------------------------------------------------------
0453     if ~Sess(s).U(u).P(1).h, return, end
0454 
0455     % basis functions
0456     %--------------------------------------------------------------
0457     bf    = SPM.xBF.bf;
0458     pst   = ([1:size(bf,1)] - 1)*dt;
0459 
0460     % orthogonalised expansion of parameteric variable
0461     %--------------------------------------------------------------
0462     str   = 'which parameter';
0463     p     = spm_input(str,'+1','m',{Sess(s).U(u).P.name});
0464     P     = Sess(s).U(u).P(p).P;
0465     q     = [];
0466     for i = 0:Sess(s).U(u).P(p).h;
0467         q = [q spm_en(P).^i];
0468     end
0469     q     = spm_orth(q);
0470 
0471 
0472     % parameter estimates for this effect
0473     %--------------------------------------------------------------
0474     j     = Sess(s).col(Sess(s).Fc(u).i);
0475     B     = beta(j);
0476 
0477     % reconstruct trial-specific responses
0478     %--------------------------------------------------------------
0479     Y     = zeros(size(bf,1),size(q,1));
0480     uj    = Sess(s).U(u).P(p).i;
0481     for i = 1:size(P,1)
0482         U      = sparse(1,uj,q(i,:),1,size(Sess(s).U(u).u,2));
0483         X      = kron(U,bf);
0484         Y(:,i) = X*B;
0485     end
0486     [P j] = sort(P);
0487     Y     = Y(:,j);
0488 
0489     % plot
0490     %--------------------------------------------------------------
0491     figure(Fgraph)
0492     subplot(2,2,3)
0493     surf(pst,P,Y')
0494     shading flat
0495     title(Sess(s).U(u).name{1},'FontSize',12)
0496     xlabel('PST {secs}')
0497     ylabel(Sess(s).U(u).P(p).name)
0498     zlabel(['responses',XYZstr])
0499     axis square
0500 
0501     % plot
0502     %--------------------------------------------------------------
0503     subplot(2,2,4)
0504     [j i] = max(mean(Y,2));
0505     plot(P,Y(i,:),'LineWidth',4,'Color',Col(2,:))
0506     str   = sprintf('response at %0.1fs',i*dt);
0507     title(str,'FontSize',12)
0508     xlabel(Sess(s).U(u).P(p).name)
0509     axis square
0510     grid on
0511 
0512 
0513 % modeling evoked responses based on Sess
0514 %----------------------------------------------------------------------
0515 case 'Volterra Kernels'
0516 
0517     % Parameter estimates and basis functions
0518     %------------------------------------------------------
0519     bf    = SPM.xBF.bf/dt;
0520     pst   = ([1:size(bf,1)] - 1)*dt;
0521 
0522     % second order kernel
0523     %--------------------------------------------------------------
0524     if u > length(Sess(s).U)
0525 
0526         % Parameter estimates and kernel
0527         %------------------------------------------------------
0528         j     = Sess(s).col(Sess(s).Fc(u).i);
0529         B     = beta(j);
0530         i     = 1;
0531         Y     = 0;
0532         for p = 1:size(bf,2)
0533         for q = 1:size(bf,2)
0534                  Y = Y + B(i)*bf(:,p)*bf(:,q)';
0535             i = i + 1;
0536         end
0537         end
0538 
0539         % plot
0540         %------------------------------------------------------
0541         figure(Fgraph)
0542         subplot(2,2,3)
0543         imagesc(pst,pst,Y)
0544         axis xy
0545         axis image
0546 
0547         title('2nd order Kernel','FontSize',12);
0548         xlabel('perstimulus time {secs}')
0549         ylabel('perstimulus time {secs}')
0550 
0551         subplot(2,2,4)
0552         plot(pst,Y)
0553         axis square
0554         grid on
0555 
0556         title(Sess(s).Fc(u).name,'FontSize',12);
0557         xlabel('perstimulus time {secs}')
0558 
0559 
0560     % first  order kernel
0561     %--------------------------------------------------------------
0562     else
0563         j     = Sess(s).col(Sess(s).Fc(u).i(1:size(bf,2)));
0564         B     = beta(j);
0565         Y     = bf*B;
0566 
0567         % plot
0568         %------------------------------------------------------
0569         figure(Fgraph)
0570         subplot(2,1,2)
0571         plot(pst,Y)
0572         grid on
0573         axis square
0574 
0575         title({'1st order Volterra Kernel' Sess(s).Fc(u).name},...
0576             'FontSize',12);
0577         xlabel('perstimulus time {secs}')
0578         ylabel(['impluse response',XYZstr])
0579     end
0580 
0581 end
0582 
0583 
0584 %-call Plot UI
0585 %----------------------------------------------------------------------
0586 spm_results_ui('PlotUi',gca)
0587 
0588 % Complete return values
0589 r_st = mars_struct('fillafromb', def_r_st, struct(...
0590     'Y', Y, 'y', y, 'beta', beta, 'Bcov', Bcov, ...
0591     'cbeta', cbeta, 'PSTH', PSTH));
Generated on Wed 11-May-2022 16:26:09 by m2html © 2003-2019