rawTable.sas

/****************************************************************

* Program Name: macro.raw.table.sas *

* Program dir : /udd/nhron/macrodevelopment/vitc/ *

* Program Date: Mar., 2008 *

* *

* Note : *

* If the input data has missing values in *

* one of the visits, the pair of data will *

* deleted for the analaysis *

* Programmed by: Marion McPhee and Rong Chen *

****************************************************************/

%macro rawTable(dsn=, /* input data */

q =, /* surrogate measure of intake */

r =, /* reference measure of intake */

m = /* biomarker */

);

data main;

set &dsn end=eof ;

if nmiss(&q, &m, &r)^=0 then do;

id1=id;

nmis+1;

end;

if eof then do;

call symput('nmis', trim(left(nmis)));

end;

run;

proc sort data=main;

by id descending id1;

run;

data main;

set main end=eof;

by id descending id1;

retain id2;

if first.id then id2=id1;

if id2=id then delete;

if eof and &nmis>0 then do;

%put Note: &nmis pair of records are deleted because of missing value of the factors in at least one of them;

end;

drop id1 id2;

/* change coding of visit */

visit=visit-1;

run;

proc sort data=main;

by id visit;

run;

proc transpose data=main out=temp1(drop=_NAME_) prefix=&q;

by id;

var &q;

run;

proc transpose data=main out=temp2(drop=_NAME_) prefix=&r;

by id;

var &r;

run;

proc transpose data=main out=temp3(drop=_NAME_) prefix=&m;

by id;

var &m;

run;

data temp;

merge temp1 temp2 temp3;

by id;

run;

/* obtain total# of observations per visit */

data _null_;

set temp end=eof;

if eof then do;

call symput('totn', trim(left(_n_)));

end;

run;

data raw1(drop=xbar1 xbar2 xbar3 xdbar1-xdbar3 xbarj11-xbarj13

xbarj21-xbarj23 sumx11-sumx13 sumx21-sumx23)

raw2(keep=xbar1-xbar3 sumxi1-sumxi3 sumxdi1-sumxdi3 xdbar1-xdbar3

xbarj11-xbarj13 xbarj21-xbarj23);

set temp end=eof;

retain sumxi1 sumxi2 sumxi3 sumxdi1-sumxdi3 sumx11-sumx13 sumx21-sumx23 nn 0;

array a(3) &q.1 &r.1 &m.1;

array b(3) &q.2 &r.2 &m.2;

array d(3) qi ri mi;

array sumxi(3) sumxi1 sumxi2 sumxi3;

array xbari(3) xbari1 xbari2 xbari3;

array sumxdi(3) sumxdi1-sumxdi3;

array sumx1(3) sumx11 sumx12 sumx13;

array sumx2(3) sumx21 sumx22 sumx23;

n=&totn;

nn=nn+1;

one=1;

do k=1 to 3;

*a(k)=a(k)/100;

*b(k)=b(k)/100;

d(k)=b(k)-a(k);

xbari(k)=((a(k)+b(k))/2);

sumxi(k)=sumxi(k)+((a(k)+b(k))/2);

sumxdi(k)=sumxdi(k)+d(k);

sumx1(k)=sumx1(k)+a(k);

sumx2(k)=sumx2(k)+b(k);

end;

output raw1;

if eof then do;

array xbar(3) xbar1 xbar2 xbar3;

array xdbar(3) xdbar1-xdbar3;

array xbarj1(3) xbarj11 xbarj12 xbarj13;

array xbarj2(3) xbarj21 xbarj22 xbarj23;

do k=1 to 3;

xbar(k)=sumxi(k)/n;

xdbar(k)=sumxdi(k)/n;

xbarj1(k)=sumx1(k)/n;

xbarj2(k)=sumx2(k)/n;

end;

output raw2;

end;

run;

proc corr cov data=raw1;

var &q.1 &q.2 &r.1 &r.2;

run;

data parta

partb;

set raw1 end=eof;

retain xbar1 xbar2 xbar3 xdbar1-xdbar3 xbarj11-xbarj13 xbarj21-xbarj23;

retain cov_a_xbar1-cov_a_xbar54 cov_b_xbar1-cov_b_xbar54

cov_c_xbar1-cov_c_xbar54 varxbar1-varxbar6

cov_xbark_xbarl1-cov_xbark_xbarl27 0;

if _n_=1 then set raw2;

retain varx1 varx2 varx3 sumaa1-sumaa9 sumbb1-sumbb9 cova1-cova9 0;

retain sumaaa1-sumaaa27 covb1-covb9 sumbbb1-sumbbb27 sumx1-sumx27

sumc11-sumc19 sumc21-sumc29 sumca1-sumca9 sumccc1-sumccc27

sumcc1-sumcc81 covc1-covc9 sumc101-sumc109 sumc201-sumc209

covca1-covca9 covcb1-covcb9 sumcca1-sumcca81 sumccb1-sumccb81

sumaaaa1-sumaaaa81 sumbbbb1-sumbbbb81 sumxx1-sumxx81

varxd1-varxd3 covc121-covc129 0;

array a(3) &q.1 &r.1 &m.1;

array b(3) &q.2 &r.2 &m.2;

array xbarj1(3) xbarj11 xbarj12 xbarj13;

array xbarj2(3) xbarj21 xbarj22 xbarj23;

array xbari(3) xbari1 xbari2 xbari3;

array xbar(3) xbar1 xbar2 xbar3;

array xdbar(3) xdbar1 xdbar2 xdbar3;

array sumaa(3,3) sumaa1-sumaa9;

array varx(3) varx1 varx2 varx3;

array vara(3,3) vara1-vara9;

array varb(3,3) varb1-varb9;

array varxbar(3,2) varxbar1-varxbar6;

array cova(3,3) cova1-cova9;

array covb(3,3) covb1-covb9;

array covc(3,3) covc1-covc9;

array covca(3,3) covca1-covca9;

array covcb(3,3) covcb1-covcb9;

array sumaaa(3,3,3) sumaaa1-sumaaa27;

array covaa(3,3,3) covaa1-covaa27;

array sumbbb(3,3,3) sumbbb1-sumbbb27;

array covbb(3,3,3) covbb1-covbb27;

array covcc(3,3,3,3) covcc1-covcc81;

array covcca(3,3,3,3) covcca1-covcca81;

array covccb(3,3,3,3) covccb1-covccb81;

array d(3) &q &r &m;

array sumbb(3,3) sumbb1-sumbb9;

array varxd(3) varxd1-varxd3;

array sumx(3,3,3) sumx1-sumx27;

array sumc1(3,3) sumc11-sumc19;

array sumc2(3,3) sumc21-sumc29;

array sumc(3,3) sumca1-sumca9;

array sumc10(3,3) sumc101-sumc109;

array sumc20(3,3) sumc201-sumc209;

array sumc0(3,3) sumca01-sumca09;

array sumcc(3,3,3,3) sumcc1-sumcc81;

array sumcca(3,3,3,3) sumcca1-sumcca81;

array sumccb(3,3,3,3) sumccb1-sumccb81;

array covac(3,3,3) covac1-covac27;

array sumaaaa(3,3,3,3) sumaaaa1-sumaaaa81;

array sumbbbb(3,3,3,3) sumbbbb1-sumbbbb81;

array sumxx(3,3,3,3) sumxx1-sumxx81;

array covaaaa(3,3,3,3) covaaaa1-covaaaa81;

array covbbbb(3,3,3,3) covbbbb1-covbbbb81;

array covaacc(3,3,3,3) covaacc1-covaacc81;

array cova_sigmat(3,3) cova_sigmat1-cova_sigmat9;

