Skip to content

Instantly share code, notes, and snippets.

@mnarayan

mnarayan/ParSpace1.m

Last active November 3, 2020 18:28
Show Gist options
  • Save mnarayan/dfca55584b930891b24bda94da2932db to your computer and use it in GitHub Desktop.
Save mnarayan/dfca55584b930891b24bda94da2932db to your computer and use it in GitHub Desktop.
Download ZIP
Scripts for "Quantifying and addressing parameter indeterminacy in the classical twin design."
Raw
ParSpace1.m
%%%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
Raw
ParSpace2.m
%%%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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment