Home > src > sosCorrelTest.m

# sosCorrelTest

## PURPOSE

- test for equivalence of two correlations

## SYNOPSIS

This is a script file.

## DESCRIPTION

     - test for equivalence of two correlations

     copyright 2009-2012 Blair Armstrong, Christine Watson, David Plaut

        This file is part of SOS

        SOS is free software: you can redistribute it and/or modify
        it for academic and non-commercial purposes
        under the terms of the GNU General Public License as published by
        the Free Software Foundation, either version 3 of the License, or
        (at your option) any later version.  For commercial or for-profit
        uses, please contact the authors (sos@cnbc.cmu.edu).

        SOS is distributed in the hope that it will be useful,
        but WITHOUT ANY WARRANTY; without even the implied warranty of
        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
        GNU General Public License for more details.

## CROSS-REFERENCE INFORMATION

This function calls:

- sosCorrelTest - test for equivalence of two correlations
- verbosePrint - master print function that regulates cmd line output based on verbosity flags

This function is called by:

- genericStatTest - parent class for statistical test objects
- sosCorrelTest - test for equivalence of two correlations

## SUBFUNCTIONS

- function obj = sosCorrelTest(varargin)
- function [userHypothesis, prob, label] = runTest(obj,varargin)
- function [userHypothesis, prob, label] =
- function constructMatchCorrelztest(obj,sosObj,varargin)
- function userHypothesis = returnNaN(~,~)
- function flag = validTestType(str)