array covb_sigmas(3,3) covb_sigmas1-covb_sigmas9;

array covb_sigmad(3,3) covb_sigmad1-covb_sigmad9;

array covb_sigmat(3,3) covb_sigmat1-covb_sigmat9;

array covc_sigmat(3,3) covc_sigmat1-covc_sigmat9;

array cov_a_xbar(3,3,3,2) cov_a_xbar1-cov_a_xbar54;

array cov_b_xbar(3,3,3,2) cov_b_xbar1-cov_b_xbar54;

array cov_c_xbar(3,3,3,2) cov_c_xbar1-cov_c_xbar54;

array cov_xbark_xbarl(3,3,2) cov_xbark_xbarl1-cov_xbark_xbarl18;

array cov_sigmasqs_xbar(3,2) cov_sigmasqs_xbar1-cov_sigmasqs_xbar6;

array cov_sigmasqd_xbar(3,2) cov_sigmasqd_xbar1-cov_sigmasqd_xbar6;

array cov_sigmasqt_xbar(3,2) cov_sigmasqt_xbar1-cov_sigmasqt_xbar6;

array cov_betaq_xbar(3,2) cov_betaq_xbar1-cov_betaq_xbar6;

array var_betaq_xbar(2) var_betaq_xbar1 var_betaq_xbar2;

array cov_qbarj_betaqrbarj(2) cov_qbarj_betaqrbarj1 cov_qbarj_betaqrbarj2;

array var_alphaq(2) var_alphaq1 var_alphaq2;

array var_betam_xbar(2) var_betam_xbar1 var_betam_xbar2;

array cov_mbarj_betamrbarj(2) cov_mbarj_betamrbarj1 cov_mbarj_betamrbarj2;

array var_alpham(2) var_alpham1 var_alpham2;

array sumc12a(3,3) sumc12a1-sumc12a9;

array sumc12b(3,3) sumc12b1-sumc12b9;

array sumc12c(3,3) sumc12c1-sumc12c9;

array covc12(3,3) covc121-covc129;

do k=1 to 3;

varx(k)=varx(k)+(xbari(k)-xbar(k))**2/n;

varxd(k)=varxd(k)+(d(k)-xdbar(k))**2/n;

/* A22.1 */

varxbar(k,1)=varxbar(k,1)+(a(k)-xbarj1(k))**2/n**2;

varxbar(k,2)=varxbar(k,2)+(b(k)-xbarj2(k))**2/n**2;

do l1=1 to 3;

sumaa(k,l1)=sumaa(k,l1)+((xbari(k)-xbar(k))**2*(xbari(l1)-xbar(l1))**2)/n;

sumbb(k,l1)=sumbb(k,l1)+(d(k)-xdbar(k))**2*(d(l1)-xdbar(l1))**2/n;

sumc1(k,l1)=(a(k)-xbarj1(k))*(a(l1)-xbarj1(l1));

sumc2(k,l1)=(b(k)-xbarj2(k))*(b(l1)-xbarj2(l1));

sumc(k,l1)=sumc1(k,l1)+sumc2(k,l1);

sumc12a(k,l1)=(a(k)-xbarj1(k))*(b(l1)-xbarj2(l1));

sumc12b(k,l1)=(a(l1)-xbarj1(l1))*(b(k)-xbarj2(k));

sumc12c(k,l1)=sumc12a(k,l1)+sumc12b(k,l1);

covc12(k,l1)=covc12(k,l1)+sumc12c(k,l1)/(2*n);

cova(k,l1)=cova(k,l1)+((xbari(k)-xbar(k))*(xbari(l1)-xbar(l1)))/n;

covb(k,l1)=covb(k,l1)+((d(k)-xdbar(k))*(d(l1)-xdbar(l1)))/n;

covca(k,l1)=covca(k,l1)+sumc1(k,l1)/n;

covcb(k,l1)=covcb(k,l1)+sumc2(k,l1)/n;

do l2=1 to 3;

sumaaa(k,l1,l2)=sumaaa(k,l1,l2)+((xbari(k)-xbar(k))**2*(xbari(l1)-xbar(l1))*

(xbari(l2)-xbar(l2)))/n;

sumbbb(k,l1,l2)=sumbbb(k,l1,l2)+((d(k)-xdbar(k))**2*(d(l1)-xdbar(l1))

*(d(l2)-xdbar(l2)))/n;

sumc1(k,l2)=(a(k)-xbarj1(k))*(a(l2)-xbarj1(l2));

sumc2(k,l2)=(b(k)-xbarj2(k))*(b(l2)-xbarj2(l2));

sumc(k,l2)=sumc1(k,l2)+sumc2(k,l2);

*sumcc(k,l1,l2)=sumcc(k,l1,l2)+sumc(k,l1)*sumc(k,l2)/(4*n);

sumx(k,l1,l2)=sumx(k,l1,l2)+((xbari(k)-xbar(k))*(xbari(l1)-xbar(l1)))

*sumc(k,l2)/(2*n);

end;

end;

end;

do k=1 to 3;

do l=1 to 3;

sumc10(k,l)=(a(k)-xbarj1(k))*(a(l)-xbarj1(l));

sumc20(k,l)=(b(k)-xbarj2(k))*(b(l)-xbarj2(l));

sumc0(k,l)=sumc10(k,l)+sumc20(k,l);

end;

end;

do k1=1 to 3;

do l1=1 to 3;

do k2=1 to 3;

do l2=1 to 3;

sumcca(k1,l1,k2,l2)=sumcca(k1,l1,k2,l2)+(sumc10(k1,l1)*sumc10(k2,l2))/n;

sumccb(k1,l1,k2,l2)=sumccb(k1,l1,k2,l2)+(sumc20(k1,l1)*sumc20(k2,l2))/n;

sumcc(k1,l1,k2,l2)=sumcc(k1,l1,k2,l2)+sumc0(k1,l1)*sumc0(k2,l2)/(4*n);

sumaaaa(k1,l1,k2,l2)=sumaaaa(k1,l1,k2,l2)+(xbari(k1)-xbar(k1))*

(xbari(l1)-xbar(l1))*(xbari(k2)-xbar(k2))*(xbari(l2)-xbar(l2))/n;

sumbbbb(k1,l1,k2,l2)=sumbbbb(k1,l1,k2,l2)+(d(k1)-xdbar(k1))

*(d(l1)-xdbar(l1))*(d(k2)-xdbar(k2))*(d(l2)-xdbar(l2))/n;

sumxx(k1,l1,k2,l2)=sumxx(k1,l1,k2,l2)+(xbari(k1)-xbar(k1))*

(xbari(l1)-xbar(l1))*sumc(k2,l2)/(2*n);

end;

end;

end;

end;

do k=1 to 3;

do l=1 to 3;

/* A23.13 */

cov_xbark_xbarl(k,l,1)=cov_xbark_xbarl(k,l,1)+(a(k)-xbarj1(k))

*(a(l)-xbarj1(l))/n**2;

cov_xbark_xbarl(k,l,2)=cov_xbark_xbarl(k,l,2)+(b(k)-xbarj2(k))

*(b(l)-xbarj2(l))/n**2;

do m=1 to 3;

/* A23.4 */

cov_a_xbar(k,l,m,1)=cov_a_xbar(k,l,m,1)+(xbari(k)-xbar(k))

*(xbari(l)-xbar(l))*(a(m)-xbarj1(m))/n;

cov_a_xbar(k,l,m,2)=cov_a_xbar(k,l,m,2)+(xbari(k)-xbar(k))

*(xbari(l)-xbar(l))*(b(m)-xbarj2(m))/n;

* put k= l= m= xbari(k)= xbar(k)= xbari(l)= xbar(l)= a(m)=

