/* Computations to support banking a time series to 45 degrees. See Wicklin (2016), "Banking to 45 degrees: Aspect ratios for time series plots" The DO Loop blog, published 20JAN2016 http://blogs.sas.com/content/iml/2016/01/20/banking-to-45-aspect-ratio-time-series.html */ /* Data sets from Cleveland (1993) Visualizing Data */ /* Download from http://www.ats.ucla.edu/stat/r/examples/vizdata/ */ libname vizdata "C:\Users\frwick\Documents\My SAS Files\vizdata"; proc iml; /* function that returns the median absolute slope for a time series with points (x_i, y_i). The cooerdinates are passed in as vectors. */ start MedianAbsSlope(x, y); dx = dif(x) / range(x); /* scaled change in horizontal directions */ dy = dif(y) / range(y); /* scaled change in vertical directions */ m = median( abs(dy/dx) ); /* median of absolute slopes */ return( 1/m ); /* value that makes median equal unity */ finish; /* Compute aspect ratios suggested by Cleveland for time series plots. Three computational methods for "banking to 45 degrees" are: method="ms": median slope method="ao": average orientation method="awo": average weighted orientation */ /* in SAS, aspect ratio is height/width = (vertical units)/(horizontal units) */ start AspectRatio(x, y, method="ms", interval={0.01 100}) global(g_dx, g_dy); if upcase(method)="MS" then return MedianAbsSlope(x, y); /* The prefix "g_" indicates global variables used for finding root */ g_dx = dif(x) / range(x); /* scaled change in horizontal directions */ g_dy = dif(y) / range(y); /* scaled change in vertical directions */ if upcase(method)="AO" then return froot("MeanAbsOrientF", interval); else if upcase(method)="AWO" then do; return froot("MeanWtAbsOrientF", interval); end; return ( . ); /* invalid method */ finish; /* Function that returns the difference between pi/4 (45 degrees) and the weighted mean of abs(angles) for line segments. Find the aspect ratio that makes this function zero. */ start MeanWtAbsOrientF(a) global(g_dx, g_dy); theta = abs(atan(a * g_dy / g_dx)); /* ABS angles, given aspect ratio */ length = sqrt(g_dx##2 + a##2 * g_dy##2); m = sum(theta#length) / sum(length); /* weighted mean abs angle */ return(m - constant("pi")/4); /* target value = pi/4 */ finish; /* Function that returns the difference between pi/4 (45 degrees) and the mean of abs(angles) for line segments. Find the aspect ratio that makes this function zero. */ start MeanAbsOrientF(a) global(g_dx, g_dy); theta = abs(atan(a * g_dy / g_dx)); /* ABS angles, given aspect ratio */ m = mean(theta); /* mean abs angle */ return(m - constant("pi")/4); /* target value = pi/4 */ finish; /************************************************************************/ /* Melanoma Data, see Cleveland (1993, p. 158) */ /************************************************************************/ title "Melanoma Data"; use vizdata.Melanoma; read all var "Year" into x; read all var "Incidence" into y; close; dx = dif(x) / range(x); /* scaled change in horizontal directions */ dy = dif(y) / range(y); /* scaled change in vertical directions */ m = median( abs(dy/dx) ); /* median of absolute slopes */ print m[L="Median Slope"] (1/m)[L="Aspect Ratio"]; /* create a series plot */ aspect = round(1/m, 0.1); call symputx("ratio", aspect); title2 "Aspect Ratio = &ratio"; %let width = 500; ods graphics / width=&width.px height=%sysevalf(50+ &width*&ratio)px; call series(x, y) label={"Year" "Incidence"} procopt=("aspect="+char(aspect)) option="markers"; /*--------------------------------------------------------------------------*/ /* look at distribution of absolute slopes */ AbsSlopes = abs(dy/dx); title2 "Distribution of Absolute Slopes"; call histogram(AbsSlopes) other="refline 2.7 / axis=x;"; create absslopes from absslopes[c="AbsSlopes"]; append from absslopes; close; submit median=m; proc univariate data=absslopes; histogram / href=&median odstitle="Distribution of Absolute Slopes"; cdfplot / vref=50 odstitle="Cumulative Distribution of Absolute Slopes"; run; endsubmit; /*--------------------------------------------------------------------------*/ /* compute aspect ratios for three methods */ m = j(3, 2); m[1,1] = AspectRatio(x, y, "ms"); m[2,1] = AspectRatio(x, y, "ao"); m[3,1] = AspectRatio(x, y, "awo"); m[,2] = 1 / m[,1]; print m[L="Banking Methods (Melanoma Data)" c = {"Aspect Ratio" "1/(Aspect Ratio)"} r={"Median Slope" "Average Orientation" "Average Weighted-Orient"}]; aspect = round(m[3,1], 0.001); call symputx("ratio", aspect); title2 "Aspect Ratio = &ratio"; %let width = 640; ods graphics / width=&width.px height=%sysevalf(100 + &width*&ratio)px; call series(x, y) procopt=("aspect="+char(aspect)) option="markers"; /************************************************************************/ /* CO2 measurements, see Cleveland (1993, p. 159) and Herr and Agrawala */ /************************************************************************/ title "CO2 Measurements"; use vizdata.carbon; read all var "Positions" into x; read all var "Data" into y; close; m = j(3, 2); m[1,1] = AspectRatio(x, y, "ms"); m[2,1] = AspectRatio(x, y, "ao"); m[3,1] = AspectRatio(x, y, "awo"); m[,2] = 1 / m[,1]; print m[L="Banking Methods (Cleveland)" c = {"Aspect Ratio" "1/(Aspect Ratio)"} r={"Median Slope" "Average Orient" "Average Weighted-Orient"}]; aspect = round(m[3,1], 0.001); call symputx("ratio", aspect); title2 "Aspect Ratio = &ratio"; %let width = 640; ods graphics / width=&width.px height=%sysevalf(100 + &width*&ratio)px; call series(x, y) procopt=("aspect="+char(aspect));