data cholesterol (index=(deathcause));
	set sashelp.heart;
	where deathcause ne '' and cholesterol ne .;
  	keep deathcause cholesterol;   
run;

*calculate cumulative percents by cause of death;
proc rank data = cholesterol out = chol_percent percent; *groups=101;
	by deathcause;
	ranks chol_percent;
	var cholesterol;
run;

* transform percents into groups;
data chol_grouped (index=(deathcole=(deathcause cholesterol)));
	set chol_percent; *chol_percent_sorted;	
	select ;
		when (25 < chol_percent <= 75) group = 50; /*central 50%*/
		when (5 < chol_percent <= 95) group = 90; /*central 90%*/
		otherwise group = 100; /*central 100%*/
	end;
run;


* output mean and median (and other basic measures) in long format;
proc univariate data=cholesterol ; *noprint;
by deathcause;
var cholesterol;
*output out = chol_sumstat 	p5=chol_p05 q1=chol_p25 median=chol_p50 q3=chol_p75 p95=chol_p95 mean=chol_mean ;
*ods output quantiles = chol_quant (rename=(estimate=cholesterol) index=(deathcol=(deathcause cholesterol)));
ods output BasicMeasures = chol_MLoc (rename=(LocValue=cholesterol) 
									drop=VarName VarMeasure VarValue
									index=(deathco=(deathcause cholesterol)));
run;


* kde stands for kernel desnity estimation and the procedure does just that;
proc kde data = cholesterol;
by deathcause;
univar cholesterol / NGRID=401 unistats percentiles plots=none /*(plots=(density HISTDENSITY)*/
     out = chol_dens (rename=(value=cholesterol) drop=var );
run;

* assign to each density the group that corresponds to its percentile; 
* the code’s done in such a way that the extraneous, i.e. non-density, rows are purged;
* add in 'markers' for mean and median;
data chol_den_6 (where = (density ne . or LocMeasure ne '') drop = group chol_percent density2 mirror2);
  	merge chol_dens chol_grouped chol_mloc (where=(LocMeasure IN('Mean','Median')));
  	by deathcause cholesterol;
	mirror=-density;	
	if first.deathcause then group2=group;	
	retain group2;
	group2=coalesce(group,group2);	/* coalesce returns (first) non-missing value*/
	density2=lag(density);
	mirror2=lag(mirror);
	if LocMeasure = 'Mean' then do; density_mean = density2; mirror_mean = mirror2; end;
		else if LocMeasure = 'Median' then do; density_median = density2; mirror_median = mirror2; end;
run;

* assign shades of blue to groups, i.e., assign cycle of group colors;
proc template;
   define Style styles.violin; 
      parent = styles.listing;
      style GraphColors from graphcolors /
           "gdata1" = cxBDD7E7 /* lighest shade*/
           "gdata2" = cx6BAED6
		   "gdata3" = cx3182BD;	/* darkest shade */	
   end;
run;

* set path, dpi, and style;
%let gpath='C:\';
%let dpi=200;
ods listing gpath=&gpath image_dpi=&dpi style=styles.violin;


* create violin plot with shaded quantile bounds and markers for median and mean;
ods graphics / width=4in height=6in imagename='Violin_Horz_Group';
title 'Cholesterol Densities by Death Cause';
title2 'Violin plot with shaded quantile bounds'; 
proc sgpanel data=chol_den_6 nocycleattrs noautolegend;
	panelby deathcause / layout=rowlattice ROWHEADERPOS= LEFT onepanel novarname noborder;
  	band x=cholesterol upper=density lower=mirror / group=group2 fill outline
    lineattrs=(pattern=solid color=black) ; *transparency=0.5;
	highlow x=cholesterol high=density_median low=mirror_median / LINEATTRS=(color=cx08519C);
  	highlow x=cholesterol high=density_mean low=mirror_mean / LINEATTRS=(color=cxDD1C77);  	
  	colaxis label='Cholesterol' refticks; *grid;
  	rowaxis display=none; *(nolabel);
run;