* xbarj1(m)= cov_a_xbar(k,l,m,1)= cov_a_xbar(k,l,m,2)= ;

/* A23.5 */

cov_b_xbar(k,l,m,1)=cov_b_xbar(k,l,m,1)+(d(k)-xdbar(k))*

(d(l)-xdbar(l))*(a(m)-xbarj1(m))/n;

cov_b_xbar(k,l,m,2)=cov_b_xbar(k,l,m,2)+(d(k)-xdbar(k))*

(d(l)-xdbar(l))*(b(m)-xbarj2(m))/n;

/* a23.6 */

cov_c_xbar(k,l,m,1)=cov_c_xbar(k,l,m,1)+((a(k)-xbarj1(k))*

(a(l)-xbarj1(l))+(b(k)-xbarj2(k))*(b(l)-xbarj2(l)))

*(a(m)-xbarj1(m))/(2*n);

cov_c_xbar(k,l,m,2)=cov_c_xbar(k,l,m,2)+((a(k)-xbarj1(k))*

(a(l)-xbarj1(l))+(b(k)-xbarj2(k))*(b(l)-xbarj2(l)))

*(b(m)-xbarj2(m))/(2*n);

end;

end;

end;

output parta;

if eof then do;

do k=1 to 3;

do l=1 to 3;

do m=1 to 3;

do j=1 to 2;

cov_a_xbar(k,l,m,j)=cov_a_xbar(k,l,m,j);

cov_b_xbar(k,l,m,j)=cov_b_xbar(k,l,m,j);

cov_c_xbar(k,l,m,j)=cov_c_xbar(k,l,m,j);

end;

end;

end;

end;

do k=1 to 3;

do l1=1 to 3;

vara(k,l1)=1/n*(sumaa(k,l1)-cova(k,l1)**2);

varb(k,l1)=1/n*(sumbb(k,l1)-covb(k,l1)**2);

covc(k,l1)=(covca(k,l1)+covcb(k,l1))/2;

do l2=1 to 3;

covaa(k,l1,l2)=(sumaaa(k,l1,l2)-cova(k,l1)*cova(k,l2))/n;

covbb(k,l1,l2)=(sumbbb(k,l1,l2)-covb(k,l1)*covb(k,l2))/n;

covac(k,l1,l2)=(sumx(k,l1,l2)-cova(k,l1)*covc(k,l2))/n;

end;

end;

end;

do k1=1 to 3;

do l1=1 to 3;

do k2=1 to 3;

do l2=1 to 3;

covcca(k1,l1,k2,l2)=(sumcca(k1,l1,k2,l2)-covca(k1,l1)*covca(k2,l2))/n;

covccb(k1,l1,k2,l2)=(sumccb(k1,l1,k2,l2)-covcb(k1,l1)*covcb(k2,l2))/n;

covcc(k1,l1,k2,l2)=(sumcc(k1,l1,k2,l2)-covc(k1,l1)*covc(k2,l2))/n;

covaaaa(k1,l1,k2,l2)=(sumaaaa(k1,l1,k2,l2)-cova(k1,l1)*cova(k2,l2))/n;

covbbbb(k1,l1,k2,l2)=(sumbbbb(k1,l1,k2,l2)-covb(k1,l1)*covb(k2,l2))/n;

covaacc(k1,l1,k2,l2)=(sumxx(k1,l1,k2,l2)-cova(k1,l1)*covc(k2,l2))/n;

end;

end;

end;

end;

delta1=varb(1,2)/covb(1,2)**2;

delta2=covcc(2,3,2,3)/covc(2,3)**2;

delta3=varb(2,3)/covb(2,3)**2;

delta4=covcc(1,1,1,1)/covc(1,1)**2;

delta5=-2*covbb(2,1,3)/(covb(1,2)*covb(2,3));

delta6=-2*covcc(1,1,2,3)/(covc(1,1)*covc(2,3));

delta=varb(1,2)/covb(1,2)**2+covcc(2,3,2,3)/covc(2,3)**2+varb(2,3)/

covb(2,3)**2+covcc(1,1,1,1)/covc(1,1)**2-2*covbb(2,1,3)/(covb(1,2)

*covb(2,3))-2*covcc(1,1,2,3)/(covc(1,1)*covc(2,3));

sigmasq_s=cova(2,3)**2*covb(1,2)/(cova(1,3)*covb(2,3));

sigmasq_d=covb(1,2)*covb(2,3)/covb(1,3);

sigmasq_t=sigmasq_s+sigmasq_d/4;

/* A2 */

lambda_tq=covb(1,2)*covc(2,3)/(covb(2,3)*covc(1,1));

lnlambda_tq=log(lambda_tq);

var_lnlambdatq=varb(1,2)/covb(1,2)**2+covcc(2,3,2,3)/covc(2,3)**2

+varb(2,3)/covb(2,3)**2+covcc(1,1,1,1)/covc(1,1)**2

-2*covbb(2,1,3)/(covb(1,2)*covb(2,3))-2*covcc(1,1,2,3)

/(covc(1,1)*covc(2,3));

se_lnlambdatq=sqrt(var_lnlambdatq);

se_lambdatq=se_lnlambdatq*lambda_tq;

c1=lnlambda_tq-1.96*se_lnlambdatq;

c2=lnlambda_tq+1.96*se_lnlambdatq;

ci1=exp(c1);

ci2=exp(c2);

/* A4 */

var_sigmasq_s=sigmasq_s**2*(4*vara(2,3)/cova(2,3)**2+varb(1,2)/covb(1,2)**2

+vara(1,3)/cova(1,3)**2+varb(2,3)/covb(2,3)**2

-4*covaa(3,1,2)/(cova(1,3)*cova(2,3))-2*covbb(2,1,3)

/(covb(1,2)*covb(2,3)));

/* A5 */

var_sigmasq_d=sigmasq_d**2*(varb(1,2)/covb(1,2)**2+varb(2,3)/covb(2,3)**2

+varb(1,3)/covb(1,3)**2+2*covbb(2,1,3)/(covb(1,2)*covb(2,3))

-2*covbb(1,2,3)/(covb(1,2)*covb(1,3))-2*covbb(3,1,2)

/(covb(1,3)*covb(2,3)));

/* A6 */

var_sigmasq_t=var_sigmasq_s+var_sigmasq_d/4;

do k=1 to 3;

do l=1 to 3;

/* A7 */

cova_sigmat(k,l)=sigmasq_s*(2*covaaaa(k,l,2,3)/cova(2,3)-

covaaaa(k,l,1,3)/cova(1,3));

/* A7.2 */

covb_sigmas(k,l)=sigmasq_s*(covbbbb(k,l,1,2)/covb(1,2)-

covbbbb(k,l,2,3)/covb(2,3));

/* A7.3 */

covb_sigmad(k,l)=sigmasq_d*(covbbbb(k,l,1,2)/covb(1,2)+

covbbbb(k,l,2,3)/covb(2,3)-covbbbb(k,l,1,3)/covb(1,3));

/* A7.1 */

covb_sigmat(k,l)=covb_sigmas(k,l)+covb_sigmad(k,l)/4;

/* A7.4 */

covc_sigmat(k,l)=sigmasq_s*2*(covaacc(k,l,2,3)/cova(2,3)-

covaacc(k,l,1,3)/cova(1,3));

end;

end;

/* A8 */

beta_q=covb(1,2)*covc(2,3)/(covb(2,3)*sigmasq_t);

