* ttest.equal_variances.gs
* Description:
* This script plots t-test values and shades areas that are significant at
* a specified level. This script assumes that the variance of the two samples
* is the same. If you are uncertain about this assumption, you should
* use a f-test to determine if you can make this assumption.
* For step-by-step instructions on how to run this script,
* see http://ccr.meteor.wisc.edu/model/grads/grads_statsig.html
****************************************************************************
* This script needs to be used in concert with the web-page specified below.
* This site allows you to compute cut-off t-values, meaning that a t-value
* BELOW the cutoff at a given grid cell means that the difference between the
* two means is NOT significant and a t-test score ABOVE the value means that
* the difference IS significant.
* I.e., if ttest (grid cell) > ttest (cutoff), then significant
* If ttest(grid cell) < ttest (cutoff), then NOT significant
*
* TO CALCULATE CUTOFF T-VALUE, goto:
* http://members.aol.com/johnp71/pdfs.html
* USE "Student t"
* ENTER 'n' DEGREES OF FREEDOM (= N1+N2-2) AND SIGNIFICANCE LEVEL 'p'
* 0.01 = 99% LEVEL, 0.05 = 95% LEVEL, 0.10 = 90% LEVEL
* THEN CLICK 'Calc t' BUTTON. The cutoff 't' value will appear in the box
* in the 't' box.
*****************************************************************************
* OPEN DATA FILES (PUT PATH AND FILE NAMES HERE)
sdfopen /grove4/selin/foam_indyr/ha.F1.6k.djfi486605.TS1.nc
sdfopen /grove4/selin/foam_indyr/ha.F1.0k.djfi486605.TS1.nc
* ENABLE OUTPUT FILE
enable print sig
set lon -180 180
* ENTER THE SAMPLE SIZES OF YOUR TWO SAMPLES BELOW
*n1, n2 are number of observations for each experiment
*degrees of freedom = n1 + n2 -2
n1 = 120
n2 = 120
* In the formulas for x1, x2, s1, and s2, you must set the time
* increment for each of the formulas according to the sample size.
* x = average, s = standard deviation
x1 = ave(ts1.1,t=1,t=120)
s1 = sqrt(ave(pow(ts1.1-x1,2),t=1,t=120)*(n1/(n1-1)))
x2 = ave(ts1.2,t=1,t=120)
s2 = sqrt(ave(pow(ts1.2-x2,2),t=1,t=120)*(n2/(n2-1)))
* compute t test statistic
ssqr = ( (n1-1)*s1*s1 + (n2-1)*s2*s2 )/(n1 + n2 - 2)
denom = sqrt(ssqr*((1.0/n1) + (1.0/n2)))
num = x1 - x2
ttest = num/denom
* Take absolute value for a two-tailed test
ttest = abs(ttest)
* Use url above to calculate the cutoff t.
* In this case, at the .01 (99%) significance level and df=238,
* the cutoff t value is 2.597
set grid off
set grads off
set gxout grfill
* set clevs equal to the cutoff t so that all grid cells with values above
* the cutoff-t will be shaded grey (they are signif)
* and plot the significant areas
set clevs 2.597
set ccols 0 15
d ttest
* run the color bar script
run cbar.gs
* contour difference of two runs over significance
set gxout contour
set cint 0.25
d x1-x2
draw title FOAM 6K-0K TS1 DJF Delta (contour), Sig at 99% (shaded)
print