%%%THIS FILE GRAPHS PARAMETER SPACE FROM TWIN DATA SETS
%%%WHENEVER CV(DZ)/CV(MZ)>1/2
%
% 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 [] = ParSpace1(cvmz, cvdz, rmin, n1,c1,a2,c2,n3,c3,a4,c4)
%
%Where:
% cvmz = covariance bw MZ twins
% cvdz = covariance bw DZ twins
% rmin = minimum r (correlation between DZ twins for non-additive genetic effect)
% n1 = non-additive parameter on row 1 (see Table I in paper cited above)
% c1 = common environment parameter on row 1 (maximum C parameter value)
% a2 = additive genetic parameter on row 2 (maximum A parameter value)
% (and so forth...)
%
% To get Figure III from paper, you would type: ParSpace1(.42,.28,0,.14,.28,.28,.14,.187,.233,.28,.14)
%
% 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 [] = ParSpace1(cvmz, cvdz, rmin, n1,c1,a2,c2,n3,c3,a4,c4)

%Create three important points on x-axis
x3=0-n3;
x1=0-n1;
x2=a2;


%find the slope and intercept for back border
beta1=(c3-c1)/(x3-x1);
beta0=c3-beta1*x3;

%Create 1000 evenly spaced data pairs for x 
xhat2=linspace(x3,x1,1000);
xhat1=zeros(1,1000)+x2;

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

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

%Create values for setting x and y axis lengths
yspace=(c1-c2)*.01;
xspace=(a4+n3)*.01;

%These commands keeps the figure window the same dimensions and coming up
%in same place, and also set the x & y axis lengths
set(gcf,'Position',[250,170,700,560])
set(gca,'FontSize',13)
set(gca,'YLim', [c2-10*yspace c1+10*yspace])
set(gca,'XLim', [x3-10*xspace x2+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([x3-.001,x1-.001],[c3,c1],'-','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(x2,c2,'\oplus','FontSize',20,'Color',[0,0,0],'HorizontalAlignment','center','VerticalAlignment','middle','Fontweight','bold')
text(x3,c3,'\bullet','FontSize',30,'Color',[0,0,0],'HorizontalAlignment','center','VerticalAlignment','middle')
text(x1,c1,'\bullet','FontSize',30,'Color',[0,0,0],'HorizontalAlignment','center','VerticalAlignment','middle')

%print the boundaries within brackets
text(x2+.01,c2+.004,['(',num2str(a2*100,'%2.2g'), ', 0, ',num2str(c2*100,'%2.2g'),')  '],'FontSize',11,'HorizontalAlignment','left','Fontweight','bold')
text(x1+.01,c1+.004,['(0, ',num2str(n1*100,'%2.2g'),', ',num2str(c1*100,'%2.2g'),')  '],'FontSize',11,'HorizontalAlignment','left','Fontweight','bold')
text(x3+.01,c3,['(0, ',num2str(n3*100,'%2.2g'),', ',num2str(c3*100,'%2.2g'),')  '],'FontSize',11,'HorizontalAlignment','left','Color',[1,1,1],'Fontweight','bold')

%Make box for r (grey to black)
g=[x2-20*xspace,x2-30*xspace];
f=linspace((c2+50*yspace),(c2+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(x2-18*xspace,c2+50*yspace,'\12\times {\i r\Hat} = 0.25','HorizontalAlignment','left')
stext(x2-18*xspace,c2+75*yspace,['\12\times {\i r\Hat} = ',num2str(rmin,'%3.3g')],'HorizontalAlignment','left')

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

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

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