octave_packages/statistics-1.1.3/repanova.m

   1 ## Copyright (C) 2011 Kyle Winfree <kyle.winfree@gmail.com>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{pval}, @var{table}, @var{st}] =} repanova (@var{X}, @var{cond})
  18 ## @deftypefnx {Function File} {[@var{pval}, @var{table}, @var{st}] =} repanova (@var{X}, @var{cond}, ['string' | 'cell'])
  19 ## Perform a repeated measures analysis of variance (Repeated ANOVA).
  20 ## X is formated such that each row is a subject and each column is a condition.
  21 ##
  22 ## condition is typically a point in time, say t=1 then t=2, etc
  23 ## condition can also be thought of as groups.
  24 ##
  25 ## The optional flag can be either 'cell' or 'string' and reflects
  26 ## the format of the table returned.  Cell is the default.
  27 ##
  28 ## NaNs are ignored using nanmean and nanstd.
  29 ##
  30 ## This fuction does not currently support multiple columns of the same
  31 ## condition!
  32 ## @end deftypefn
  33
  34 function [p, table, st] = repanova(varargin)
  35
  36 switch nargin
  37      case 0
  38           error('Too few inputs.');
  39      case 1
  40                 X = varargin{1};
  41                 for c = 1:size(X, 2)
  42                         condition{c} = ['time', num2str(c)];
  43                 end
  44                 option = 'cell';
  45      case 2
  46           X = varargin{1};
  47           condition = varargin{2};
  48           option = 'cell';
  49      case 3
  50           X = varargin{1};
  51           condition = varargin{2};
  52           option = varargin{3};
  53      otherwise
  54           error('Too many inputs.');
  55 end
  56      % Find the means of the subjects and measures, ignoring any NaNs
  57      u_subjects = nanmean(X,2);
  58      u_measures = nanmean(X,1);
  59      u_grand = nansum(nansum(X)) / (size(X,1) * size(X,2));
  60      % Differences between rows will be reflected in SS subjects, differences
  61      % between columns will be reflected in SS_within subjects.
  62      N = size(X,1); % number of subjects
  63      J = size(X,2); % number of samples per subject
  64      SS_measures = N * nansum((u_measures - u_grand).^2);
  65      SS_subjects = J * nansum((u_subjects - u_grand).^2);
  66      SS_total = nansum(nansum((X - u_grand).^2));
  67      SS_error = SS_total - SS_measures - SS_subjects;
  68      df_measures = J - 1;
  69      df_subjects = N - 1;
  70      df_grand = (N*J) - 1;
  71      df_error = df_grand - df_measures - df_subjects;
  72      MS_measures = SS_measures / df_measures;
  73      MS_subjects = SS_subjects / df_subjects;
  74      MS_error = SS_error / df_error; % variation expected as a result of sampling error alone
  75      F = MS_measures / MS_error;
  76      p = 1 - fcdf(F, df_measures, df_error); % Probability of F given equal means.
  77
  78      if strcmp(option, 'string')
  79           table  = [sprintf('\nSource\tSS\tdf\tMS\tF\tProb > F'), ...
  80                     sprintf('\nSubject\t%g\t%i\t%g', SS_subjects, df_subjects, MS_subjects), ...
  81                     sprintf('\nMeasure\t%g\t%i\t%g\t%g\t%g', SS_measures, df_measures, MS_measures, F, p), ...
  82                     sprintf('\nError\t%g\t%i\t%g', SS_error, df_error, MS_error), ...
  83                     sprintf('\n')];
  84      else
  85           table  = {'Source', 'Partial SS', 'df', 'MS', 'F', 'Prob > F'; ...
  86                     'Subject', SS_subjects, df_subjects, MS_subjects, '', ''; ...
  87                     'Measure', SS_measures, df_measures, MS_measures, F, p};
  88      end
  89
  90      st.gnames = condition'; % this is the same struct format used in anova1
  91      st.n = repmat(N, 1, J);
  92      st.source = 'anova1'; % it cannot be assumed that 'repanova' is a supported source for multcompare
  93      st.means = u_measures;
  94      st.df = df_error;
  95      st.s = sqrt(MS_error);
  96 end
  97
  98 % This function was created with guidance from the following websites:
  99 % http://courses.washington.edu/stat217/rmANOVA.html
 100 % http://grants.hhp.coe.uh.edu/doconnor/PEP6305/Topic%20010%20Repeated%20Measures.htm