RadialMaze

RadialMaze - Compute simple statistics for the 8-arm radial maze task.

 Compute the number of reference and working memory errors, rewards found, etc.

  USAGE

    [counts,indices] = RadialMaze(data,baitedArms,<options>)

    data           list of [rat day trial configuration group arm] tuples
    baitedArms     list of baited arms (from 1 to 8) for each (rat,configuration)
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'plots'       analyses to plot (default = {'individuals','groups'})
     'measures'    behavioral measures to compute/plot (default = {'perf'})
     'anova'       perform ANOVAs for selected measures (default = 'on')
    =========================================================================

  OUTPUT


    counts     a N x 12 matrix, where each line lists one case
                (rat,day,trial,config,group,measure1,...,measureN)
                    where the measures are nRefErrors, nWorkErrors, nRewardsFound,
              nVisits, performance, nPenaltyErrors, nRefErrors, and nAdjustedErrors
    indices     in the same format (N x 12 matrix) as counts


  NOTE

    The option 'plots' determines whether individual and/or group analyses should
    be plotted. The option 'measures' determines which kinds of analyses should be
    plotted: possible values include 'ref','work','rewards','visits','perf','penalty',
    and 'adjusted'.

0001 function [counts,indices] = RadialMaze(data,baitedArms,varargin)
0002 
0003 %RadialMaze - Compute simple statistics for the 8-arm radial maze task.
0004 %
0005 % Compute the number of reference and working memory errors, rewards found, etc.
0006 %
0007 %  USAGE
0008 %
0009 %    [counts,indices] = RadialMaze(data,baitedArms,<options>)
0010 %
0011 %    data           list of [rat day trial configuration group arm] tuples
0012 %    baitedArms     list of baited arms (from 1 to 8) for each (rat,configuration)
0013 %    <options>      optional list of property-value pairs (see table below)
0014 %
0015 %    =========================================================================
0016 %     Properties    Values
0017 %    -------------------------------------------------------------------------
0018 %     'plots'       analyses to plot (default = {'individuals','groups'})
0019 %     'measures'    behavioral measures to compute/plot (default = {'perf'})
0020 %     'anova'       perform ANOVAs for selected measures (default = 'on')
0021 %    =========================================================================
0022 %
0023 %  OUTPUT
0024 %
0025 %
0026 %    counts     a N x 12 matrix, where each line lists one case
0027 %                (rat,day,trial,config,group,measure1,...,measureN)
0028 %                    where the measures are nRefErrors, nWorkErrors, nRewardsFound,
0029 %              nVisits, performance, nPenaltyErrors, nRefErrors, and nAdjustedErrors
0030 %    indices     in the same format (N x 12 matrix) as counts
0031 %
0032 %
0033 %  NOTE
0034 %
0035 %    The option 'plots' determines whether individual and/or group analyses should
0036 %    be plotted. The option 'measures' determines which kinds of analyses should be
0037 %    plotted: possible values include 'ref','work','rewards','visits','perf','penalty',
0038 %    and 'adjusted'.
0039 
0040 % Copyright (C) 2008-2011 by Gabrielle Girardeau & Michaël Zugaro
0041 %
0042 % This program is free software; you can redistribute it and/or modify
0043 % it under the terms of the GNU General Public License as published by
0044 % the Free Software Foundation; either version 3 of the License, or
0045 % (at your option) any later version.
0046 
0047 % We use a few globals to simplify the code
0048 global showANOVAs;
0049 global measures plots;
0050 global nTrials nDays nArms nConfigurations nRats nGroups;
0051 global nVisitsDayOne;
0052 global p0 pInf;
0053 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0054 global alpha;
0055 
0056 alpha = 0.05; % Significance level for complementary tests
0057 
0058 % Column names for the 'data', 'counts' and 'indices' matrices
0059 RAT = 1;
0060 DAY = 2;
0061 TRIAL = 3;
0062 CONF = 4;
0063 GROUP = 5;
0064 
0065 % Column names for the 'data' matrix
0066 ARM = 6;
0067 
0068 % Column names for the 'counts' and 'indices' matrices
0069 REF = 6;
0070 WORK = 7;
0071 REWARDS = 8;
0072 VISITS = 9;
0073 PERF = 10;
0074 PENALTY = 11;
0075 ADJUSTED = 12;
0076 
0077 % Constants
0078 nDays = max(data(:,DAY));
0079 nTrials = 3;
0080 nArms = 8;
0081 nGroups = max(data(:,GROUP));
0082 nConfigurations = max(data(:,CONF));
0083 nRats = length(unique(data(:,RAT)));
0084 
0085 % Used later to estimate chance levels for performance
0086 nVisitsDayOne = round(Accumulate(data(:,DAY),1)/nRats/nTrials);
0087 nVisitsDayOne = nVisitsDayOne(1);
0088 
0089 % Default values
0090 allMeasures = {'ref','work','rewards','visits','perf','penalty','adjusted'};
0091 allPlots = {'groups','individuals'};
0092 defaultMeasures = {'perf'};
0093 defaultPlots = {'groups','individuals'};
0094 listMeasures = false;
0095 listPlots = false;
0096 measures = {};
0097 plots = {};
0098 showANOVAs = true;
0099 
0100 if nargin < 2,
0101     error('Incorrect number of parameters (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).');
0102 end
0103 
0104 if mod(length(varargin),2) ~= 0,
0105   error('Incorrect number of parameters (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).');
0106 end
0107 
0108 % Parse options
0109 for i = 1:2:length(varargin),
0110     if ~ischar(varargin{i}),
0111         error(['Parameter ' num2str(i+firstIndex) ' is not a property (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).']);
0112     end
0113     switch(lower(varargin{i})),
0114         case 'measures',
0115             p = lower(varargin{i+1});
0116             if ~isa(p,'cell'),
0117                 error('Value associated with property ''measures'' is not a cell array of strings (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).');
0118             end
0119             for j = 1:length(p),
0120                 if ~isastring(p{j},allMeasures{:}),
0121                     error('Incorrect value for property ''measures'' (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).');
0122                 end
0123             end
0124             measures = {measures{:},p{:}};
0125             listMeasures = true;
0126         case 'plots',
0127             f = lower(varargin{i+1});
0128             if ~isa(f,'cell'),
0129                 error('Value associated with property ''plots'' is not a cell array of strings (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).');
0130             end
0131             for j = 1:length(f),
0132             if ~isastring(f{j},allPlots{:}),
0133                     error('Incorrect value for property ''plots'' (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).');
0134                 end
0135             end
0136             plots = {plots{:},f{:}};
0137             listPlots = true;
0138         case 'anova',
0139             doIt = lower(varargin{i+1});
0140             if ~isastring(doIt,'on','off'),
0141                 error('Incorrect value for property ''anova'' (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).');
0142             end
0143             showANOVAs = strcmp(doIt,'on');
0144         otherwise,
0145             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help RadialMaze">RadialMaze</a>'' for details).']);
0146     end
0147 end
0148 
0149 if ~listMeasures, measures = defaultMeasures; end
0150 if ~listPlots, plots = defaultPlots; end
0151 
0152 % Compute stats for each rat and each configuration
0153 counts = Counts(data,baitedArms);
0154 indices = ComputeIndices(counts);
0155 
0156 %EstimatePerformanceChanceLevels;
0157 p0 = 0.481869891541284; 
0158 pInf = 0.306735773403144;
0159 
0160 h = ComputeANOVAs(indices);
0161 
0162 if any(ismember('individuals',plots)),
0163     PlotIndividualStats(counts);
0164 end
0165 if any(ismember('groups',plots)),
0166     PlotGroupStats(indices,h);
0167 end
0168 
0169 %-----------------------------------------------------------------------------------------------------------------
0170 %  Count errors (etc.)
0171 %
0172 %  INPUT
0173 %    data, baitedArms
0174 %
0175 %  OUTPUT
0176 %    counts, a N x 12 matrix, where each line lists one case (rat,day,trial,config,group,measure1,...,measureN)
0177 %    where the measures are nRefErrors, nWorkErrors, nRewardsFound, nVisits, performance, nPenaltyErrors,
0178 %    nRefErrors, and nAdjustedErrors
0179 %-----------------------------------------------------------------------------------------------------------------
0180 
0181 function counts = Counts(data,baitedArms)
0182 
0183 global nTrials nDays nArms nConfigurations nGroups;
0184 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0185 
0186 nRats = max(data(:,RAT));
0187 
0188 % Create output array and prefill it with NaNs so missing data will be readily visible
0189 % 5 variables (rat, day, trial, config, group) and 7 measures (ref, work, etc.)
0190 %  counts = nan(nRats*nDays*nTrials*nConfigurations,5+7);
0191 
0192 % Let's go...
0193 line = 1;
0194 for rat = 1:nRats,
0195 
0196     % Which group does this rat belong to?
0197     group = data(data(:,RAT)==rat,GROUP);
0198     group = group(1);
0199 
0200     for configuration = 1:nConfigurations,
0201 
0202         % Determine baited/unbaited arms for this rat and this configuration
0203         b = baitedArms(rat,configuration,:);
0204         baitedArm1 = b(1);
0205         baitedArm2 = b(2);
0206         baitedArm3 = b(3);
0207         unbaitedArms = setdiff(1:8,[baitedArm1 baitedArm2 baitedArm3]);
0208         unbaitedArm1 = unbaitedArms(1);
0209         unbaitedArm2 = unbaitedArms(2);
0210         unbaitedArm3 = unbaitedArms(3);
0211         unbaitedArm4 = unbaitedArms(4);
0212         unbaitedArm5 = unbaitedArms(5);
0213 
0214         for day = 1:nDays,
0215             for trial = 1:nTrials,
0216 
0217                 % Select the block of data corresponding the current rat, config, trial and day
0218                 block = data(data(:,RAT)==rat&data(:,CONF)==configuration&data(:,DAY)==day&data(:,TRIAL)==trial,:);
0219 
0220                 % Total number of visits
0221                 nVisits = size(block,1);
0222 
0223                 % Missing data: skip
0224                 if nVisits == 0,
0225                     continue;
0226                 end
0227 
0228                 % Set values of the variables for this line in 'counts'
0229                 counts(line,1:REF-1) = [rat,day,trial,configuration,group];
0230 
0231                 % Initially set all measures to zero (except number of visits!)
0232                 counts(line,REF:ADJUSTED) = zeros(1,7);
0233                 counts(line,VISITS) = nVisits;
0234 
0235                 % No arm visited yet
0236                 visitedBaitedArm1 = false;
0237                 visitedBaitedArm2 = false;
0238                 visitedBaitedArm3 = false;
0239                 visitedUnbaitedArm1 = false;
0240                 visitedUnbaitedArm2 = false;
0241                 visitedUnbaitedArm3 = false;
0242                 visitedUnbaitedArm4 = false;
0243                 visitedUnbaitedArm5 = false;
0244 
0245                 for visit = 1:nVisits,
0246                     % Test successively visited arms to count the number of working memory errors, reference memory errors, etc.
0247                     % In particular, check for repeated visits in baited arms (returning in an already visited baited arm is considered a working memory error)
0248                     currentArm = block(visit,ARM);
0249                     switch currentArm,
0250                         case baitedArm1,
0251                             if visitedBaitedArm1,
0252                                 % Returning in an already visited baited arm is considered as a working memory error
0253                                 counts(line,WORK) = counts(line,WORK)+1;
0254                             else
0255                                 visitedBaitedArm1 = true;
0256                                 counts(line,REWARDS) = counts(line,REWARDS)+1;
0257                             end
0258                         case baitedArm2,
0259                             if visitedBaitedArm2,
0260                                 counts(line,WORK) = counts(line,WORK)+1;
0261                             else
0262                                 visitedBaitedArm2 = true;
0263                                 counts(line,REWARDS) = counts(line,REWARDS)+1;
0264                             end
0265                         case baitedArm3,
0266                             if visitedBaitedArm3,
0267                                 counts(line,WORK) = counts(line,WORK)+1;
0268                             else
0269                                 visitedBaitedArm3 = true;
0270                                 counts(line,REWARDS) = counts(line,REWARDS)+1;
0271                             end
0272                         case unbaitedArm1,
0273                             if visitedUnbaitedArm1,
0274                                 counts(line,WORK) = counts(line,WORK)+1;
0275                             else
0276                                 visitedUnbaitedArm1 = true;
0277                                 counts(line,REF) = counts(line,REF)+1;
0278                             end
0279                         case unbaitedArm2,
0280                             if visitedUnbaitedArm2,
0281                                 counts(line,WORK) = counts(line,WORK)+1;
0282                             else
0283                                 visitedUnbaitedArm2 = true;
0284                                 counts(line,REF) = counts(line,REF)+1;
0285                             end
0286                         case unbaitedArm3,
0287                             if visitedUnbaitedArm3,
0288                                 counts(line,WORK) = counts(line,WORK)+1;
0289                             else
0290                                 visitedUnbaitedArm3 = true;
0291                                 counts(line,REF) = counts(line,REF)+1;
0292                             end
0293                         case unbaitedArm4,
0294                             if visitedUnbaitedArm4,
0295                                 counts(line,WORK) = counts(line,WORK)+1;
0296                             else
0297                                 visitedUnbaitedArm4 = true;
0298                                 counts(line,REF) = counts(line,REF)+1;
0299                             end
0300                         case unbaitedArm5,
0301                             if visitedUnbaitedArm5,
0302                                 counts(line,WORK) = counts(line,WORK)+1;
0303                             else
0304                                 visitedUnbaitedArm5 = true;
0305                                 counts(line,REF) = counts(line,REF)+1;
0306                             end
0307                     end
0308                 end
0309                 counts(line,PERF) = counts(line,REWARDS)/counts(line,VISITS);
0310                 % Adjusted errors = visits in non-rewarded arms + unvisited baited arms
0311                 counts(line,ADJUSTED) = counts(line,REF)+(3-counts(line,REWARDS));
0312                 % Penalty errors = 5 if the rat did not find all the rewards
0313                 if counts(line,REWARDS)~= 3,
0314                     counts(line,PENALTY) = 5;
0315                 else
0316                     counts(line,PENALTY) = counts(line,REF);
0317                 end
0318                 line = line + 1;
0319             end
0320         end
0321     end
0322 end
0323 
0324 
0325 %-----------------------------------------------------------------------------------------------------------------
0326 %  Transform counts/proportions for ANOVA (square root transform, arcsine square root transform)
0327 %
0328 %  INPUT
0329 %    counts
0330 %
0331 %  OUTPUT
0332 %    indices, in the same format (N x 12 matrix) as counts
0333 %-----------------------------------------------------------------------------------------------------------------
0334 
0335 function indices = ComputeIndices(counts)
0336 
0337 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0338 
0339 indices = counts;
0340 indices(:,[REF WORK REWARDS VISITS PENALTY ADJUSTED]) = ComputeCountIndex(counts(:,[REF WORK REWARDS VISITS PENALTY ADJUSTED]));
0341 indices(:,PERF) = ComputeProportionIndex(counts(:,PERF));
0342 
0343 function index = ComputeCountIndex(count)
0344 
0345 index = sqrt(count+1);
0346 
0347 function index = ComputeProportionIndex(proportion)
0348 
0349 index = asin(sqrt(proportion))/(pi/2);
0350 
0351 %-----------------------------------------------------------------------------------------------------------------
0352 %  Individual plots: plot error counts for each rat in each configuration
0353 %
0354 %  INPUT
0355 %    counts
0356 %-----------------------------------------------------------------------------------------------------------------
0357 
0358 function PlotIndividualStats(counts)
0359 
0360 global nTrials nDays nConfigurations measures plots testRats controlRats;
0361 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0362 
0363 if any(ismember('ref',measures)),
0364     DoPlotIndividualStats(counts,REF,'# ref mem errors',[0 5]);
0365 end
0366 if any(ismember('work',measures)),
0367     DoPlotIndividualStats(counts,WORK,'# work mem errors',[0 5]);
0368 end
0369 if any(ismember('visits',measures)),
0370     DoPlotIndividualStats(counts,VISITS,'# visits',[0 10]);
0371 end
0372 if any(ismember('rewards',measures)),
0373     DoPlotIndividualStats(counts,REWARDS,'# rewards',[0 3]);
0374 end
0375 if any(ismember('perf',measures)),
0376     DoPlotIndividualStats(counts,PERF,'performance',[0 1]);
0377 end
0378 if any(ismember('penalty',measures)),
0379     DoPlotIndividualStats(counts,PENALTY,'# penalty errors',[0 5]);
0380 end
0381 if any(ismember('adjusted',measures)),
0382     DoPlotIndividualStats(counts,ADJUSTED,'# adjusted errors',[0 5]);
0383 end
0384 
0385 
0386 function DoPlotIndividualStats(counts,measure,measureStr,yLim)
0387 
0388 global nRats nConfigurations nDays;
0389 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0390 
0391 p0 = 0.481869891541284; 
0392 pInf = 0.306735773403144;
0393 
0394 colors = 'rbkgcmy'; % Hopefully there will not be more than 7 groups!
0395 for configuration = 1:nConfigurations,
0396     figure;
0397     for rat = 1:nRats,
0398         SquareSubplot(nRats,rat);
0399         d = counts(counts(:,RAT)==rat&counts(:,CONF)==configuration,[DAY measure GROUP]);
0400         if isempty(d), continue; end
0401         y = Accumulate(d(:,1),d(:,2));
0402         n = Accumulate(d(:,1),1);
0403         plot(y./n,colors(d(1,3)));
0404         xlabel(['Rat #' int2str(rat) ' - Configuration #' int2str(configuration)]);
0405         ylabel(measureStr);
0406         xlim([0 nDays+1]);set(gca,'xtick',1:nDays);ylim(yLim);
0407         PlotIntervals([pInf,p0],'rectangles','h');
0408     end
0409 end
0410 
0411 %-----------------------------------------------------------------------------------------------------------------
0412 %  Group plots: plot error indices for all rats in each configuration
0413 %
0414 %  INPUT
0415 %    indices
0416 %    h   post-hoc ttest results from ANOVAs
0417 %-----------------------------------------------------------------------------------------------------------------
0418 
0419 function PlotGroupStats(indices,h)
0420 
0421 global measures;
0422 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0423 
0424 % Reference memory errors
0425 if any(ismember('ref',measures)),
0426     DoPlotGroupStats(indices,REF,'index ref mem errors',[0 5],h);
0427 end
0428 if any(ismember('work',measures)),
0429     DoPlotGroupStats(indices,WORK,'index work mem errors',[0.5 2],h);
0430 end
0431 if any(ismember('rewards',measures)),
0432     DoPlotGroupStats(indices,REWARDS,'index rewards errors',[0.3 1],h);
0433 end
0434 if any(ismember('visits',measures)),
0435     DoPlotGroupStats(indices,VISITS,'index visits errors',[0.3 1],h);
0436 end
0437 if any(ismember('perf',measures)),
0438     DoPlotGroupStats(indices,PERF,'index performance errors',[0.3 1],h);
0439 end
0440 if any(ismember('penalty',measures)),
0441     DoPlotGroupStats(indices,PENALTY,'index penalty errors',[0.3 1],h);
0442 end
0443 if any(ismember('adjusted',measures)),
0444     DoPlotGroupStats(indices,ADJUSTED,'index adjusted errors',[0.3 1],h);
0445 end
0446 
0447 
0448 function DoPlotGroupStats(indices,measure,measureStr,yLim,h)
0449 
0450 global nTrials nDays nConfigurations nRats nGroups;
0451 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0452 global p0 pInf;
0453 
0454 colors = 'rbkgcmy'; % Hopefully there will not be more than 7 groups!
0455 for configuration = 1:nConfigurations,
0456     figure;hold on;
0457     for group = 1:nGroups,
0458         d = indices(indices(:,GROUP)==group&indices(:,CONF)==configuration,[DAY measure]);
0459         if isempty(d), continue; end
0460         [m,s] = Accumulate(d(:,1),d(:,2));
0461         n = Accumulate(d(:,1),1);
0462         sem = s./sqrt(n);
0463         p = errorbar(m,sem,colors(group));
0464         set(p,'marker','.','markersize',24);
0465         xlabel(['Configuration #' int2str(configuration)]);
0466         ylabel(['index ' measureStr]);
0467         xlim([0 nDays+1]);set(gca,'xtick',1:nDays);ylim(yLim);
0468 
0469         if ~isempty(h),
0470             x = find(h(:,1+group*2,measure));
0471             plot(x,(m(x)-sem(x))-.025,'k+','markersize',10);
0472             if group == 1 & configuration == 1,
0473                 x = find(h(:,1,measure));
0474                 plot(x,(m(x)+sem(x))+.025,'k*','markersize',10);
0475             end
0476         end
0477     end
0478     if measure == PERF,
0479         PlotIntervals([pInf,p0],'rectangles','h');
0480     end
0481 end
0482 
0483 
0484 %-----------------------------------------------------------------------------------------------------------------
0485 %  Compute ANOVAs
0486 %
0487 %  INPUT
0488 %    indices
0489 %-----------------------------------------------------------------------------------------------------------------
0490 
0491 function h = ComputeANOVAs(indices)
0492 
0493 global measures;
0494 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0495 
0496 if any(ismember('ref',measures)),
0497     h(:,:,REF) = DoComputeANOVAs(indices,REF,'ref mem errors');
0498 end
0499 if any(ismember('work',measures)),
0500     h(:,:,WORK) = DoComputeANOVAs(indices,WORK,'work mem errors');
0501 end
0502 if any(ismember('visits',measures)),
0503     h(:,:,VISITS) = DoComputeANOVAs(indices,VISITS,'visits');
0504 end
0505 if any(ismember('rewards',measures)),
0506     h(:,:,REWARDS) = DoComputeANOVAs(indices,REWARDS,'rewards');
0507 end
0508 if any(ismember('perf',measures)),
0509     h(:,:,PERF) = DoComputeANOVAs(indices,PERF,'performance');
0510 end
0511 if any(ismember('penalty',measures)),
0512     h(:,:,PENALTY) = DoComputeANOVAs(indices,PENALTY,'penalty errors');
0513 end
0514 if any(ismember('adjusted',measures)),
0515     h(:,:,ADJUSTED) = DoComputeANOVAs(indices,ADJUSTED,'adjusted errors');
0516 end
0517 
0518 
0519 function h = DoComputeANOVAs(indices,measure,measureStr)
0520 
0521 global showANOVAs;
0522 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0523 
0524 [p,t,stats,terms,h] = ComputeTwoWayANOVA(indices,measure);
0525 if showANOVAs,
0526     f = statdisptable(t,['ANOVA for ' measureStr],'Analysis of Variance');
0527     set(f,'numbertitle','off');
0528 end
0529 
0530 %-----------------------------------------------------------------------------------------------------------------
0531 %  Compute day x group ANOVA on indices for a given measure (ref errors, work errors, etc.)
0532 %  !!! Only for configuration 1 !!!
0533 %
0534 %  INPUT
0535 %    indices
0536 %
0537 %  OUTPUT
0538 %    p,t,stats,terms   ANOVA results
0539 %    h                 = [H0,p0,H1,p1,...Hn,pn] where [H0,p0] is the result of a t-test comparing groups 1 and 2
0540 %                      (one row per day), and [Hi,pi] is the result of a t-test comparing group i to the upper
0541 %                      chance level
0542 %-----------------------------------------------------------------------------------------------------------------
0543 
0544 function [p,t,stats,terms,h] = ComputeTwoWayANOVA(indices,measure)
0545 
0546 global RAT DAY TRIAL CONF GROUP ARM REF WORK REWARDS VISITS PERF PENALTY ADJUSTED;
0547 global p0 pInf;
0548 global nDays nGroups alpha;
0549 
0550 i = indices(:,CONF) == 1;
0551 [p,t,stats,terms] = anovan(indices(i,measure),{indices(i,GROUP) indices(i,DAY)},'model','full','varnames',{'group','day'},'display','off');
0552 
0553 h = zeros(nDays,2+2*nGroups);
0554 switch measure,
0555     case { REF, WORK, VISITS, PENALTY, ADJUSTED },
0556         tail = 'right';
0557     case { PERF, REWARDS },
0558         tail = 'left';
0559 end
0560 
0561 % Compare groups 1 and 2, and each group to chance level
0562 if all(p<0.05),
0563     for day = 1:nDays,
0564         d1 = indices(i&indices(:,DAY)==day&indices(:,GROUP)==1,measure);
0565         d2 = indices(i&indices(:,DAY)==day&indices(:,GROUP)==2,measure);
0566         if ~isempty(d1) && ~isempty(d2),
0567 %              keyboard
0568 %              [h(day,1),h(day,2)] = ttest2(d1,d2,alpha,tail,'equal');
0569             [h(day,1),h(day,2)] = ttest2(d1,d2,alpha,tail);
0570         end
0571         for group = 1:nGroups,
0572             d0 = indices(i&indices(:,DAY)==day&indices(:,GROUP)==group,measure);
0573             if ~isempty(d0),
0574                 [h(day,1+2*group),h(day,2+2*group)] = ttest(d0,p0,alpha,'right');
0575             end
0576         end
0577     end
0578 end
0579 
0580 %-----------------------------------------------------------------------------------------------------------------
0581 %  DEPRECATED CODE
0582 %-----------------------------------------------------------------------------------------------------------------
0583 %  Compute repeated measures ANOVA on indices for a given measure (ref errors, work errors, etc.)
0584 %  This can be used to compare e.g. test rats and control rats on the same configuration, or test rats on both
0585 %  configurations
0586 %
0587 %  INPUT
0588 %    data1           'indices' for group 1
0589 %    data2           'indices' for group 2
0590 %-----------------------------------------------------------------------------------------------------------------
0591 
0592 function ComputeRepeatedMeasuresANOVA(data1,data2)
0593 
0594 global nTrials nDays;
0595 
0596 group1 = ones(size(data1));
0597 n = size(data1,2);
0598 days1 = repmat((1:nDays)',1,n);
0599 trials1 = repmat((1:nTrials),nDays,n/3);
0600 rats1 = repmat(ceil((1:n)/3),nDays,1);
0601 
0602 group2 = 2*ones(size(data2));
0603 n = size(data2,2);
0604 days2 = repmat((1:nDays)',1,n);
0605 trials2 = repmat((1:nTrials),nDays,n/3);
0606 rats2 = repmat(ceil((1:n)/3),nDays,1);
0607 
0608 data = [data1(:);data2(:)];
0609 groups = [group1(:);group2(:)];
0610 days = [days1(:);days2(:)];
0611 trials = [trials1(:);trials2(:)];
0612 rats = [rats1(:);rats2(:)];
0613 
0614 nans = isnan(data(:,1));
0615 data(nans,:) = [];
0616 groups(nans,:) = [];
0617 days(nans,:) = [];
0618 trials(nans,:) = [];
0619 rats(nans,:) = [];
0620 
0621 d = [data groups rats days];
0622 
0623 RMAOV1MS(d);
0624 
0625 %-----------------------------------------------------------------------------------------------------------------
0626 %  Estimate chance levels for performance
0627 %-----------------------------------------------------------------------------------------------------------------
0628 
0629 function EstimatePerformanceChanceLevels
0630 
0631 global p0 pInf;
0632 
0633 p0 = mean(bootstrp(10000,@SimulateOneTrial,0));
0634 pInf = mean(bootstrp(10000,@SimulateOneTrial,Inf));
0635 
0636 % nWME = # working memory errors
0637 % If nWME is Inf, then the rat can choose any arm on each trial, independent of previous choices
0638 % Otherwise, the performance is computed in two steps: 1) the rat never makes working memory errors
0639 % and 2) a fixed number of working memory errors are added to the performance computed in 1).
0640 
0641 function index = SimulateOneTrial(nWME)
0642 
0643 global nVisitsDayOne;
0644 
0645 if isinf(nWME),
0646     % Draw n integers from 1 to 8 (= arm number)
0647     n = nVisitsDayOne;
0648     r = floor(Clip(rand(n,1),0,1-eps)*8)+1;
0649     % How many 'arms' came out from the list [1 2 3]?
0650     found = any(r==1) + any(r==2) + any(r==3);
0651 else
0652     % Find arm #1
0653     r = randperm(8);
0654     n1 = find(r==1|r==2|r==3);n1 = n1(1);
0655     % Find arm #2
0656     r = randperm(8-n1);
0657     n2 = find(r==1|r==2);n2 = n2(1);
0658     % Find arm #3
0659     r = randperm(8-n1-n2);
0660     n3 = find(r==1);
0661     found = 3;
0662     n = n1+n2+n3+nWME;
0663 end
0664 % Corresponding proportion index
0665 index = ComputeProportionIndex(found/n);