var_beta_q=beta_q**2*(varb(1,2)/covb(1,2)**2+covcc(2,3,2,3)/covc(2,3)**2

+varb(2,3)/covb(2,3)**2+var_sigmasq_t/sigmasq_t**2

-2*covbb(2,1,3)/(covb(1,2)*covb(2,3))-2*covb_sigmat(1,2)

/(covb(1,2)*sigmasq_t)-2*covc_sigmat(2,3)

/(covc(2,3)*sigmasq_t)+2*covb_sigmat(2,3)

/(covb(2,3)*sigmasq_t));

/* A9 */

beta_m=covc(2,3)/sigmasq_t;

var_beta_m=beta_m**2*(covcc(2,3,2,3)/covc(2,3)**2+var_sigmasq_t/sigmasq_t**2

-2*covc_sigmat(2,3)/(covc(2,3)*sigmasq_t));

/* A10.1 */

var_betaq_sigmad=(beta_q**2*sigmasq_d)**2*(9*varb(1,2)/covb(1,2)**2

+4*covcc(2,3,2,3)/

covc(2,3)**2+varb(2,3)/covb(2,3)**2+4*var_sigmasq_t/sigmasq_t**2

+varb(1,3)/covb(1,3)**2-6*covbb(2,1,3)/(covb(1,2)*covb(2,3))

-12*covb_sigmat(1,2)/(covb(1,2)*sigmasq_t)-6*covbb(1,2,3)

/(covb(1,2)*covb(1,3))-8*covc_sigmat(2,3)/(covc(2,3)*sigmasq_t)

+4*covb_sigmat(2,3)/(covb(2,3)*sigmasq_t)

+2*covbb(3,1,2)/(covb(2,3)*covb(1,3))+4*covb_sigmat(1,3)

/(covb(1,3)*sigmasq_t));

/* A10.2 */

covb_betaq_sigmat=(beta_q**2*sigmasq_d)*(3*covbb(1,1,2)/covb(1,2)

-covbbbb(1,1,2,3)/covb(2,3)

-2*covb_sigmat(1,1)/sigmasq_t-covbb(1,1,3)/covb(1,3));

/* A10 */

sigmasq_eq=(covb(1,1)-beta_q**2*sigmasq_d)/2;

var_sigmasq_eq=(varb(1,1)+var_betaq_sigmad-2*covb_betaq_sigmat)/4;

/* A11 */

sigmasq_er=(covb(2,2)-sigmasq_d)/2;

var_sigmasq_er=(varb(2,2)+var_sigmasq_d-2*covb_sigmad(2,2))/4;

/* A12.1 */

var_betasqm_sigmasqd=beta_m**4*sigmasq_d**2*(2*covcc(2,3,2,3)/covc(2,3)**2

+4*var_sigmasq_t/sigmasq_t**2+varb(1,2)/covb(1,2)**2

+varb(2,3)/covb(2,3)**2+varb(1,3)/covb(1,3)**2

-8*covc_sigmat(2,3)/(covc(2,3)*sigmasq_t)-4*covb_sigmat(1,2)

/(covb(1,2)*sigmasq_t)-4*covb_sigmat(2,3)/(covb(2,3)*sigmasq_t)

+4*covb_sigmat(1,3)/(covb(1,3)*sigmasq_t)+2*covbb(2,1,3)

/(covb(1,2)*covb(2,3))-2*covbb(1,2,3)/(covb(1,2)*covb(1,3))

-2*covbb(3,1,2)/(covb(1,3)*covb(2,3)));

/* A12.2 */

covb_betasqm_sigmasqd=beta_m**2*sigmasq_d*(covbbbb(1,2,3,3)/covb(1,2)

+covbb(3,2,3)/covb(2,3)-covbb(3,1,3)/covb(1,3)

-2*covb_sigmat(3,3)/sigmasq_t);

/* A12 */

sigmasq_em=(covb(3,3)-beta_m**2*sigmasq_d)/2;

var_sigmasq_em=(varb(3,3)+var_betasqm_sigmasqd-2*covb_betasqm_sigmasqd)/4;

/* A13.1 */

var_betasqq_sigmasqs=(beta_q**2*sigmasq_s)**2*(9*varb(1,2)/covb(1,2)**2

+4*covcc(2,3,2,3)/covc(2,3)**2+4*vara(2,3)/cova(2,3)**2

+9*varb(2,3)/covb(2,3)**2+4*var_sigmasq_t/sigmasq_t**2

+vara(1,3)/cova(1,3)**2-18*covbb(2,1,3)/(covb(1,2)*covb(2,3))

-12*covb_sigmat(1,2)/(covb(1,2)*sigmasq_t)+8*covaacc(2,3,2,3)

/(cova(2,3)*covc(2,3))-8*covc_sigmat(2,3)/(covc(2,3)*sigmasq_t)

-4*covaacc(1,3,2,3)/(cova(1,3)*covc(2,3))-8*cova_sigmat(2,3)/(cova(2,3)

*sigmasq_t)-4*covaa(3,1,2)/(cova(1,3)*cova(2,3))

+12*covb_sigmat(2,3)/(covb(2,3)*sigmasq_t)+4*cova_sigmat(1,3)

/(cova(1,3)*sigmasq_t));

/* A13.2 */

cova11_betasqq_sigmasqs=beta_q**2*sigmasq_s*(2*covaacc(1,1,2,3)

/covc(2,3)+2*covaaaa(1,1,2,3)/cova(2,3)-2*cova_sigmat(1,1)

/sigmasq_t -covaa(1,1,3)/cova(1,3));

/* A13.3 */

cova11_sigmasqeq=-1.0*beta_q**2*sigmasq_d*(covaacc(1,1,2,3)/covc(2,3)

-cova_sigmat(1,1)/sigmasq_t);

/* A13.5 */

cov_betasqq_sigmasqs_b11=beta_q**2*sigmasq_s*(3*covbb(1,1,2)/covb(1,2)

-3*covbbbb(1,1,2,3)/covb(2,3)-2*covb_sigmat(1,1)/sigmasq_t);

/* A13.7 */

cov_lnbetaq_lnsigmasqs=varb(1,2)/covb(1,2)**2-2*covbb(2,1,3)/

(covb(1,2)*covb(2,3))+2*covac(2,2,3)/(covc(2,3)*cova(2,3))

-covaacc(1,3,2,3)/(cova(1,3)*covc(2,3))+varb(2,3)/covb(2,3)**2

-2*cova_sigmat(2,3)/(cova(2,3)*sigmasq_t)-covb_sigmat(1,2)

/(covb(1,2)*sigmasq_t)+cova_sigmat(1,3)/(cova(1,3)*sigmasq_t)

+covb_sigmat(2,3)/(covb(2,3)*sigmasq_t);

/* A13.8 */

cov_lnbetaq_lnsigmasqd= varb(1,2)/covb(1,2)**2-covbb(1,2,3)/(covb(1,2)*covb(1,3))

-varb(2,3)/covb(2,3)**2+covbb(3,1,2)/(covb(1,3)*covb(2,3))

-covb_sigmat(1,2)/(covb(1,2)*sigmasq_t)-covb_sigmat(2,3)

/(covb(2,3)*sigmasq_t)+covb_sigmat(1,3)/(covb(1,3)*sigmasq_t);

/* A13.9 */

cov_lnsigmasqs_lnsigmasqd=varb(1,2)/covb(1,2)**2-covbb(1,2,3)/(covb(1,2)*covb(1,3))

-varb(2,3)/covb(2,3)**2+covbb(3,1,2)/(covb(1,3)*covb(2,3));

/* A13.6 */

cov_betaq_sigmas_betaq_sigmad=beta_q**4*sigmasq_s*sigmasq_d

