Spearman Rank Correlation Trend Evaluation in EnviroInsite

The Spearman Rank Correlation trend test source code in EnviroInsite is a C# adaptation of the C++ source code provided with the Numerical Recipes reference book. The Microsoft .NET code was used for sorting and the critical values of s_r for sample size of 60 and below were taken from an online source at http://www.york.ac.uk/depts/maths/tables/spearman.pdf. For sample sizes greater than 60, the student's t distribution described of the critical value was calculated as implemented in the Numerical Recipes code.

public string SpearmanRankCorrelationEval(int PercentSig, List<DateTime> dts, List<double> vals)

{

if (dts.Count <= 3 || (PercentSig == 1 && dts.Count <= 4))

return "Inconclusive";

List<double> dts_tck = new List<double>(dts.Count);

int i = 0;

foreach (DateTime dt in dts)

{

dts_tck.Add( dt.Ticks );

i++;

}

double d = 0;

double zd = 0;

double probd = 0;

double rs = 0;

double probrs = 0;

spear(ref dts_tck, ref vals, ref d, ref zd, ref probd, ref rs, ref probrs);

double dPercentSig = (double)PercentSig/100;

int n = dts.Count;

if (n > 60)

{

if (probrs < dPercentSig)

{

if (rs < 0)

return "Negative Trend";

else

return "Positive Trend";

}

else

return "Inconclusive";

}

else

{

double[] s_cr = null;

if (PercentSig == 10)

s_cr = new double[]{ -1, -1, -1, 1.0, 0.8, 0.657, 0.571, 0.524, 0.483, 0.455,

0.427, 0.406, 0.385, 0.367, 0.354, 0.341, 0.328, 0.317, 0.309, 0.299,

0.292, 0.284, 0.278, 0.271, 0.265, 0.259, 0.255, 0.250, 0.245, 0.240,

0.236, 0.232, 0.229, 0.225, 0.222, 0.219, 0.216, 0.212, 0.210, 0.207,

0.204, 0.202, 0.199, 0.197, 0.194, 0.192, 0.190, 0.188, 0.186, 0.184,

0.182, 0.180, 0.179, 0.177, 0.175, 0.174, 0.172, 0.171, 0.169, 0.168};

else if(PercentSig == 5)

s_cr = new double[]{ -1, -1, -1, 1.0, 0.90, 0.829, 0.714, 0.643, 0.600, 0.564,

0.536, 0.503, 0.484, 0.464, 0.443, 0.429, 0.414, 0.401, 0.391, 0.380,

0.370, 0.361, 0.353, 0.344, 0.337, 0.331, 0.324, 0.317, 0.312, 0.306,

0.301, 0.298, 0.291, 0.287, 0.283, 0.279, 0.275, 0.271, 0.267, 0.264,

0.261, 0.257, 0.254, 0.251, 0.248, 0.246, 0.243, 0.240, 0.238, 0.235,

0.233, 0.231, 0.228, 0.226, 0.224, 0.222, 0.220, 0.218, 0.216, 0.214};

else if (PercentSig == 1)

s_cr = new double[]{ -1, -1, -1, -1, 1.0, 0.943, 0.893, 0.833, 0.783, 0.745,

0.709, 0.671, 0.648, 0.622, 0.604, 0.582, 0.566, 0.550, 0.535, 0.520,

0.508, 0.496, 0.486, 0.476, 0.466, 0.457, 0.448, 0.440, 0.433, 0.425,

0.418, 0.412, 0.405, 0.399, 0.394, 0.388, 0.383, 0.378, 0.373, 0.368,

0.364, 0.359, 0.355, 0.351, 0.347, 0.343, 0.340, 0.336, 0.333, 0.329,

0.326, 0.323, 0.320, 0.317, 0.314, 0.311, 0.308, 0.306, 0.303, 0.301};

if ( s_cr[n-1] < 0 || Math.Abs(rs) < s_cr[n-1])

{

return "Inconclusive";

}

else

{

if (rs > 0)

return "Positive Trend";

else

return "Negative Trend";

}

void spear(ref List<double> data1, ref List<double> data2, ref double d,

ref double zd, ref double probd, ref double rs, ref double probrs)

{

/* Given two data arrays, data1[0..n-1] and data2[0..n-1], this routine returns their

sum squared difference of ranks as D, the number of standard deviations by which D

deviates from its null-hypothesis expected value as zd, the two-side p-value of this

deviation as probd, Spearmans rank correlation as rs, and the two-sided p-value of its

deviation from zero as probrs. The external routines crank and sort2 are used. A small

value of either probd or proprs indicates a significant correlation (rs positive) or

anticorrelation (rs negative)*/

double bet = 0;

int j;

int n = data1.Count;

double vard, t, fac, en3n, en, df, aved;

double sf = 0;

double sg = 0;

List<double> wksp1 = new List<double>();

List<double> wksp2 = new List<double>();

for (j = 0; j<n; j++)

{

wksp1.Add( data1[j] );

wksp2.Add( data2[j] );

}

sort2(ref wksp1, ref wksp2);

crank(ref wksp1, ref sf);

sort2(ref wksp2, ref wksp1);

crank(ref wksp2, ref sg);

d = 0;

for(j = 0; j<n; j++)

d += Math.Pow(wksp1[j]-wksp2[j],2);

en = n;

en3n = en*en*en-en;

aved = en3n / 6.0 - (sf+ sg)/12;

fac = (1.0 - sf/en3n) * (1.0-sg/en3n);

if (fac > 0.0)

{

vard = ((en - 1) * en * en * Math.Pow(en + 1, 2) / 36) * fac;

zd = (d - aved) / Math.Sqrt(vard);

probd = erfcc(Math.Abs(zd) / 1.4142136);

rs = (1 - (6 / en3n) * (d + (sf + sg) / 12)) / Math.Sqrt(fac);

fac = (rs + 1) * (1 - rs);

t = rs * Math.Sqrt((en - 2) / fac);

df = en - 2;

probrs = betai(0.5 * df, 0.5, df / (df + t * t));

}

else

{

rs = 0;

probrs = 0.0;

}

void crank( ref List<double> w, ref double s)

// Given a sorted array w[0..n-1], replaces the elements by their rank, including midranking

// of ties and returns as s the sum of f^3-f, where f is the number of element in each tie.

{

int j=1;

int ji, jt;

double t, rank;

int n = w.Count;

s = 0;

while(j < n)

{

if(w[j]!= w[j-1])

{

w[j-1]=j;

++j;

}

else

{

for(jt=j+1; jt <= n && w[jt-1] == w[j-1]; jt++);

rank = 0.5 * (j+jt-1);

for(ji=j; ji<=(jt-1); ji++)

w[ji-1] = rank;

t = jt - j;

s += (t*t*t-t);

j = jt;

}

if (j == n) w[n - 1] = n;

}

void sort2( ref List<double> arr, ref List<double> brr)

{

List<pr> prs = new List<pr>(arr.Count);

int i = 0;

foreach(double x in arr)

{

pr p = new pr();

p.x_ = x;

p.y_ = brr[i];

prs.Add(p);

i++;

}

prs.Sort();

i = 0;

foreach(pr p in prs)

{

arr[i] = p.x_;

brr[i] = p.y_;

i++;

}