%%%THIS FILE GRAPHS PARAMETER SPACE FROM TWIN DATA SETS
%%%WHENEVER CV(DZ)<1/2*CV(MZ)
%
% All parameter sets within this space are mathematically equally likely
% but are not necessarily biologically equally likely
%
% By Matt Keller
% Nov 26, 2004
% 
% For more explanation, see:
% Keller & Coventry (2005). Quantifying and Addressing Parameter
% Indeterminacy in Twin Studies. Twin Research and Human Genetics, 8(3), pp??
%
% function [] = ParSpace2(cvmz, cvdz, rmin, n5,c5,a6,n6,n7,c7,a8,n8)
%
%Where:
% cvmz = covariance bw MZ twins
% cvdz = covariance bw DZ twins
% rmin = minimum r (correlation between DZ twins for non-additive genetic effect)
% n5 = non-additive parameter on row 5 (see Table I in paper cited above)
% c5 = common environment parameter on row 5 (maximum C parameter value)
% a6 = additive genetic parameter on row 6 (maximum A parameter value)
% (and so forth...)
%
% To get Figure II from paper, you would type: ParSpace2(.47,.17,0,.29,.17,.35,.12,.39,.08,.225,.242)
%
% A figure should be created on your desktop
% This function also saves a jpg file called "ParSpace1.jpg" in your default Matlab folder
% as well as a enhanced meta file "ParSpace1.emf" because the jpg file is
% sometimes messed up


function [] = ParSpace2(cvmz,cvdz,rmin, n5,c5,a6,n6,n7,c7,a8,n8)

%Create four important points on x-axis
x7=0-n7;
x5=0-n5;
x8=a8-n8;
x6=a6-n6;

%Create values for setting x and y axis lengths
yspace=c5*.01;
xspace=(x6+n7)*.01;


%find the slope and intercept for back border
beta1=(c7-c5)/(x7-x5);
beta0=c7-beta1*x7;

%Create 1000 evenly spaced data pairs for x 
xhat2=linspace(x7,x5,1000);
xhat1=linspace(x8,x6,1000);

%Create 1000 evenly spaced data pairs for y 
yhat2=beta1*xhat2+beta0;
yhat1=zeros(1,1000);

%Create color for lines depending on r
shade=linspace(0,.8,1000);


%These commands keeps the figure window the same dimensions and coming up
%in same place
set(gcf,'Position',[250,170,700,560])
set(gca,'FontSize',13)
set(gca,'YLim', [0 c5+10*yspace])
set(gca,'XLim', [x7-10*xspace x6+25*xspace])


%Create the space
hold on
for t = 1:1000
plot([xhat2(1,t),xhat1(1,t)],[yhat2(1,t),yhat1(1,t)],'-','LineWidth',2,'Color',[shade(1,t),shade(1,t),shade(1,t)])
end
plot([x7-.001,x5-.001],[c7,c5],'-','LineWidth',3,'Color',[1,1,1])


%labels
sxlabel('\16\times\i V\Hat_A - V\Hat_{NA}')
sylabel('\16\times\i V\Hat_C')
stitle(['\14\times Space of Mathematically Equally Likely Parameters Given {\i{CV}\Hat_{MZ}} = ', num2str(cvmz,'%3.3g'), ' {\i{CV}\Hat_{DZ}} = ', num2str(cvdz,'%3.3g')])

%create dots on boundaries
text(x8,0,'\oplus','FontSize',20,'Color',[0,0,0],'HorizontalAlignment','center','VerticalAlignment','middle','Fontweight','bold')
text(x6,0,'\bullet','FontSize',30,'Color',[0,0,0],'HorizontalAlignment','center','VerticalAlignment','middle')
text(x5,c5,'\bullet','FontSize',30,'Color',[0,0,0],'HorizontalAlignment','center','VerticalAlignment','middle')
text(x7,c7,'\bullet','FontSize',30,'Color',[0,0,0],'HorizontalAlignment','center','VerticalAlignment','middle')

%print the boundaries within brackets
text(x8-.03,.004,['(',num2str(a8*100,'%2.2g'),', ',num2str(n8*100,'%2.2g'),', 0)  '],'FontSize',11,'HorizontalAlignment','right','Fontweight','bold')
text(x6+.01,.004,['(',num2str(a6*100,'%2.2g'),', ',num2str(n6*100,'%2.2g'),', 0)  '],'FontSize',11,'HorizontalAlignment','left','Fontweight','bold')
text(x5+.01,c5+.004,['(0, ',num2str(n5*100,'%2.2g'),', ',num2str(c5*100,'%2.2g'),')  '],'FontSize',11,'HorizontalAlignment','left','Fontweight','bold')
text(x7+.01,c7,['(0, ',num2str(n7*100,'%2.2g'),', ',num2str(c7*100,'%2.2g'),')  '],'FontSize',11,'HorizontalAlignment','left','Color',[1,1,1],'Fontweight','bold')

%Make box for r (grey to black)
g=[x6-25*xspace,x6-35*xspace];
f=linspace(50*yspace,75*yspace,51);
shade2=linspace(0,.8,50);
for i=1:50
    fill([g(1),g(1),g(2),g(2)],[f(i),f(i+1),f(i+1),f(i)],[shade2(i),shade2(i),shade2(i)],'EdgeColor','none')
end


%Place legends for r and parameter estimates
stext(x6-23*xspace,50*yspace,'\12\times {\i r\Hat} = 0.25','HorizontalAlignment','left')
stext(x6-23*xspace,75*yspace,['\12\times {\i r\Hat} = ',num2str(rmin,'%3.3g')],'HorizontalAlignment','left')

text(x6-30*xspace,c5,'\fontsize{30}\bullet\fontsize{20}\bf\oplus','Color',[0,0,0],'HorizontalAlignment','right','VerticalAlignment','middle')
stext(x6-28*xspace,c5,'\11\times Boundaries ({\i V\Hat_A , V\Hat_{NA} , V\Hat_C}) \mult 100','HorizontalAlignment','left','VerticalAlignment','top')

text(x6-30*xspace,(9*c5)/10,'\oplus','FontSize',20,'Color',[0,0,0],'HorizontalAlignment','right','VerticalAlignment','middle','Fontweight','bold')
text(x6-28*xspace,(9*c5)/10,'Classical Twin Design Estimates','FontSize',11,'Fontname','times','HorizontalAlignment','left')
hold off


print -djpeg ParSpace2               %This saves a jpg file called "ParSpace1.jpg" in your default Matlab folder
print -dmeta ParSpace2