*(4*var_beta_q/beta_q**2+2*cov_lnbetaq_lnsigmasqs

+2*cov_lnbetaq_lnsigmasqd+cov_lnsigmasqs_lnsigmasqd);

/* A13.4 */

cov_betasqq_sigmasqs_sigmasqeq=(cov_betasqq_sigmasqs_b11

-cov_betaq_sigmas_betaq_sigmad)/2;

/* A13 */

var_sigmasq_q=vara(1,1)+var_betasqq_sigmasqs+var_sigmasq_eq/4

-2*cova11_betasqq_sigmasqs-cova11_sigmasqeq

+cov_betasqq_sigmasqs_sigmasqeq;

/* A14.1 */

cova22_sigmasqs=sigmasq_s*(2*covaa(2,2,3)/cova(2,3)-covaaaa(1,3,2,2)

/cova(1,3));

/* A14.3 */

covb22_sigmasqs=sigmasq_s*(covbb(2,1,2)/covb(1,2)-covbb(2,2,3)/covb(2,3));

/* A14.4 */

cov_sigmasqs_sigmasqd=sigmasq_s*sigmasq_d*cov_lnsigmasqs_lnsigmasqd;

/* A14.2 */

cov_sigmasqs_sigmasqer=(covb22_sigmasqs-cov_sigmasqs_sigmasqd)/2;

/* A14 */

var_sigmasq_r=vara(2,2)+var_sigmasq_s+var_sigmasq_er/4

-2*cova22_sigmasqs+cov_sigmasqs_sigmasqer;

/* A15.1 */

var_betasqm_sigmasqs=beta_m**4*sigmasq_s**2*(4*covcc(2,3,2,3)

/covc(2,3)**2+4*vara(2,3)/cova(2,3)**2+varb(1,2)/covb(1,2)**2

+4*var_sigmasq_t/sigmasq_t**2+vara(1,3)/cova(1,3)**2

+varb(2,3)/covb(2,3)**2+8*covaacc(2,3,2,3)/(cova(2,3)*covc(2,3))

-8*covc_sigmat(2,3)/(covc(2,3)*sigmasq_t)-4*covaacc(1,3,2,3)

/(cova(1,3)*covc(2,3))-8*cova_sigmat(2,3)/(cova(2,3)*sigmasq_t)

-4*covaa(3,1,2)/(cova(1,3)*cova(2,3)) -4*covb_sigmat(1,2)

/(covb(1,2)*sigmasq_t)-2*covbb(2,1,3)/(covb(1,2)*covb(2,3))

+4*cova_sigmat(1,3)/(cova(1,3)*sigmasq_t)+4*covb_sigmat(2,3)

/(covb(2,3)*sigmasq_t));

/* A15.2 */

cova33_betasqm_sigmasqs=beta_m**2*sigmasq_s*(2*covac(3,2,3)/covc(2,3)

+2*covaa(3,2,3)/cova(2,3)-2*cova_sigmat(3,3)/sigmasq_t

-covaa(3,1,3)/cova(1,3));

/* A15.3 */

cova33_sigmasqem=-1*beta_m**2*sigmasq_d*(covac(3,2,3)/covc(2,3)

-cova_sigmat(3,3)/sigmasq_t);

/* A15.5 */

cov_betasqm_sigmasqs_b33=beta_m**2*sigmasq_s*(covbbbb(1,2,3,3)/covb(1,2)

-2*covb_sigmat(3,3)/sigmasq_t-covbb(3,2,3)/covb(2,3));

/* A15.8 */

cov_betam_sigmasqs=beta_m*sigmasq_s*(2*covaacc(2,3,2,3)/(covc(2,3)

*cova(2,3))-covaacc(2,3,1,3)/(covc(2,3)*cova(1,3))

-2*cova_sigmat(2,3)/(cova(2,3)*sigmasq_t)-covb_sigmat(1,2)

/(covb(1,2)*sigmasq_t)+cova_sigmat(1,3)/(cova(1,3)*sigmasq_t)

+covb_sigmat(2,3)/(covb(2,3)*sigmasq_t));

/* A15.9 */

cov_sigmasqs_sigmasqd=sigmasq_s*sigmasq_d*cov_lnsigmasqs_lnsigmasqd;

/* A15.10 */

cov_betam_sigmasqd=beta_m*sigmasq_d*(-1*covb_sigmat(1,2)

/(covb(1,2)*sigmasq_t)-covb_sigmat(2,3)/(covb(2,3)*sigmasq_t)

+covb_sigmat(1,3)/(covb(1,3)*sigmasq_t));

/* A15.7 */

cov_betam_sigmas_betam_sigmad=beta_m**4*sigmasq_s*sigmasq_d

*(4*var_beta_m/beta_m**2+4*cov_betam_sigmasqs

/(beta_m*sigmasq_s)+4*cov_betam_sigmasqd/(beta_m*sigmasq_d)

+cov_sigmasqs_sigmasqd/(sigmasq_s*sigmasq_d));

/* A15.4 */

cov_betasqm_sigmasqs_sigmasqem=(cov_betasqm_sigmasqs_b33

-cov_betam_sigmas_betam_sigmad)/2;

/* A15 */

var_sigmasqm=vara(3,3)+var_betasqm_sigmasqs+var_sigmasq_em/4

-2*cova33_betasqm_sigmasqs-cova33_sigmasqem

+cov_betasqm_sigmasqs_sigmasqem;

/* A16.3 */

sigmasq_q=cova(1,1)-beta_q**2*sigmasq_s-sigmasq_eq/2;

sigmasq_r=cova(2,2)-sigmasq_s-sigmasq_er/2;

sigmasq_m=cova(3,3)-beta_m**2*sigmasq_s-sigmasq_em/2;

sigma_qr=cova(1,2)-beta_q*sigmasq_s;

rhoqr=sigma_qr/(sigmasq_q*sigmasq_r)**0.5;

cova12_betaq_sigmasqs=beta_q*sigmasq_s*(covaacc(1,2,2,3)

/covc(2,3)+2*covaa(2,1,3)/cova(2,3)-covaa(1,2,3)/cova(1,3)

-cova_sigmat(1,2)/sigmasq_t);

/* A16.2 */

var_betaq_sigmasqs=beta_q**2*sigmasq_s**2*(4*varb(1,2)/covb(1,2)**2

+covcc(2,3,2,3)/covc(2,3)**2+4*vara(2,3)/cova(2,3)**2

+4*varb(2,3)/covb(2,3)**2+vara(1,3)/cova(1,3)**2

+var_sigmasq_t/sigmasq_t**2-8*covbb(2,1,3)/(covb(1,2)*covb(2,3))

-4*covb_sigmat(1,2)/(covb(1,2)*sigmasq_t)+4*covaacc(2,3,2,3)

/(covc(2,3)*cova(2,3))-2*covaacc(2,3,1,3)/(covc(2,3)*cova(1,3))

-2*covc_sigmat(2,3)/(covc(2,3)*sigmasq_t)-4*covaa(3,2,1)

/(cova(2,3)*cova(1,3))-4*cova_sigmat(2,3)

/(cova(2,3)*sigmasq_t)+4*covb_sigmat(2,3)

/(covb(2,3)*sigmasq_t)+2*cova_sigmat(1,3)

/(cova(1,3)*sigmasq_t));

/* A16.1 */

var_lnsigma_qr=(vara(1,2)+var_betaq_sigmasqs-2*cova12_betaq_sigmasqs)

/sigma_qr**2;

/* A16.4 */

var_lnsigma_q=var_sigmasq_q/(4*sigmasq_q**2);

/* A16.5 */

var_lnsigma_r=var_sigmasq_r/(4*sigmasq_r**2);

/* A16.8 */