## SOURCE CODE

    0001 % - test for equivalence of two correlations
    0002 %
    0003 % copyright 2009-2012 Blair Armstrong, Christine Watson, David Plaut
    0004 %
    0005 %    This file is part of SOS
    0006 %
    0007 %    SOS is free software: you can redistribute it and/or modify
    0008 %    it for academic and non-commercial purposes
    0009 %    under the terms of the GNU General Public License as published by
    0010 %    the Free Software Foundation, either version 3 of the License, or
    0011 %    (at your option) any later version.  For commercial or for-profit
    0012 %    uses, please contact the authors (sos@cnbc.cmu.edu).
    0013 %
    0014 %    SOS is distributed in the hope that it will be useful,
    0015 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
    0016 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0017 %    GNU General Public License for more details.
    0018
    0019 %    You should have received a copy of the GNU General Public License
    0020 %    along with SOS (see COPYING.txt).
    0021 %    If not, see <http://www.gnu.org/licenses/>.
    0022
    0023
    0024 classdef sosCorrelTest < handle
    0025     % runs user-specified z-tests on equivalence of a sample and a specified
    0026     % target (population) correlation coefficient,
    0027     % and tests user hypotheses on their outcomes
    0028     % NOTE: This uses a z-score approximation when testing the equivalence
    0029     %   of the sample and population correlation coefficients.  Using
    0030     %   sample sizes < 10 is not recommended, and a sample size of at least
    0031     %   4 is required to have sufficient df for the test.
    0032     %   Additionally, you should avoid trying to match to a correlation of
    0033     %   1.0, beceause that is undefined in the equation.
    0034     %   Further, note that for a same p-value, closer matches are required
    0035     %   nearer to a correlation of 1.0 than to a correlation of 0, by
    0036     %   virtue of the nonlinearity which is part of the stats test.
    0037     %
    0038     %PROPERTIES
    0039     %     sosObj % sos object the test is associated with
    0040     %     name % string label to associate with the test
    0041     %     s1 % sample1
    0042     %     s2 % sample2
    0043     %     s1ColName % column name of data in sample1
    0044     %     s2ColName % column name of data in sample2
    0045     %     s1Col   % column index in sample1
    0046     %     s2Col   % column index in sample2
    0047     %     type    % type of ztest (currently only 'matchCorrel')
    0048     %     runSpecificTest % handle to specific test method
    0049     %     desiredpvalCondition % indicates the desired outcome of the test, either NaN, <= or >=
    0050     %     desiredpvalConditionHandle  % handle to the function that will evaluate whether the desiredpValCondition was met
    0051     %     desiredpval % desired p-val to test against
    0052     %     targValue % target correlation value for the test
    0053     %     tail % tail of test.  left -> t-vals < 0, right = t-vals > 0; both = pos. and neg. z-score will reject the null
    0054     %     lastp % last p-value that was calcualted
    0055     %
    0056     %METHODS:
    0057     %   sosCorrelTest(sosObj,varargin)  %constructs an sosCorrelTest object of the specified type.
    0058     %   constructMatchCorrelztest(obj,sosObj,varargin)  %initialize object
    0059     %   [userHypothesis, prob, label] = runTest(obj,varargin) %runs the test
    0060     %   [userHypothesis, prob, label] =  runMatchCorrelztest(obj,varargin) % runs a matchCorrel test
    0061     %
    0062     %METHODS (Static)
    0063     %   userHypothesis = returnNaN(~,~)  % returns NaN
    0064     %   flag = validTestType(str) % returns 1 if the name of type of test is 'matchCorrel'; error otherwise.
    0065
    0066
    0067     %% PROPERTIES
    0068     properties
    0069         sosObj % sos object the test is associated with
    0070         name % string label to associate with the test
    0071         s1 % sample1
    0072         s2 % sample2
    0073         s1ColName % column name of data in sample1
    0074         s2ColName % column name of data in sample2
    0075         s1Col   % column index in sample1
    0076         s2Col   % column index in sample2
    0077         type    % type of ttest (paired,independent,
    0078         runSpecificTest % handle to specific test method
    0079         desiredpvalCondition % indicates the desired outcome of the test, either NaN, <= or >=
    0080         desiredpvalConditionHandle  % handle to the function that will evaluate whether the desiredpValCondition was met
    0081         desiredpval % desired p-val to test against
    0082         targVal % target value for single-sample ttest
    0083         tail % tail of test.  left -> t-vals < 0, right = t-vals > 0; both = pos. and neg. t will reject the null
    0084         lastp % last p-value that was calcualted
    0085         label % label denoting what is being tested
    0086     end % properties
    0087
    0088     methods
    0089
    0090         %% sosCorrelTest CONSTRUCTOR
    0091         function obj = sosCorrelTest(varargin)
    0092             %constructs an sosCorrelTest object of the specified type.
    0093             %
    0094             %PARAMETERS:
    0095             % sosObj - sosObject test will be associated with
    0096             % 'type'/string, type of test - ('matchCorrel' only currently)
    0097             % 'sample1'/sample - a sample associated with sosObj.
    0098             % ... Other parameters as required by the constructor for the
    0099             % specific type of t-test requested.  See those constructors
    0100             % for additional info.
    0101             %
    0102             % Returns an sosCorrelTest object.
    0103
    0104
    0105             % perform basic checks for universally required components of
    0106             % the z-test
    0107             p = inputParser;
    0108             p.addParamValue('sosObj', 'null',...
    0109                 @(sosObj)strcmp(class(sosObj),'sos'));
    0110             p.addParamValue('type','null', ...
    0111                 @(type)sosCorrelTest.validTestType(type));
    0112             p.addParamValue('sample1','null', ...
    0113                             @(sample1)strcmp(class(sample1),'sample'));
    0114             p.addParamValue('name','noname',@(name)ischar(name));
    0115
    0116             %keep all of the parameters for specific t-tests
    0117             p.KeepUnmatched = true;
    0118             p.parse(varargin{:});
    0119
    0120
    0121
    0122             %basic checks are passed, construct the specific type of t-test
    0123             %that has been requested
    0124
    0125            sosObj = p.Results.sosObj; %#ok<PROP>
    0126
    0127             if(strcmp(p.Results.type,'matchCorrel'))
    0128                 obj.constructMatchCorrelztest(sosObj,varargin{:}) %#ok<PROP>
    0129             else
    0130                error(['Specified <type>: ',p.Results.type, ...
    0131                         ' is not supported']);
    0132             end
    0133
    0134             obj.lastp = NaN;
    0135             obj.sosObj = p.Results.sosObj;
    0136
    0137             % add the test's name.
    0138             if any(strcmp(p.UsingDefaults,'name'))
    0139
    0140                  numTests = length(obj.sosObj.sosstattests);
    0141
    0142              obj.name = ['ztest_',num2str(numTests+1)];
    0143              else
    0144                  obj.name = p.Results.name;
    0145              end
    0146
    0147         end % Constructor
    0148
    0149
    0150         %% [userHypothesis, prob, label] = runTest(obj,varargin) METHOD
    0151         function [userHypothesis, prob, label] = runTest(obj,varargin)
    0152             %runs the ztest
    0153             %
    0154             % PARAMETERS:
    0155             % 'reportStyle'/'short'|'full'/'none' - style of report to be
    0156             %           printed.  Either none, short or full
    0157             %
    0158             % RETURNS:
    0159             %   userHypothesis - whether the user's hypothesis has been met or not. NaN if no user hypothesis for the test.
    0160             %   prob - p-value from the ttest
    0161             %   label - string label denoting what was tested.
    0162
    0163             [userHypothesis, prob, label] = obj.runSpecificTest(varargin);
    0164             obj.lastp = prob;
    0165         end
    0166
    0167
    0168         %% [userHypothesis, prob, label] =  runMatchCorrelztest(varargin) METHOD
    0169         function [userHypothesis, prob, label] = ...
    0170                 runMatchCorrelztest(obj,varargin)
    0171             %runs a matchCorrel test
    0172             %
    0173             % PARAMETERS:
    0174             % 'reportStyle'/'short'|'full' - style of report to be printed.  Either short or long
    0175             %
    0176             % RETURNS:
    0177             %   userHypothesis - whether the user's hypothesis has been met or not. NaN if no user hypothesis for the test.
    0178             %   prob - p-value from the ttest
    0179             %   label - string label denoting what was tested.
    0180
    0181
    0182             varargin = varargin{1};
    0183
    0184             p = inputParser;
    0185
    0186             p.addParamValue('reportStyle','short', ...
    0187                     @(reportStyle)any([strcmp(reportStyle,'short') ...
    0188                                     strcmp(reportStyle,'full') ...
    0189                                     strcmp(reportStyle,'none')]));
    0190             p.parse(varargin{:});
    0191
    0192
    0193             reportStyle = p.Results.reportStyle;
    0194
    0195
    0196
    0197             % Run the stats test
    0198             sumx1sq = sum((obj.s1.zdata{obj.s1Col}).^2);
    0199             sumx1 = sum((obj.s1.zdata{obj.s1Col}));
    0200             n1 = length(obj.s1.zdata{obj.s1Col});
    0201
    0202             sumx2sq = sum((obj.s2.zdata{obj.s2Col}).^2);
    0203             sumx2 = sum((obj.s2.zdata{obj.s2Col}));
    0204             n2 = length(obj.s2.zdata{obj.s2Col});
    0205
    0206             sumx1x2 = (obj.s1.zdata{obj.s1Col})' * (obj.s2.zdata{obj.s2Col});
    0207
    0208             ssx1 = sumx1sq - ((sumx1)^2)/n1;
    0209             ssx2 = sumx2sq - ((sumx2)^2)/n2;
    0210
    0211             spx1x2 = sumx1x2 - sumx1*sumx2/n1;
    0212
    0213             % defining overrides...
    0214             if(ssx1 == 0 && ssx2 == 0) % no variance in either condition, define correlation as being perfect in this case
    0215                 cor = 1;
    0216             elseif(ssx1 == 0 || ssx2 == 0)
    0217                 cor = 0;
    0218             else % compute correlation normally
    0219                 cor = spx1x2/(sqrt(ssx1)*sqrt(ssx2));
    0220             end
    0221
    0222             % convert correlation to z-score
    0223
    0224             r= cor;
    0225
    0226             if n1 > 4
    0227                 % we can compute the z-score
    0228
    0229                 z = (0.5*log((1+r)/(1-r)) - 0.5*log((1+obj.targVal)/(1-obj.targVal))) ...
    0230                         /(1/sqrt(n1-3));
    0231                 [h,prob,ci,stats] = ztest(z,0,1);     %#ok<ASGLU>
    0232             else
    0233                 %we can't convert sinze formula has sqrt of N-3 in denom
    0234
    0235                 %special override cases...
    0236             end
    0237
    0238
    0239 % From ttest code.  Can remove.
    0240 %             % if prob was NaN, it's because the denominator in the t-test
    0241 %             % was 0.  So just check the means manually to determine
    0242 %             % prob
    0243 %             if isnan(prob)
    0244 %              m1 = mean(obj.s1.data{obj.s1Col});
    0245 %              m2 = mean(obj.s2.data{obj.s2Col});
    0246 %
    0247 %              if m1==m2
    0248 %                  prob = 1;
    0249 %              else
    0250 %                  prob = 0;
    0251 %              end
    0252 %
    0253 %             end
    0254
    0255             userHypothesis = obj.desiredpvalConditionHandle(...
    0256                                   prob,obj.desiredpval);
    0257
    0258              label = [obj.s1.name, '{',obj.s1ColName, '}-', ...
    0259                     obj.s2.name, '{',obj.s2ColName, '}'];
    0260
    0261              if (isnan(userHypothesis))
    0262                  printHyp = 'N/A';
    0263              elseif userHypothesis == 1
    0264                  printHyp = 'PASS';
    0265              elseif userHypothesis == 0
    0266                  printHyp = 'FAIL';
    0267              end
    0268
    0269
    0270              if (strcmp(reportStyle,'short'))
    0271                 verbosePrint([' UserHyp: ', printHyp, '; ', label, ': ', ...
    0272                      'z[',obj.type,'](',num2str(n1),') = ', ...
    0273                      num2str(stats), ', p = ', num2str(prob), ...
    0274                      ' p-des: ',num2str(obj.desiredpval)], ...
    0275                      'sosCorrelTest_runMatchCorrelztest');
    0276              elseif (strcmp(reportStyle,'full'))
    0277                  verbosePrint([' UserHyp: ', printHyp , ...
    0278                      '; ', label, ': ', ...
    0279                      'z[',obj.type,'](',num2str(n1),') = ', ...
    0280                      num2str(stats), ', p = ', num2str(prob), ...
    0281                      ' p-des: ',num2str(obj.desiredpval), ...
    0282                      ' cor = ', num2str(r),' targCor = ', num2str(obj.targVal)], ...
    0283                      'sosCorrelTest_runMatchCorrelztest');
    0284              end
    0285         end
    0286
    0287
    0288
    0289
    0290         %% constructMatchCorrelztest(sosObj,varargin) METHOD
    0291         function constructMatchCorrelztest(obj,sosObj,varargin)
    0292             %initialize an matchCorrel z-test object.
    0293             % Assumes equal variance.
    0294             %
    0295             % PARAMETERS:
    0296             % Required:
    0297             %   sosObj - sos Object test is to be associated with
    0298             %   sample1 - a sample object
    0299             %   sample2 - a sample object
    0300             %   s1ColName - name of data column in sample1
    0301             %   s2ColName - name of data column in sample2
    0302             %
    0303             % Optional:
    0304             %   desiredpvalCondition/string - desired condition for the ttest
    0305             %       pval. Either it should exceed (=>) some value, be less
    0306             %       than '<=' some condition, or be 'N/A' if there is no
    0307             %       desired condition.  Default is N/A.  Note that the
    0308             %       ordering of '=' and '<' is important, so though '<=' is
    0309             %       valid, '=<' is not.
    0310             %   desiredpval - desired p-value to evaluate the condition
    0311             %       against.  Defaults to 0.05
    0312             %   tail - tail of test - left/right/both
    0313             %   targVal - the target correlation value to match (must be in
    0314             %   range (-1,1) (note round, not square brackets, as -1 and 1
    0315             %   are NOTE allowed)
    0316             %
    0317             % RETURNS:
    0318             %   Configured sosCorrelTest object.
    0319
    0320
    0321             p = inputParser;
    0322
    0323             p.addRequired('sosObj', ...
    0324                 @(sosObj)strcmp(class(sosObj),'sos'));
    0325             p.addParamValue('type','null', ...
    0326                 @(type)sosCorrelTest.validTestType(type));
    0327             p.addParamValue('sample1','null', ...
    0328                         @(sample1)strcmp(class(sample1),'sample'));
    0329             p.addParamValue('sample2','null', ...
    0330                         @(sample2)strcmp(class(sample2),'sample'));
    0331             %NaN will fail by default
    0332             p.addParamValue('s1ColName',NaN, ...
    0333                 @(s1ColName)ischar(s1ColName));
    0334             p.addParamValue('s2ColName',NaN, ...
    0335                 @(s2ColName)ischar(s2ColName));
    0336             p.addParamValue('desiredpvalCondition','N/A', ...
    0337                 @(desiredpvalCondition)any([strcmp(desiredpvalCondition,'<='), ...
    0338                             strcmp(desiredpvalCondition,'=>'), ...
    0339                             strcmp(desiredpvalCondition,'N/A')]));
    0340             p.addParamValue('desiredpval', 0.05, ...
    0341                 @(desiredpval)validateattributes(desiredpval, ...
    0342                     {'numeric'}, ...
    0343                     {'scalar', 'positive', '>=', 0, '<=', 1}));
    0344             p.addParamValue('tail','both', ...
    0345                 @(tail)any([strcmp(tail,'both'), strcmp(tail,'left'),strcmp(tail,'right')]));
    0346             p.addParamValue('name','noname',@(name)ischar(name)); % dealt with in the main constructor
    0347             p.addParamValue('targVal',0.0, ...
    0348                 @(targVal)validateattributes(targVal, ...
    0349                     {'numeric'}, ...
    0350                     {'scalar', '>=', -1, '<=', 1}));
    0351             p.parse(sosObj,varargin{:});
    0352
    0353
    0354             sample1 = p.Results.sample1;
    0355             sample2 = p.Results.sample2;
    0356             s1ColName = p.Results.s1ColName; %#ok<PROP>
    0357             s2ColName = p.Results.s2ColName; %#ok<PROP>
    0358             obj.desiredpvalCondition = p.Results.desiredpvalCondition;
    0359
    0360
    0361             if strcmp(obj.desiredpvalCondition,'N/A')
    0362                 obj.desiredpvalConditionHandle = @sosCorrelTest.returnNaN;
    0363             elseif strcmp(obj.desiredpvalCondition,'<=')
    0364                 obj.desiredpvalConditionHandle = @le;
    0365             elseif strcmp(obj.desiredpvalCondition,'=>')
    0366                 obj.desiredpvalConditionHandle = @ge;
    0367             end
    0368
    0369             % can't have a desiredpval without a desiredpvalcondition
    0370             if strcmp(obj.desiredpvalCondition,'N/A')
    0371                 obj.desiredpval = NaN;
    0372             else
    0373                 obj.desiredpval = p.Results.desiredpval;
    0374             end
    0375
    0376
    0377             % perform additional checks on the objects
    0378
    0379             present1 = sosObj.containsSample(sample1);
    0380             if (present1 == 0 )
    0381                 error('sos sosObject does not contain sample1');
    0382             end
    0383
    0384             present2 = sosObj.containsSample(sample2);
    0385             if (present2 == 0 )
    0386                 error('sos sosObject does not contain sample2');
    0387             end
    0388
    0389             col1 = sample1.colName2colNum(s1ColName); %#ok<PROP>
    0390             if(col1 == -1)
    0391                 error('<s1ColName> not a column of data in <sample1>');
    0392             end
    0393
    0394             col2 = sample1.colName2colNum(s2ColName); %#ok<PROP>
    0395             if(col2 == -1)
    0396                 error('<s2ColName> not a column of data in <sample2>');
    0397             end
    0398
    0399             if isempty(sample1.data)
    0400                 error('sample 1 does not contain items - did you fill it yet?');
    0401             end
    0402
    0403             if isempty(sample2.data)
    0404                 error('sample 2 does not contain items - did you fill it yet?');
    0405             end
    0406             %all variables check out, create the stats test
    0407
    0408
    0409             if (length(sample1.data{col1}) ~= length(sample2.data{col2}))
    0410                 error('Sample 1 and Sample 2 must have the same number of observations for a paired comparison');
    0411             end
    0412
    0413             if sample1.n < 4
    0414                 % see note in the calculation for why the following is the
    0415                 % case
    0416                 error('The sosCorrelTest requires at least 4 items per sample');
    0417             end
    0418
    0419
    0420             obj.s1 = sample1;
    0421             obj.s2 = sample2;
    0422             obj.s1ColName = s1ColName; %#ok<PROP>
    0423             obj.s2ColName = s2ColName; %#ok<PROP>
    0424             obj.s1Col = col1;
    0425             obj.s2Col = col2;
    0426             obj.type = p.Results.type;
    0427             obj.tail = p.Results.tail;
    0428             obj.targVal = p.Results.targVal;
    0429
    0430             obj.runSpecificTest = @obj.runMatchCorrelztest;
    0431
    0432             obj.label = [obj.s1.name, '{',obj.s1ColName, '}-', ...
    0433                     obj.s2.name, '{',obj.s2ColName, '}', ...
    0434                     ':z[',obj.type,']'];
    0435
    0436
    0437         end
    0438
    0439
    0440     end
    0441
    0442
    0443     methods (Static)
    0444
    0445         %% userHypothesis = returnNaN(~,~) STATIC FUNCTION
    0446         function userHypothesis = returnNaN(~,~)
    0447             % returns NaN
    0448             userHypothesis = NaN;
    0449         end
    0450
    0451         %%  flag = validTestType(str) STATIC FUNCTION
    0452         function flag = validTestType(str)
    0453             % returns 1 if the name of type of test is 'matchCorrel'; error otherwise.
    0454
    0455             flag = 0; %#ok<NASGU>
    0456
    0457             if(ischar(str) == false)
    0458                 error('<Type> must be "matchCorrel"');
    0459             end
    0460
    0461             if (strcmp(str,'matchCorrel'))
    0462                 flag = 1;
    0463             else
    0464                 error('<Type> must be "matchCorrel"');
    0465             end
    0466
    0467         end
    0468     end
    0469
    0470 end
    0471