sosCorrelTest

 - 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.

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