cova12_betasqq_sigmasqs=beta_q**2*sigmasq_s*(2*covaacc(1,2,2,3)

/covc(2,3)+2*covaa(2,1,3)/cova(2,3)-covaa(1,2,3)/cova(1,3)

-2*cova_sigmat(1,2)/sigmasq_t);

/* A16.9 */

cova12_sigmasq_eq=-1*beta_q**2*sigmasq_d*(2*covaacc(1,2,2,3)/covc(2,3)

-2*cova_sigmat(1,2)/sigmasq_t)/2;

/* A16.10 */

cova11_betaq_sigmasqs=beta_q*sigmasq_s*(covaacc(1,1,2,3)/covc(2,3)

+2*covaaaa(1,1,2,3)/cova(2,3)-covaa(1,1,3)/cova(1,3)

-cova_sigmat(1,1)/sigmasq_t);

/* A17.1 */

cov_betaq_sigmasqs=beta_q*sigmasq_s*(varb(1,2)/covb(1,2)**2

-covbb(2,1,3)/(covb(1,2)*covb(2,3))+2*covaacc(2,3,2,3)

/(covc(2,3)*cova(2,3))-covaacc(2,3,1,3)/(covc(2,3)*cova(1,3))

-covbb(2,1,3)/(covb(2,3)*covb(1,2))+varb(2,3)/covb(2,3)**2

-2*cova_sigmat(2,3)/(cova(2,3)*sigmasq_t)-covb_sigmat(1,2)

/(covb(1,2)*sigmasq_t)+cova_sigmat(1,3)/(cova(1,3)*sigmasq_t)

+covb_sigmat(2,3)/(covb(2,3)*sigmasq_t));

/* A17 */

cov_betaqsigmasqs_betaqsigmasqs=beta_q**3*sigmasq_s**2

*(2*var_beta_q/beta_q**2+var_sigmasq_s/sigmasq_s**2

+3*cov_betaq_sigmasqs/(beta_q*sigmasq_s));

/* A18.1 */

cov_betaqsigmasqs_b11=beta_q*sigmasq_s*(2*covbb(1,1,2)/covb(1,2)

-2*covbbbb(1,1,2,3)/covb(2,3)-covb_sigmat(1,1)/sigmasq_t);

/* A18.3 */

cov_betaq_sigmasqd=beta_q*sigmasq_d*(varb(1,2)/covb(1,2)**2

-covbb(1,2,3)/(covb(1,2)*covb(1,3))-varb(2,3)/covb(2,3)**2

+covbb(3,2,1)/(covb(2,3)*covb(1,3))-covb_sigmat(1,2)

/(covb(1,2)*sigmasq_t)-covb_sigmat(2,3)/(covb(2,3)*sigmasq_t)

+covb_sigmat(1,3)/(covb(1,3)*sigmasq_t));

/* A18.2 */

cov_betaqsigmasqs_betaqsigmasqd=beta_q**3*sigmasq_s*sigmasq_d

*(2*var_beta_q/beta_q**2+2*cov_betaq_sigmasqs

/(beta_q*sigmasq_s)+cov_betaq_sigmasqd/(beta_q*sigmasq_d)

+cov_sigmasqs_sigmasqd/(sigmasq_s*sigmasq_d));

/* A18 */

cov_betaqsigmasqs_sigmasqeq=(cov_betaqsigmasqs_b11

-cov_betaqsigmasqs_betaqsigmasqd)/2;

/* A16.7 */

cov_sigmaqr_sigmasqq=covaa(1,1,2)-cova12_betasqq_sigmasqs

-cova12_sigmasq_eq/2-cova11_betaq_sigmasqs

+cov_betaqsigmasqs_betaqsigmasqs+cov_betaqsigmasqs_sigmasqeq/2;

/* A16.6 */

cov_lnsigmaqr_lnsigmaq=cov_sigmaqr_sigmasqq/(2*sigma_qr*sigmasq_q);

/* A19.3 */

cova22_betaqsigmasqs=beta_q*sigmasq_s*(covaacc(2,2,2,3)/covc(2,3)

+2*covaa(2,2,3)/cova(2,3)-covaaaa(2,2,1,3)/cova(1,3)

-cova_sigmat(2,2)/sigmasq_t);

/* A19.4 */

cov_betaqsigmasqs_sigmasqs=beta_q*sigmasq_s**2*(cov_betaq_sigmasqs

/(beta_q*sigmasq_s)+var_sigmasq_s/sigmasq_s**2);

/* A19.6 */

cov_betaqsigmasqs_b22=beta_q*sigmasq_s*(2*covbb(2,1,2)/covb(1,2)

-2*covbb(2,3,2)/covb(2,3)-covb_sigmat(2,2)/sigmasq_t);

/* A19.7 */

cov_betaqsigmasqs_sigmasqd=beta_q*sigmasq_s*sigmasq_d

*(cov_lnbetaq_lnsigmasqd+cov_lnsigmasqs_lnsigmasqd);

/* A19.5 */

cov_betaqsigmasqs_sigmasqer=(cov_betaqsigmasqs_b22

-cov_betaqsigmasqs_sigmasqd)/2;

/* A19.1 */

cov_sigmaqr_sigmasqr=covaa(2,1,2)-cova_sigmat(1,2)-cova22_betaqsigmasqs

+cov_betaqsigmasqs_sigmasqs+cov_betaqsigmasqs_sigmasqer/2;

/* A19 */

cov_lnsigmaqr_lnsigmar=cov_sigmaqr_sigmasqr/(2*sigma_qr*sigmasq_r);

/* A20.2 */

cova22_betasqqsigmasqs=beta_q**2*sigmasq_s*(2*covaacc(2,2,2,3)

/covc(2,3)+2*covaa(2,2,3)/cova(2,3)-covaaaa(2,2,1,3)/cova(1,3)

-2*cova_sigmat(2,2)/sigmasq_t);

cov_betasqqsigmasqs_sigmasqs=beta_q**2*sigmasq_s**2

*(2*cov_betaq_sigmasqs/(beta_q*sigmasq_s)+var_sigmasq_s

/sigmasq_s**2);

/* A20.4 */

cov_betasqqsigmasqs_b22=beta_q**2*sigmasq_s*(3*covbb(2,2,1)/covb(1,2)

-3*covbb(2,2,3)/covb(2,3)-2*covb_sigmat(2,2)/sigmasq_t);

/* A20.5 */

cov_betasqqsigmasqs_sigmasqd=beta_q**2*sigmasq_s*sigmasq_d

*(2*cov_lnbetaq_lnsigmasqd+cov_lnsigmasqs_lnsigmasqd);

/* A20.3 */

cov_betasqqsigmasqs_sigmasqer=(cov_betasqqsigmasqs_b22

-cov_betasqqsigmasqs_sigmasqd)/2;

/* A20.6 */

cova22_sigmasqeq=beta_q**2*sigmasq_d*(cova_sigmat(2,2)/sigmasq_t

-covaacc(2,2,2,3)/covc(2,3));

/* A20.8 */

cov_betasqqsigmasqd_sigmasqs=beta_q**2*sigmasq_s*sigmasq_d*

(2*cov_lnbetaq_lnsigmasqs+cov_lnsigmasqs_lnsigmasqd);

/* A20.7 */

cov_sigmasqs_sigmasqeq=(covb_sigmas(1,1)-cov_betasqqsigmasqd_sigmasqs)

/2;

/* A20.10 */

covb22_betasqqsigmasqd=beta_q**2*sigmasq_d*(3*covbb(2,2,1)/covb(1,2)

-covbb(2,2,3)/covb(2,3)-covbbbb(2,2,1,3)/covb(1,3)

-2*covb_sigmat(2,2)/sigmasq_t);

