main.js (5471B)
1 /** 2 * @license Apache-2.0 3 * 4 * Copyright (c) 2018 The Stdlib Authors. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 19 'use strict'; 20 21 // MODULES // 22 23 var isNumberArray = require( '@stdlib/assert/is-number-array' ).primitives; 24 var isTypedArrayLike = require( '@stdlib/assert/is-typed-array-like' ); 25 var setReadOnly = require( '@stdlib/utils/define-read-only-property' ); 26 var quantileFactory = require( './../../base/dists/normal/quantile' ).factory; 27 var cdfFactory = require( './../../base/dists/normal/cdf' ).factory; 28 var atanh = require( '@stdlib/math/base/special/atanh' ); 29 var tanh = require( '@stdlib/math/base/special/tanh' ); 30 var tCDF = require( './../../base/dists/t/cdf' ); 31 var sqrt = require( '@stdlib/math/base/special/sqrt' ); 32 var min = require( '@stdlib/math/base/special/min' ); 33 var print = require( './print.js' ); // eslint-disable-line stdlib/no-redeclare 34 var pcorr = require( './pcorr.js' ); 35 var validate = require( './validate.js' ); 36 37 38 // VARIABLES // 39 40 var normQuantile = quantileFactory( 0.0, 1.0 ); 41 var normCDF = cdfFactory( 0.0, 1.0 ); 42 43 44 // MAIN // 45 46 /** 47 * Computes a Pearson product-moment correlation test between paired samples. 48 * 49 * @param {NumericArray} x - first data array 50 * @param {NumericArray} y - second data array 51 * @param {Options} [options] - function options 52 * @param {number} [options.alpha=0.05] - significance level 53 * @param {string} [options.alternative='two-sided'] - alternative hypothesis (`two-sided`, `less` or `greater`) 54 * @param {number} [options.rho=0.0] - correlation under H0 55 * @throws {TypeError} x argument has to be a typed array or array of numbers 56 * @throws {TypeError} y argument has to be a typed array or array of numbers 57 * @throws {Error} x and y must be arrays of the same length 58 * @throws {Error} x and y must contain at least four elements 59 * @throws {TypeError} options has to be simple object 60 * @throws {TypeError} must provide valid options 61 * @returns {Object} test result object 62 * 63 * @example 64 * var x = [ 2, 4, 3, 1, 2, 3 ]; 65 * var y = [ 3, 2, 4, 1, 2, 4 ]; 66 * var out = pcorrTest( x, y ); 67 */ 68 function pcorrTest( x, y, options ) { 69 var method; 70 var alpha; 71 var cint; 72 var opts; 73 var pval; 74 var stat; 75 var alt; 76 var err; 77 var out; 78 var rho; 79 var val; 80 var df; 81 var se; 82 var n; 83 var r; 84 var z; 85 86 if ( !isTypedArrayLike( x ) && !isNumberArray( x ) ) { 87 throw new TypeError( 'invalid argument. First argument `x` must be a numeric array. Value: `' + x + '`.' ); 88 } 89 if ( !isTypedArrayLike( y ) && !isNumberArray( y ) ) { 90 throw new TypeError( 'invalid argument. Second argument `y` must be a numeric array. Value: `' + y + '`.' ); 91 } 92 n = x.length; 93 if ( n !== y.length ) { 94 throw new Error( 'invalid arguments. Arguments `x` and `y` must be arrays of the same length' ); 95 } 96 opts = {}; 97 if ( options ) { 98 err = validate( opts, options ); 99 if ( err ) { 100 throw err; 101 } 102 } 103 if ( opts.alpha === void 0 ) { 104 alpha = 0.05; 105 } else { 106 alpha = opts.alpha; 107 } 108 if ( n < 4 ) { 109 throw new Error( 'not enough observations. `x` and `y` must contain at least four observations.' ); 110 } 111 if ( opts.rho === void 0 ) { 112 rho = 0.0; 113 } else { 114 rho = opts.rho; 115 } 116 if ( opts.alternative === void 0 ) { 117 alt = 'two-sided'; 118 } else { 119 alt = opts.alternative; 120 } 121 122 r = pcorr( x, y ); 123 z = atanh( r ); 124 se = 1.0 / sqrt( n - 3 ); 125 if ( rho === 0.0 ) { 126 // Use t-test for H0: rho = 0.0 vs H1: rho != 0.0... 127 method = 't-test for Pearson correlation coefficient'; 128 df = n - 2; 129 stat = sqrt( df ) * r / sqrt( 1.0 - (r*r) ); 130 switch ( alt ) { 131 case 'greater': 132 pval = 1.0 - tCDF( stat, df ); 133 break; 134 case 'less': 135 pval = tCDF( stat, df); 136 break; 137 case 'two-sided': 138 default: 139 pval = 2.0 * min( tCDF( stat, df), 1.0 - tCDF( stat, df )); 140 break; 141 } 142 } else { 143 // Use large-sample normality to calculate p-value based on Fisher's z transform... 144 method = 'Fisher\'s z transform test for Pearson correlation coefficient'; 145 stat = ( z - atanh( rho ) ) * sqrt( n - 3 ); 146 switch ( alt ) { 147 case 'greater': 148 pval = normCDF( -stat ); 149 break; 150 case 'less': 151 pval = 1.0 - normCDF( -stat ); 152 break; 153 case 'two-sided': 154 default: 155 pval = 2.0 * min( normCDF( -stat ), 1.0 - normCDF( -stat )); 156 break; 157 } 158 } 159 160 switch ( alt ) { 161 case 'greater': 162 cint = [ tanh( z - ( se*normQuantile( 1.0 - alpha ) ) ), 1.0 ]; 163 break; 164 case 'less': 165 cint = [ -1.0, tanh( z + ( se*normQuantile( 1.0 - alpha ) ) ) ]; 166 break; 167 case 'two-sided': 168 default: 169 val = se * normQuantile( 1.0 - ( alpha/2.0 ) ); 170 cint = [ tanh( z - val ), tanh( z + val ) ]; 171 break; 172 } 173 174 out = {}; 175 setReadOnly( out, 'rejected', pval <= alpha ); 176 setReadOnly( out, 'alpha', alpha ); 177 setReadOnly( out, 'pValue', pval ); 178 setReadOnly( out, 'statistic', stat ); 179 setReadOnly( out, 'ci', cint ); 180 setReadOnly( out, 'alternative', alt ); 181 setReadOnly( out, 'method', method ); 182 setReadOnly( out, 'nullValue', rho ); 183 setReadOnly( out, 'pcorr', r ); 184 setReadOnly( out, 'print', print ); 185 return out; 186 } 187 188 189 // EXPORTS // 190 191 module.exports = pcorrTest;