eeg_anova_script

 eeg_anova_script - example task ERP analysis

0001 
0002 % eeg_anova_script - example task ERP analysis
0003 
0004 %  Analysis script for ERP components in
0005 %  fixed and variable task waveforms for PTSD-PET project
0006 
0007 % Define sample rate & epoch parameters
0008 fprintf('Defining ERP parameters\n');
0009 sample_rate = 2.5;
0010 epoch_start = -200;
0011 epoch_end = 1500;
0012 points = 680;
0013 times = linspace(epoch_start,(epoch_end - sample_rate),points);
0014 
0015 % Define component timing & rows
0016 comp = [145 195 250 385 550 750 950 1150 1350];
0017 comp_rows = (comp - epoch_start) / sample_rate;
0018 
0019 % Load data (680 rows x 124 cols)
0020 path = 'd:\data_emse\ptsdpet\link14hz\';
0021 for i=1:10,
0022     file = sprintf('%sc%02do_link14hz.txt',path,i);
0023     var  = sprintf('c%02do_link14hz',i);
0024     if ~exist(var,'var')
0025         fprintf('loading %s\n',file);
0026         load(file);
0027     end
0028     file = sprintf('%sc%02dt_link14hz.txt',path,i);
0029     var  = sprintf('c%02dt_link14hz',i);
0030     if ~exist(var,'var')
0031         fprintf('loading %s\n',file);
0032         load(file);
0033     end
0034     file = sprintf('%sp%02do_link14hz.txt',path,i);
0035     var  = sprintf('p%02do_link14hz',i);
0036     if ~exist(var,'var')
0037         fprintf('loading %s\n',file);
0038         load(file);
0039     end
0040     file = sprintf('%sp%02dt_link14hz.txt',path,i);
0041     var  = sprintf('p%02dt_link14hz',i);
0042     if ~exist(var,'var')
0043         fprintf('loading %s\n',file);
0044         load(file);
0045     end
0046 end
0047 
0048 % Select data rows for each component
0049 
0050 cont = zeros(1,2,124);
0051 ptsd = zeros(1,2,124);
0052 
0053 for i = 1:max(size(comp))
0054     c = comp(i);
0055     t = comp_rows(i);
0056     for j=1:10,
0057         sub = sprintf('c%02do_link14hz',j);
0058         fprintf('Component %d, selecting subject %s\n',c,sub);
0059         data = eval(sub);
0060         cont(j,1,:) = data(t,:);
0061 
0062         sub = sprintf('c%02dt_link14hz',j);
0063         fprintf('Component %d, selecting subject %s\n',c,sub);
0064         data = eval(sub);
0065         cont(j,2,:) = data(t,:);
0066         
0067         sub = sprintf('p%02do_link14hz',j);
0068         fprintf('Component %d, selecting subject %s\n',c,sub);
0069         data = eval(sub);
0070         ptsd(j,1,:) = data(t,:);
0071         
0072         sub = sprintf('p%02dt_link14hz',j);
0073         fprintf('Component %d, selecting subject %s\n',c,sub);
0074         data = eval(sub);
0075         ptsd(j,2,:) = data(t,:);
0076     end
0077     
0078     % Combine component data into single matrix
0079     data = [cont;ptsd]; % 20x2x124
0080     
0081     % Output ANOVA analysis.
0082     file = sprintf('%scomp%d_anova.txt',path,c);
0083     ANOVA = fopen(file,'w');
0084     file = sprintf('%scomp%d_p_task.txt',path,c);
0085     TASK = fopen(file,'w');
0086     file = sprintf('%scomp%d_p_group.txt',path,c);
0087     GROUP = fopen(file,'w');
0088     file = sprintf('%scomp%d_p_taskxgroup.txt',path,c);
0089     TASKxGROUP = fopen(file,'w');
0090     
0091     formatstr = '';
0092     
0093     fprintf('Doing ANOVA on electrode:\n');
0094     
0095     x = 0;
0096     for k=1:124,
0097         
0098         % Extract data for analysis
0099         anova_data = data(:,:,k);
0100         
0101         if (x==24),
0102             x=0;
0103             fprintf('%3d\n',k);
0104         else,
0105             x = x + 1;
0106             fprintf('%3d ',k);
0107         end
0108         fprintf(ANOVA,'\n\nELECTRODE %3d\n',k);
0109         
0110         % Calculate & output summary stats
0111         avg.cont = mean(anova_data( 1:10,:));
0112         avg.ptsd = mean(anova_data(10:20,:));
0113         sd.cont = std(anova_data( 1:10,:));
0114         sd.ptsd = std(anova_data(10:20,:));
0115         se.cont = sd.cont / sqrt(10);
0116         se.ptsd = sd.ptsd / sqrt(10);
0117         fprintf(ANOVA,'%-12s\t%-12s\t%12s\t%12s\t%12s\n','Group','Task','Mean','StDev','SE');
0118         fprintf(ANOVA,'%-12s\t%-12s\t%12s\t%12s\t%12s\n','-----','----','----','-----','--');
0119         fprintf(ANOVA,'%-12s\t%-12s\t%12.6f\t%12.6f\t%12.6f\n','Cont','Fixed',   avg.cont(1),sd.cont(1),se.cont(1));
0120         fprintf(ANOVA,'%-12s\t%-12s\t%12.6f\t%12.6f\t%12.6f\n','Cont','Variable',avg.cont(2),sd.cont(2),se.cont(2));
0121         fprintf(ANOVA,'%-12s\t%-12s\t%12.6f\t%12.6f\t%12.6f\n','PTSD','Fixed',   avg.ptsd(1),sd.ptsd(1),se.ptsd(1));
0122         fprintf(ANOVA,'%-12s\t%-12s\t%12.6f\t%12.6f\t%12.6f\n\n','PTSD','Variable',avg.ptsd(2),sd.ptsd(2),se.ptsd(2));
0123         
0124         % Run 2way ANOVA
0125         [p, tab] = anova2(anova_data,10,'off');
0126         
0127         % Collate the p values
0128         task(k) = p(1);
0129         group(k) = p(2);
0130         taskxgroup(k) = p(3);
0131         
0132         % Output ANOVA table
0133         tab(2,1) = {'Task'};
0134         tab(3,1) = {'Group'};
0135         tab(4,1) = {'Task x Group'};
0136         for m=1:6,
0137             for n=1:6,
0138                 if n==1, fprintf(ANOVA,'%-12s\t',tab{m,n});
0139                 else
0140                     if m==1, fprintf(ANOVA,'%12s\t',tab{m,n});
0141                     else     fprintf(ANOVA,'%12.6f\t',tab{m,n});
0142                     end
0143                 end
0144             end
0145             fprintf(ANOVA,'\n');
0146         end
0147         
0148         % This string is used to output all the p values below
0149         formatstr = strcat(formatstr,'%12.6f\t');
0150     
0151     end
0152     fclose(ANOVA);
0153     
0154     % Output the p values
0155     fprintf(TASK,formatstr,task); fclose(TASK);
0156     fprintf(GROUP,formatstr,group); fclose(GROUP);
0157     fprintf(TASKxGROUP,formatstr,taskxgroup); fclose(TASKxGROUP);
0158     
0159 end