/* A20.11 */

cov_sigmasqd_betasqqsigmasqd=beta_q**2*sigmasq_d*(var_sigmasq_d

/sigmasq_d+2*cov_betaq_sigmasqd/beta_q);

/* A20.9 */

cov_sigmasqeq_sigmasqer=(covbbbb(1,1,2,2)-covb_sigmad(1,1)

-covb22_betasqqsigmasqd+cov_sigmasqd_betasqqsigmasqd)/4;

/* A20.1 */

cov_sigmasqq_sigmasqr=covaaaa(1,1,2,2)-cova_sigmat(1,1)

-cova22_betasqqsigmasqs+cov_betasqqsigmasqs_sigmasqs

+cov_betasqqsigmasqs_sigmasqer/2

-cova22_sigmasqeq/2+cov_sigmasqs_sigmasqeq/2

+cov_sigmasqeq_sigmasqer/4;

/* A20 */

cov_lnsigmaq_lnsigmar=cov_sigmasqq_sigmasqr/(4*sigmasq_q*sigmasq_r);

/* A16 */

var_rhoqr=rhoqr**2*(var_lnsigma_qr+var_lnsigma_q+var_lnsigma_r

-2*cov_lnsigmaqr_lnsigmaq-2*cov_lnsigmaqr_lnsigmar

+2*cov_lnsigmaq_lnsigmar);

/* A21.11 */

cov_sigmasqs_sigmasqt=var_sigmasq_s+cov_sigmasqs_sigmasqd/4;

/* A21.12 */

cov_sigmasqer_sigmasqt=(covb_sigmat(2,2)-cov_sigmasqs_sigmasqd

-var_sigmasq_d/4)/2;

/* A21.10 */

cov_sigmasqr_sigmasqt=cova_sigmat(2,2)-cov_sigmasqs_sigmasqt

-cov_sigmasqer_sigmasqt/2;

/* A21.9 */

cov_ab22_sigmasqr=cova_sigmat(2,2)-covb_sigmat(2,2)/4

-cov_sigmasqr_sigmasqt;

/* A21.8 */

cov_lnab22_lnsigmasqt=cov_ab22_sigmasqr/(sigmasq_t*(cova(2,2)

-covb(2,2)/4-sigmasq_r));

/* A21.7 */

var_lnsigmasqt=var_sigmasq_t/sigmasq_t**2;

/* A21.6 */

cov_b22_sigmasqer=(varb(2,2)-covb_sigmad(2,2))/2;

/* A21.5 */

cov_b22_sigmasqr=-1*covb_sigmas(2,2)-cov_b22_sigmasqer/2;

/* A21.3 */

cov_a22_sigmasqr=vara(2,2)-cova_sigmat(2,2);

/* A21.2 */

var_ab22_sigmasqr=vara(2,2)+varb(2,2)/16+var_sigmasq_r

-2*cov_a22_sigmasqr+cov_b22_sigmasqr/2;

/* A21.1 */

var_lnab22_sigmasqr=var_ab22_sigmasqr

/((cova(2,2)-covb(2,2)/4-sigmasq_r)**2);

/* A21 */

rho_t=(cova(2,2)-covb(2,2)/4-sigmasq_r)/sigmasq_t;

var_rho_t=rho_t**2*(var_lnab22_sigmasqr+var_lnsigmasqt

-2*cov_lnab22_lnsigmasqt);

/* A23.9 */

do m=1 to 3;

do j=1 to 2;

cov_sigmasqs_xbar(m,j)=sigmasq_s*(2*cov_a_xbar(2,3,m,j)/cova(2,3)

+cov_b_xbar(1,2,m,j)/covb(1,2)-cov_a_xbar(1,3,m,j)/cova(1,3)

-cov_b_xbar(2,3,m,j)/covb(2,3));

/* A23.10 */

cov_sigmasqd_xbar(m,j)=sigmasq_d*(cov_b_xbar(1,2,m,j)/covb(1,2)

+cov_b_xbar(2,3,m,j)/covb(2,3)-cov_b_xbar(1,3,m,j)/covb(1,3));

/* A23.8 */

cov_sigmasqt_xbar(m,j)=cov_sigmasqs_xbar(m,j)+cov_sigmasqd_xbar(m,j)

/4;

/* A23.7 */

cov_betaq_xbar(m,j)=beta_q*(cov_b_xbar(1,2,m,j)/covb(1,2)

+cov_c_xbar(2,3,m,j)/covc(2,3)-cov_b_xbar(2,3,m,j)/covb(2,3)

-cov_sigmasqt_xbar(m,j)/sigmasq_t);

/* A23.2 */

if j=1 then

var_betaq_xbar(j)=beta_q**2*varxbar(2,1)+xbarj1(2)**2*var_beta_q;

if j=2 then

var_betaq_xbar(j)=beta_q**2*varxbar(2,2)+xbarj2(2)**2*var_beta_q;

/* A23.11 */

if j=1 then

cov_qbarj_betaqrbarj(j)=beta_q*cov_xbark_xbarl(1,2,j);

if j=2 then

cov_qbarj_betaqrbarj(j)=beta_q*cov_xbark_xbarl(1,2,j);

/* A23 */

var_alphaq(j)=varxbar(1,j)+var_betaq_xbar(j)

-2 * cov_qbarj_betaqrbarj(j);

end;

end;

do j=1 to 2;

/* A24.2 */

if j=1 then

var_betam_xbar(j)=beta_m**2*varxbar(2,1)+xbarj1(2)**2*var_beta_m;

if j=2 then

var_betam_xbar(j)=beta_m**2*varxbar(2,2)+xbarj2(2)**2*var_beta_m;

/* A24.3 */

cov_mbarj_betamrbarj(j)=beta_m*cov_xbark_xbarl(2,3,j);

/* A24 */

var_alpham(j)=varxbar(3,j)+var_betam_xbar(j)-2*cov_mbarj_betamrbarj(j);

end;

mu_t1=xbarj1(2);

mu_t2=xbarj2(2);

alpha_q1=xbarj1(1)-beta_q*mu_t1;

alpha_q2=xbarj2(1)-beta_q*mu_t2;

alpha_m1=xbarj1(3)-beta_m*mu_t1;

alpha_m2=xbarj2(3)-beta_m*mu_t2;

* lambda_rq=(beta_q*sigmasq_t+sigma_qr)/covc(1,1);

lambda_rq=covc(1,2)/covc(1,1);

var_lambda_rq=(covc(2,2)-covc(1,2)**2/covc(1,1))/((n-2)*covc(1,1));

se_lambda_rq=sqrt(var_lambda_rq);

one=1;

xcovc12=covc(1,2);

xcovc12c=covc12(1,2);

xcovc23=covc(2,3);

xcovc23c=covc12(2,3);

output partb;

end;

data outp;

set partb;

array vr varxbar3 varxbar4 var_alphaq1 var_alphaq2 var_alpham1 var_alpham2

var_sigmasq_s var_sigmasq_d var_sigmasq_t var_beta_q var_beta_m

var_sigmasq_eq var_sigmasq_er var_sigmasq_em var_sigmasq_q

var_sigmasq_r var_sigmasqm var_rho_t var_rhoqr;

array se sexbar3 sexbar4 se_alphaq1 se_alphaq2 se_alpham1 se_alpham2

se_sigmasq_s se_sigmasq_d se_sigmasq_t se_beta_q se_beta_m

se_sigmasq_eq se_sigmasq_er se_sigmasq_em se_sigmasq_q

se_sigmasq_r se_sigmasqm se_rho_t se_rhoqr;

