*! Date : 26 Oct 2011 *! Version : 1.00 *! Authors : Wenfeng Gong; Jan Osterman (Duke Center for Health Policy) *! Email : gongwenf@gmail.com *! Description : Proportional Venn Diagrams /* Future ideas 1. treat problem when a population is involved in the overlap part of the other two populations. eg. pvenn pA(45) pB(25) pC(50) pAB(18) pAC(24) pBC(3) ang(120) 2. fill color in circles. */ * capture program drop pvenn pr pvenn version 11.0 preserve set more off /* data input */ syntax [varlist(default=none)] [, plabel(namelist min=1 max=3) /// pA(real 0) pB(real 0) pC(real 0) pAB(real 0) pAC(real 0) pBC(real 0) /// ptotal(real 0) ang(real 0) pos(integer 1) drawtotal(integer 1) drawlabel(integer 1) /// lc(string) lp(string) lw(string) n(integer 20000) *] local gopt "`options'" *di "options: ""`options'" local popA=`pA' local popB=`pB' local popC=`pC' local popAB=`pAB' local popAC=`pAC' local popBC=`pBC' local angle=`ang' local popt=`ptotal' local nameA="Population A" local nameB="Population B" local nameC="Population C" local i 0 foreach nam of local plabel { if `++i'==1 local nameA "`nam'" if `i'==2 local nameB "`nam'" if `i'==3 local nameC "`nam'" } local i 0 foreach v of local varlist { confirm numeric variable `v' local vtab`++i'= "`v'" *di "vtab`i'= " "`vtab`i''" } if `i'>3 { di "too many variables selected" exit } else{ if `i'>=1 { tab `vtab1', matcell(cellvals) * matrix list cellvals * di cellvals[2,1] local popA=cellvals[2,1] local popt=`r(N)' local nameA "`vtab1'" } if `i'>=2 { tab `vtab1' `vtab2', matcell(cellvals) local popB=cellvals[1,2]+cellvals[2,2] local popAB=cellvals[2,2] local nameB "`vtab2'" } if `i'==3{ tab `vtab1' `vtab3', matcell(cellvals) local popC=cellvals[1,2]+cellvals[2,2] local popAC=cellvals[2,2] local nameC "`vtab3'" tab `vtab2' `vtab3', matcell(cellvals) local popBC=cellvals[2,2] } } display "`nameA': " "`popA'" display "`nameB': " "`popB'" display "`nameC': " "`popC'" display "AB Overlap: " "`popAB'" display "AC Overlap: " "`popAC'" display "BC Overlap: " "`popBC'" display "Total population: " "`popt'" display "Draw total population rectangle (0,not; 1, yes): " "`drawtotal'" display "Rotation angle of circle B: " "`angle'" display "Relative position of circle C (1,clockwise; -1,counterclockwise): " "`pos'" if `drawtotal'==0 local popt=0 /* Set up line color, pattern, width */ if "`lc'"=="" local lc="cranberry emerald brown lavender" if "`lw'"=="" local lw="0.8 0.8 0.8 0.8" if "`lp'"=="" local lp="solid solid solid solid" forvalues z=1/4 { local xopt`z'=" n(`n') `gopt'" local conlc`z': word `z' of `lc' *di "conlc`z' " "`conlc`z''" if "`conlc`z''"~="" local xopt`z' "`xopt`z'' lc(`conlc`z'')" else local xopt`z' "`xopt`z'' lc(`conlc1')" local conlp`z': word `z' of `lp' if "`conlp`z''"~="" local xopt`z' "`xopt`z'' lp(`conlp`z'')" else local xopt`z' "`xopt`z'' lp(`conlp1')" local conlw`z': word `z' of `lw' if "`conlw`z''"~="" local xopt`z' "`xopt`z'' lw(`conlw`z'')" else local xopt`z' "`xopt`z'' lw(`conlw1')" *di "xopt`z'= " "`xopt`z''" } /* Adjust scale of parameters */ local r1= 1 local r2= ((`popB'/`popA')^(0.5))* `r1' local r3= ((`popC'/`popA')^(0.5))* `r1' local p1= _pi * `r1' * `r1' local p2= (`popB'/`popA') * `p1' local p3= (`popC'/`popA') * `p1' local p12= (`popAB'/`popA') * `p1' local p13= (`popAC'/`popA') * `p1' local p23= (`popBC'/`popA') * `p1' local popt= (`popt'/`popA') * `p1' if `popA'<=0 { local r1= 0 local r2= 1 local r3= ((`popC'/`popB')^(0.5))* `r2' local p1= 0 local p2= _pi * `r2' * `r2' local p3= (`popC'/`popB') * `p2' local p12= (`popAB'/`popB') * `p2' local p13= (`popAC'/`popB') * `p2' local p23= (`popBC'/`popB') * `p2' local popt= (`popt'/`popB') * `p2' } if `popA'<=0 & `popB'<=0 { local r1=0 local r2=0 local r3=1 local p1= 0 local p2= 0 local p3= _pi * `r3' * `r3' local p12= (`popAB'/`popC') * `p3' local p13= (`popAC'/`popC') * `p3' local p23= (`popBC'/`popC') * `p3' local popt= (`popt'/`popC') * `p3' } *di "r1= " "`r1'" *di "r2= " "`r2'" *di "r3= " "`r3'" *di "p1= " "`p1'" *di "p2= " "`p2'" *di "p3= " "`p3'" *di "p12= " "`p12'" *di "p13= " "`p13'" *di "p23= " "`p23'" /* Estimate a1 a2*/ local a1=`r1' local a2=`r2' local tp12=0 if `p12'>0 { if `r1'>= `r2' { while (abs(`tp12'-`p12')/`p12')> 0.01 { local a2= `a2'-`r2'/1000 local a1= (`r1'^2-`r2'^2+`a2'^2)^(0.5) if `a2'< ((-1)*`r2') { di "cannot find appropriate location for Circle B" exit } else { local tp12 = (`r1'^2)*(acos(`a1'/`r1'))- `a1'*((`r1'^2-`a1'^2)^(0.5)) + (`r2'^2)*(acos(`a2'/`r2'))- `a2' * ((`r2'^2-`a2'^2)^(0.5)) } } } else{ while (abs(`tp12'-`p12')/`p12')> 0.01 { local a1= `a1'-`r1'/1000 local a2= (`r2'^2-`r1'^2+`a1'^2)^(0.5) if `a1'< ((-1)*`r1') { di "cannot find appropriate location for Circle A" exit } else { local tp12 = (`r1'^2)*(acos(`a1'/`r1'))- `a1'*((`r1'^2-`a1'^2)^(0.5)) + (`r2'^2)*(acos(`a2'/`r2'))- `a2' * ((`r2'^2-`a2'^2)^(0.5)) } } } } *di "a1= " "`a1'" *di "a2= " "`a2'" *di "tp12= " "`tp12'" if `r3'>0 & `r3'!=. { /* Estimate a13 a31*/ local a13=`r1' local a31=`r3' local tp13=0 if `p13'>0 { if `r1'>= `r3' { while (abs(`tp13'-`p13')/`p13')> 0.01 { local a31= `a31'-`r3'/1000 local a13= (`r1'^2-`r3'^2+`a31'^2)^(0.5) if `a31'< ((-1)*`r3') { di "cannot find appropriate location for Circle C" exit } else { local tp13 = (`r1'^2)*(acos(`a13'/`r1'))- `a13'*((`r1'^2-`a13'^2)^(0.5)) + (`r3'^2)*(acos(`a31'/`r3'))- `a31' * ((`r3'^2-`a31'^2)^(0.5)) } } } else{ while (abs(`tp13'-`p13')/`p13')> 0.01 { local a13= `a13'-`r1'/1000 local a31= (`r3'^2-`r1'^2+`a13'^2)^(0.5) if `a13'< ((-1)*`r1') { di "cannot find appropriate location for Circle C" exit } else { local tp13 = (`r1'^2)*(acos(`a13'/`r1'))- `a13'*((`r1'^2-`a13'^2)^(0.5)) + (`r3'^2)*(acos(`a31'/`r3'))- `a31' * ((`r3'^2-`a31'^2)^(0.5)) } } } } *di "a13= " "`a13'" *di "a31= " "`a31'" *di "tp13= " "`tp13'" /* Estimate a23 a32*/ local a23=`r2' local a32=`r3' local tp23=0 if `p23'>0 { if `r2'>= `r3' { while (abs(`tp23'-`p23')/`p23')> 0.01 { local a32= `a32'-`r3'/1000 local a23= (`r2'^2-`r3'^2+`a32'^2)^(0.5) if `a32'< ((-1)*`r3') { di "cannot find appropriate location for Circle C" exit } else { local tp23 = (`r2'^2)*(acos(`a23'/`r2'))- `a23'*((`r2'^2-`a23'^2)^(0.5)) + (`r3'^2)*(acos(`a32'/`r3'))- `a32' * ((`r3'^2-`a32'^2)^(0.5)) } } } else{ while (abs(`tp23'-`p23')/`p23')> 0.01 { local a23= `a23'-`r2'/1000 local a32= (`r3'^2-`r2'^2+`a23'^2)^(0.5) if `a23'< ((-1)*`r2') { di "cannot find appropriate location for Circle C" exit } else { local tp23 = (`r2'^2)*(acos(`a23'/`r2'))- `a23'*((`r2'^2-`a23'^2)^(0.5)) + (`r3'^2)*(acos(`a32'/`r3'))- `a32' * ((`r3'^2-`a32'^2)^(0.5)) } } } } *di "a23= " "`a23'" *di "a32= " "`a32'" *di "tp23= " "`tp23'" } /* locate the center of circles */ local cx1=0 local cy1=0 local cx2=0 local cy2=0 local cx3=0 local cy3=0 local l12=0 local l13=0 local l23=0 local angle213=0 local l12= `a1'+`a2' local cx2= `cx1'+(`l12') * cos(_pi*`angle'/180) local cy2= `cy1'+(`l12') * sin(_pi*`angle'/180) if `r3'>0 & `r3'!=.{ local l13= `a13'+`a31' local l23= `a23'+`a32' if `l12'*`l13'!=0 local angle213= acos((`l12'^2+`l13'^2-`l23'^2)/(2*`l12'*`l13')) local angle2=_pi*`angle'/180-`pos'*`angle213' local cx3=`cx1'+(`l13') * cos(`angle2') local cy3=`cy1'+(`l13') * sin(`angle2') } *di "`l12' `l13' `l23'" *di "angle213= " "`angle213'" *di "angle2= " "`angle2'" *di "cx3= " "`cx3'" *di "cy3= " "`cy3'" /* locate the place for label */ local lnameA="" local lnameB="" local lnameC="" local xlabel1=`cx1' local ylabel1=`cy1' local xlabel2=`cx2' local ylabel2=`cy2' local xlabel3=`cx3' local ylabel3=`cy3' if `drawlabel'==1 { if `r1'>0 local lnameA="`nameA'" if `r2'>0 local lnameB="`nameB'" if `r3'>0 local lnameC="`nameC'" local ll12=`l12'-`r2' local ll13=`l13'-`r3' local ll21=`l12'-`r1' local ll23=`l23'-`r3' local ll31=`l13'-`r1' local ll32=`l23'-`r2' local angles1=(_pi*`angle'/180-`pos'*`angle213'/2) local temp= min(`ll12',`ll13',0)*cos(`angles1') if `temp'!=0 local xlabel1=`cx1'+ `temp' +0.3 *`temp'/abs(`temp') local temp= min(`ll12',`ll13',0)*sin(`angles1') if `temp'!=0 local ylabel1=`cy1'+ `temp' +0.3 *`temp'/abs(`temp') local angle123=0 if `l12'*`l23'!=0 local angle123= acos((`l12'^2+`l23'^2-`l13'^2)/(2*`l12'*`l23')) local angles2=(_pi*`angle'/180+`pos'*`angle123'/2) local temp= min(`ll21',`ll23',0)*cos(`angles2') if `temp'!=0 local xlabel2=`cx2'- `temp' -0.3 *`temp'/abs(`temp') local temp= min(`ll21',`ll23',0)*sin(`angles2') if `temp'!=0 local ylabel2=`cy2'- `temp' -0.3 *`temp'/abs(`temp') local angle132=0 if `l13'*`l23'!=0 local angle132= acos((`l13'^2+`l23'^2-`l12'^2)/(2*`l13'*`l23')) local angles3=_pi*`angle'/180+`pos'*(`angle213'+ `angle132'/2) local temp= min(`ll31',`ll32',0)*cos(`angles3') if `temp'!=0 local xlabel3=`cx3'+ `temp' +0.3 *`temp'/abs(`temp') local temp= min(`ll31',`ll32',0)*sin(`angles3') if `temp'!=0 local ylabel3=`cy3'+ `temp' +0.3 *`temp'/abs(`temp') *di "`ll12',`ll13'" *di "angles1:" "`angles1'" *di "xlabel1:" "`xlabel1'" *di "`ll21',`ll23'" qui local angles2=180* `angles2'/_pi *di "angles2:" "`angles2'" *di "xlabel2:" "`xlabel2'" *di "xlabel3:" "`xlabel3'" } /*draw the circles*/ foreach circles in 1 2 3{ local rangel`circles'= (`cx`circles''-`r`circles'') local rangeh`circles'= (`cx`circles''+`r`circles'') *local erangel`circles'=floor(2*`rangel`circles'') /2 *local erangeh`circles'=ceil(2*`rangeh`circles'') /2 } local minX=min(`rangel1',`rangel2', `rangel3') local maxX=max(`rangeh1',`rangeh2', `rangeh3') *di "minX= " "`minX'" *di "maxX= " "`maxX'" local Xrange=`maxX'-`minX' foreach circles in 1 2 3{ local rangeb`circles'= (`cy`circles''-`r`circles'') local ranget`circles'= (`cy`circles''+`r`circles'') *local erangeb`circles'=floor(2*`rangeb`circles'')/2 *local eranget`circles'=ceil(2*`ranget`circles'') /2 } local minY=min(`rangeb1',`rangeb2', `rangeb3') local maxY=max(`ranget1',`ranget2', `ranget3') *di "minY= " "`minY'" *di "maxY= " "`maxY'" local Yrange=`maxY'-`minY' local aspectratio=`Yrange'/`Xrange' local rangearea=`Xrange'*`Yrange' *di "rangearea: `rangearea'" if `rangearea'<=`popt' & `popt'>0 { local totalX=sqrt(`popt'/`rangearea')*`Xrange' local totalY=sqrt(`popt'/`rangearea')*`Yrange' } else { if `popt'>0 di "Margin (i.e. subjects outside all the population A, B, and C) is not large enough; total population cannot be represented as outside rectangle." local totalX=`Xrange' local totalY=`Yrange' } local totalminX=`minX'-(`totalX'-`Xrange')/2 local totalmaxX=`maxX'+(`totalX'-`Xrange')/2 local totalminY=`minY'-(`totalY'-`Yrange')/2 local totalmaxY=`maxY'+(`totalY'-`Yrange')/2 local minX=floor(2*`totalminX')/2 local maxX=ceil(2*`totalmaxX')/2 local minY=floor(2*`totalminY')/2 local maxY=ceil(2*`totalmaxY')/2 local ylabel4=`totalmaxY'-0.1 local xlabel4=`totalminX'+0.1 *di "aspectratio= " "`aspectratio'" *di "rangel1= " "`rangel1'" *di "rangel2= " "`rangel2'" *di "rangel3= " "`rangel3'" if `rangearea'<=`popt' { if `r3'>0 & `r3'!=.{ twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' "`lnameA'") xlabel(none) ylabel(none) /// xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio') legend(order(1 "`nameA'" 3 "`nameB'" 5 "`nameC'" 7 "Total Population")) /// || function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') `xopt1' /// || function y = sqrt(`r2'^2 - (x-`cx2')^2)+`cy2',range(`rangel2' `rangeh2') text(`ylabel2' `xlabel2' "`lnameB'") `xopt2' /// || function y = -sqrt(`r2'^2 - (x-`cx2')^2)+`cy2', range(`rangel2' `rangeh2') `xopt2' /// || function y = sqrt(`r3'^2 - (x-`cx3')^2)+`cy3',range(`rangel3' `rangeh3') text(`ylabel3' `xlabel3' "`lnameC'") `xopt3' /// || function y = -sqrt(`r3'^2 - (x-`cx3')^2)+`cy3',range(`rangel3' `rangeh3') `xopt3' /// || function y = `totalminX',range(`totalminY' `totalmaxY') horizontal xtitle("") ytitle("") `xopt4' /// || function y = `totalmaxX',range(`totalminY' `totalmaxY') horizontal `xopt4' /// || function y = `totalminY',range(`totalminX' `totalmaxX') `xopt4' text(`ylabel4' `xlabel4' "Total", place(se)) /// || function y = `totalmaxY',range(`totalminX' `totalmaxX') `xopt4' } else { if `r2'>0 & `r2'!=.{ twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' "`lnameA'") xlabel(none) ylabel(none) /// xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio') legend(order(1 "`nameA'" 3 "`nameB'" 5 "Total Population")) /// || function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') `xopt1' /// || function y = sqrt(`r2'^2 - (x-`cx2')^2)+`cy2',range(`rangel2' `rangeh2') text(`ylabel2' `xlabel2' "`lnameB'") `xopt2' /// || function y = -sqrt(`r2'^2 - (x-`cx2')^2)+`cy2', range(`rangel2' `rangeh2') `xopt2' /// || function y = `totalminX',range(`totalminY' `totalmaxY') xtitle("") ytitle("") horizontal `xopt4' /// || function y = `totalmaxX',range(`totalminY' `totalmaxY') horizontal `xopt4' /// || function y = `totalminY',range(`totalminX' `totalmaxX') `xopt4' text(`ylabel4' `xlabel4' "Total", place(se)) /// || function y = `totalmaxY',range(`totalminX' `totalmaxX') `xopt4' } else { twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' "`lnameA'") xlabel(none) ylabel(none) /// xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio') legend(order(1 "`nameA'" 3 "Total Population")) /// || function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') `xopt1' /// || function y = `totalminX',range(`totalminY' `totalmaxY') xtitle("") ytitle("") horizontal `xopt4' /// || function y = `totalmaxX',range(`totalminY' `totalmaxY') horizontal `xopt4' /// || function y = `totalminY',range(`totalminX' `totalmaxX') `xopt4' text(`ylabel4' `xlabel4' "Total", place(se)) /// || function y = `totalmaxY',range(`totalminX' `totalmaxX') `xopt4' } } exit } if `r3'>0 & `r3'!=.{ twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' "`lnameA'") xlabel(none) ylabel(none) /// xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio') legend(order(1 "`nameA'" 3 "`nameB'" 5 "`nameC'" )) /// || function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') `xopt1' /// || function y = sqrt(`r2'^2 - (x-`cx2')^2)+`cy2',range(`rangel2' `rangeh2') text(`ylabel2' `xlabel2' "`lnameB'") `xopt2' /// || function y = -sqrt(`r2'^2 - (x-`cx2')^2)+`cy2', range(`rangel2' `rangeh2') `xopt2' /// || function y = sqrt(`r3'^2 - (x-`cx3')^2)+`cy3',range(`rangel3' `rangeh3') text(`ylabel3' `xlabel3' "`lnameC'") `xopt3' /// || function y = -sqrt(`r3'^2 - (x-`cx3')^2)+`cy3',range(`rangel3' `rangeh3') `xopt3' } else { if `r2'>0 & `r2'!=.{ twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' "`lnameA'") xlabel(none) ylabel(none) /// xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio') legend(order(1 "`nameA'" 3 "`nameB'" )) /// || function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') `xopt1' /// || function y = sqrt(`r2'^2 - (x-`cx2')^2)+`cy2',range(`rangel2' `rangeh2') text(`ylabel2' `xlabel2' "`lnameB'") `xopt2' /// || function y = -sqrt(`r2'^2 - (x-`cx2')^2)+`cy2', range(`rangel2' `rangeh2') `xopt2' } else { twoway function y = sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') text(`ylabel1' `xlabel1' "`lnameA'") xlabel(none) ylabel(none) /// xtitle("") ytitle("") `xopt1' aspectratio(`aspectratio') legend(order(1 "`nameA'" )) /// || function y = -sqrt(`r1'^2 - (x-`cx1')^2)+`cy1', range(`rangel1' `rangeh1') `xopt1' } } end