do over vr;

se=sqrt(vr);

end;

array ee alpha_q1 alpha_q2 alpha_m1 alpha_m2 mu_t1 mu_t2

lambda_rq ;

array aa se_alphaq1 se_alphaq2 se_alpham1 se_alpham2

sexbar3 sexbar4 se_lambda_rq;

array ca ci1_alphaq1 ci1_alphaq2 ci1_alpham1 ci1_alpham2

ci1xbar3 ci1xbar4 ci1_lambda_rq;

array cb ci2_alphaq1 ci2_alphaq2 ci2_alpham1 ci2_alpham2

ci2xbar3 ci2xbar4 ci2_lambda_rq;

do over aa;

ca=ee-1.96*aa;

cb=ee+1.96*aa;

end;

put lambda_rq= se_lambda_rq= ci1_lambda_rq= ci2_lambda_rq= ;

array bb se_sigmasq_s se_sigmasq_d se_sigmasq_t se_beta_q se_beta_m

se_sigmasq_eq se_sigmasq_er se_sigmasq_em se_sigmasq_q

se_sigmasq_r se_sigmasqm se_rhoqr se_rho_t;

array ba ci1_sigmasq_s ci1_sigmasq_d ci1_sigmasq_t ci1_beta_q ci1_beta_m

ci1_sigmasq_eq ci1_sigmasq_er ci1_sigmasq_em ci1_sigmasq_q

ci1_sigmasq_r ci1_sigmasqm ci1_rhoqr ci1_rho_t;

array bc ci2_sigmasq_s ci2_sigmasq_d ci2_sigmasq_t ci2_beta_q ci2_beta_m

ci2_sigmasq_eq ci2_sigmasq_er ci2_sigmasq_em ci2_sigmasq_q

ci2_sigmasq_r ci2_sigmasqm ci2_rhoqr ci2_rho_t;

array dd sigmasq_s sigmasq_d sigmasq_t beta_q beta_m sigmasq_eq

sigmasq_er sigmasq_em sigmasq_q sigmasq_r sigmasq_m

rhoqr rho_t;

do over bb;

ba=log(dd)-1.96*(bb/dd);

bc=log(dd)+1.96*(bb/dd);

ba=exp(ba);

bc=exp(bc);

end;

beta_m_z=beta_m/se_beta_m;

beta_q_z=beta_q/se_beta_q;

pval_beta_q=2*(1-probnorm(abs(beta_q_z)));

pval_beta_m=2*(1-probnorm(abs(beta_m_z)));

file print;

put @10 'Parameter estimates from correlated measurement error raw model'/;

put @10 'Parameter'/ @10 'Type' @25 'Parameter' @47 'Estimate' @60 'SE'

@67 '95% CI' @81 'p-value'/;

put @10 'Intercept' @25 'mu t1' mu_t1 45-52 .4 sexbar3 55-62 .4

@65 '(' ci1xbar3 66-71 .4 ',' ci2xbar3 73-78 .4 ')' /;

put @25 'mu t2' mu_t2 45-52 .4 sexbar4 55-62 .4

@65 '(' ci1xbar4 66-71 .4 ',' ci2xbar4 73-78 .4 ')' /;

put @25 'alpha q1' alpha_q1 45-52 .4 se_alphaq1 55-62 .4

@65 '(' ci1_alphaq1 66-71 .4 ',' ci2_alphaq1 73-78 .4 ')' /;

put @25 'alpha q2' alpha_q2 45-52 .4 se_alphaq2 55-62 .4

@65 '(' ci1_alphaq2 66-71 .4 ',' ci2_alphaq2 73-78 .4 ')' /;

put @25 'alpha m1' alpha_m1 45-52 .4 se_alpham1 55-62 .4

@65 '(' ci1_alpham1 66-71 .4 ',' ci2_alpham1 73-78 .4 ')' /;

put @25 'alpha m2' alpha_m2 45-52 .4 se_alpham2 55-62 .4

@65 '(' ci1_alpham2 66-71 .4 ',' ci2_alpham2 73-78 .4 ')' /;

put @10 'Regression' @25 'beta q ' beta_q 45-52 .4 se_beta_q 55-62 .4

@65 '(' ci1_beta_q 66-71 .4 ',' ci2_beta_q 73-78 .4 ')'

pval_beta_q 82-86 .3 /;

put @25 'beta m ' beta_m 45-52 .4 se_beta_m 55-62 .4

@65 '(' ci1_beta_m 66-71 .4 ',' ci2_beta_m 73-78 .4 ')'

pval_beta_m 82-86 .3 /;

put @10 'variance' @25 'sigma s**2' sigmasq_s 45-52 .4 se_sigmasq_s 55-62 .4

@65 '(' ci1_sigmasq_s 66-71 .4 ',' ci2_sigmasq_s 73-78 .4 ')' /

@10 'component' ;

put @25 'sigma d**2' sigmasq_d 45-52 .4 se_sigmasq_d 55-62 .4

@65 '(' ci1_sigmasq_d 66-71 .4 ',' ci2_sigmasq_d 73-78 .4 ')' /;

put @25 'sigma t**2' sigmasq_t 45-52 .4 se_sigmasq_t 55-62 .4

@65 '(' ci1_sigmasq_t 66-71 .4 ',' ci2_sigmasq_t 73-78 .4 ')' /;

put @25 'sigma q**2' sigmasq_q 45-52 .4 se_sigmasq_q 55-62 .4

@65 '(' ci1_sigmasq_q 66-71 .4 ',' ci2_sigmasq_q 73-78 .4 ')' /;

put @25 'sigma r**2' sigmasq_r 45-52 .4 se_sigmasq_r 55-62 .4

@65 '(' ci1_sigmasq_r 66-71 .4 ',' ci2_sigmasq_r 73-78 .4 ')' /;

put @25 'sigma m**2' sigmasq_m 45-52 .4 se_sigmasqm 55-62 .4

@65 '(' ci1_sigmasqm 66-71 .4 ',' ci2_sigmasqm 73-78 .4 ')' /;

put @25 'sigma eq**2' sigmasq_eq 45-52 .4 se_sigmasq_eq 55-62 .4

@65 '(' ci1_beta_m 66-71 .4 ',' ci2_beta_m 73-78 .4 ')' /;

put @25 'sigma er**2' sigmasq_er 45-52 .4 se_sigmasq_er 55-62 .4

@65 '(' ci1_sigmasq_er 66-71 .4 ',' ci2_sigmasq_er 73-78 .4 ')' /;

put @25 'sigma em**2' sigmasq_em 45-52 .4 se_sigmasq_em 55-62 .4

@65 '(' ci1_sigmasq_em 66-71 .4 ',' ci2_sigmasq_em 73-78 .4 ')' /;

put @10 'correlation' @25 'rho qr' rhoqr 45-52 .4 se_rhoqr 55-62 .4

@65 '(' ci1_rhoqr 66-71 .4 ',' ci2_rhoqr 73-78 .4 ')' /;

put @25 'rho t' rho_t 45-52 .4 se_rho_t 55-62 .4

@65 '(' ci1_rho_t 66-71 .4 ',' ci2_rho_t 73-78 .4 ')' /;

put @10 'deattenuation'@25 'lambda rq' lambda_rq 45-52 .3 se_lambda_rq 55-62 .4

@65 '(' ci1_lambda_rq 66-71 .4 ',' ci2_lambda_rq 73-78 .4 ')' /

@10 'factor';

put @25 'lambda tq' lambda_tq 45-52 .3 se_lambdatq 55-62 .4

@65 '(' ci1 66-71 .4 ',' ci2 73-78 .4 ')' /;

%mend rawTable;