/*********************************************************** * Program Name: rawAdjstTable.sas * * Program Date: Mar., 2008 * * * * Note : * * This program assumes there is a variable * * called "visit" in the input data * * and visit has value 1 for 1st visit; * * 2: for the 2nd visit * * please make sure the input data has * * variable "visit" with value "1" and "2" * * * * If the input data has missing values in * * one of the visits, the pair of data will * * deleted for the analysis * * * * Programmed by: Marion McPhee and Rong Chen * ***********************************************************/ %macro rawAdjstTable(dsn=, /* input data */ q=, /* surrogate measure of intake */ r=, /* reference measure of intake */ m=, /* biomarker */ covar= /* covariates for adjusting the model */); %numvar(&covar); %let nvar2=%eval(&nvar*2+1); %let nvar3=%eval(&nvar+2); /* for the future use */ /* check for missing */ /* and delete the missing pairs */ data _null_; %let k=1; %let nmissf=%str(); %do %while(&k<=&nvar); %let cvar=%scan(&covar,&k); %if &k=1 %then %do; %let nmissf=&cvar; %end; %else %do; %let nmissf=&nmissf,&cvar; %end; %let k=%eval(&k+1); %end; %let nmissf=&nmissf,&q,&r,&m; run; %put &nmissf; data main; set &dsn end=eof ; if nmiss(&nmissf)^=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; %let count=1; %let cvar = %scan(&covar, &count); %do %while(&cvar ne); %if &count=1 %then %do; proc transpose data=main out=temp(drop=_NAME_) prefix=&cvar; by id; var &cvar; run; data temp; set temp; nn1=&cvar.1; nn2=&cvar.2; call symput("mn1", "&cvar.1"); call symput("mn2", "&cvar.2"); call symput("lbl1", "&cvar"); run; %end; %else %do; proc transpose data=main out=tempn(drop=_NAME_) prefix=&cvar; by id; var &cvar; run; data tempn; set tempn; nn%eval(&count*2-1)=&cvar.1; nn%eval(&count*2)=&cvar.2; call symput("mn%eval(&count*2-1)", "&cvar.1"); call symput("mn%eval(&count*2)", "&cvar.2" ); call symput("lbl&count", "&cvar"); run; data temp; merge temp tempn; by id; run; %end; %let count=%eval(&count+1); %let cvar=%scan(&covar, &count); %end; 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 temp; by id; run; proc means data=temp noprint; var nn1 - nn%eval(&nvar2-1); output out=zz mean=mn1-mn&nvar2; run; proc datasets nolist; delete temp1 temp2 temp3 temp4 tempn; run; /* count total # of observations for one visit */ data _null_; set temp end=eof; if eof then do; call symput('totn', trim(left(_n_))); end; run; proc mixed method=ml noclprint data=main; ods output covB=ricovb SolutionF=riparm r=rmtrx; class id; model ri=&covar visit /solution covb corrb; repeated /type=cs sub=id r rcorr; proc mixed method=ml noclprint data=main; ods output covB=qicovb SolutionF=qiparm r=qmtrx; class id; model qi=&covar visit /solution covb corrb; repeated /type=cs sub=id r rcorr; proc mixed method=ml noclprint data=main; ods output covB=micovb SolutionF=miparm r=mmtrx; class id; model mi=&covar visit /solution covb corrb; repeated /type=cs sub=id r rcorr; data parm1; set riparm end=eof; retain beta_ri1-beta_ri&nvar3 num 0; num=num+1; array beta_ri(*) beta_ri1-beta_ri&nvar3; beta_ri(num)=estimate; one=1; if eof then output; data parm2; set qiparm end=eof; retain beta_qi1-beta_qi&nvar3 num 0; num=num+1; array beta_qi(*) beta_qi1-beta_qi&nvar3; beta_qi(num)=estimate; one=1; if eof then output; data parm3; set miparm end=eof; retain beta_mi1-beta_mi&nvar3 num 0; num=num+1; array beta_mi(*) beta_mi1-beta_mi&nvar3; beta_mi(num)=estimate; one=1; if eof then output; data _null_; %let resid_ri1=&r.1; %let resid_ri2=&r.2; %let resid_qi1=&q.1; %let resid_qi2=&q.2; %let resid_mi1=&m.1; %let resid_mi2=&m.2; %let count=1; %let k=1; %do %while(&count<%eval(&nvar+1)); %let count=%eval(&count+1); %let resid_ri1=%cmpres(&resid_ri1)-beta_ri&count*(%cmpres(&&mn&k)-mn&k); %let resid_qi1=%cmpres(&resid_qi1)-beta_qi&count*(%cmpres(&&mn&k)-mn&k); %let resid_mi1=%cmpres(&resid_mi1)-beta_mi&count*(%cmpres(&&mn&k)-mn&k); %let k=%eval(&count*2-1); %let j=%eval((&count-1)*2); %let resid_ri2=%cmpres(&resid_ri2)-beta_ri&count*(%cmpres(&&mn&j)-mn&j); %let resid_qi2=%cmpres(&resid_qi2)-beta_qi&count*(%cmpres(&&mn&j)-mn&j); %let resid_mi2=%cmpres(&resid_mi2)-beta_mi&count*(%cmpres(&&mn&j)-mn&j); %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 mn1-mn%eval(&nvar*2) beta_ri1-beta_ri&nvar3 beta_mi1-beta_mi&nvar3 beta_qi1-beta_qi&nvar3; if _n_=1 then set zz; if _n_=1 then set parm1; if _n_=1 then set parm2; if _n_=1 then set parm3; if _n_=1 then put beta_ri1= beta_ri2= beta_mi1= beta_mi2= ; retain sumxi1 sumxi2 sumxi3 sumxdi1-sumxdi3 sumx11-sumx13 sumx21-sumx23 nn 0; visit=0; resid_ri1=&resid_ri1-beta_ri&nvar3*visit; resid_qi1=&resid_qi1-beta_qi&nvar3*visit; resid_mi1=&resid_mi1-beta_mi&nvar3*visit; visit=1; resid_ri2=&resid_ri2-beta_ri&nvar3*visit; resid_qi2=&resid_qi2-beta_qi&nvar3*visit; resid_mi2=&resid_mi2-beta_mi&nvar3*visit; array a(3) resid_qi1 resid_ri1 resid_mi1; array b(3) resid_qi2 resid_ri2 resid_mi2; array d(3) resid_qi resid_ri resid_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; proc means data=raw1 noprint; id one; var &q.1 &q.2 &r.1 &r.2 &m.1 &m.2; output out=mn mean=mqi1 mqi2 mri1 mri2 mmi1 mmi2 stderr=seqi1 seqi2 seri1 seri2 semi1 semi2; 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) resid_qi1 resid_ri1 resid_mi1; array b(3) resid_qi2 resid_ri2 resid_mi2; 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) resid_qi resid_ri resid_mi; 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; output partb; end; data beta; set partb; keep beta_m beta_q; data mix2; set main; retain beta_m beta_q; if _n_=1 then set beta; if _n_=1 then put beta_m= beta_q= ; qr_ij=&q-beta_q*(&r); mr_ij=&m-beta_m*(&r); proc mixed method=ml noclprint data=mix2; ods output covB=qrcovb SolutionF=qrparm; class id; model qr_ij=&covar visit/solution covb corrb; repeated /type=cs sub=id r rcorr; proc mixed method=ml noclprint data=mix2; ods output covB=mrcovb SolutionF=mrparm; class id; model mr_ij=&covar visit/solution covb corrb; repeated /type=cs sub=id r rcorr; data parm4; set qrparm end=eof; retain beta_qr1-beta_qr&nvar3 se_qr1-se_qr&nvar3 pval_qr1-pval_qr&nvar3 num 0; num=num+1; array beta_qr(*) beta_qr1-beta_qr&nvar3; array se_qr(*) se_qr1-se_qr&nvar3; array pval_qr(*) pval_qr1-pval_qr&nvar3; beta_qr(num)=estimate; se_qr(num)=stderr; pval_qr(num)=probt; one=1; if eof then output; proc print data=parm4; run; data parm5; set mrparm end=eof; retain beta_mr1-beta_mr&nvar3 se_mr1-se_mr&nvar3 pval_mr1-pval_mr&nvar3 num 0; num=num+1; array beta_mr(*) beta_mr1-beta_mr&nvar3; array se_mr(*) se_mr1-se_mr&nvar3; array pval_mr(*) pval_mr1-pval_mr&nvar3; beta_mr(num)=estimate; se_mr(num)=stderr; pval_mr(num)=probt; one=1; if eof then output; proc print data=parm5; run; data covar4; set qrcovb end=eof; retain qrvar1 qrvar&nvar3 qrcovar%eval(&nvar3*2+1) one; if effect='Intercept' then qrvar1=col1; if effect='Intercept' then qrcovar%eval(&nvar3*2+1)=col&nvar3; if effect='visit' then qrvar&nvar3=col&nvar3; one=1; if eof then output; data covar5; set mrcovb end=eof; retain mrvar1 mrvar&nvar3 mrcovar%eval(&nvar3*2+1) one; if effect='Intercept' then mrvar1=col1; if effect='Intercept' then mrcovar%eval(&nvar3*2+1)=col&nvar3; if effect='visit' then mrvar&nvar3=col&nvar3; one=1; if eof then output; title; %let kk=1; %let jj=1; data outp; *set partb; merge partb parm4 parm5 covar4 covar5 mn; by one; 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))); array beta_qr(*) beta_qr1-beta_qr&nvar3; array se_qr(*) se_qr1-se_qr&nvar3; array beta_mr(*) beta_mr1-beta_mr&nvar3; array se_mr(*) se_mr1-se_mr&nvar3; array pval_qr(*) pval_qr1-pval_qr&nvar3; array pval_mr(*) pval_mr1-pval_mr&nvar3; array ci1_qr(*) ci1_qr1-ci1_qr&nvar3; array ci2_qr(*) ci2_qr1-ci2_qr&nvar3; array ci1_mr(*) ci1_mr1-ci1_mr&nvar3; array ci2_mr(*) ci2_mr1-ci2_mr&nvar3; alpha_q1=beta_qr(1); se_alpha_q1=se_qr(1); alpha_q2=beta_qr(1)+beta_qr(&nvar3); var_q2=qrvar1+qrvar&nvar3+2*qrcovar%eval(&nvar3*2+1); se_alpha_q2=sqrt(var_q2); alpha_m1=beta_mr(1); se_alpha_m1=se_mr(1); alpha_m2=beta_mr(1)+beta_mr(&nvar3); var_m2=mrvar1+mrvar&nvar3+2*mrcovar%eval(&nvar3*2+1); se_alpha_m2=sqrt(var_m2); do ii=1 to dim(se_qr); ci1_qr(ii)=beta_qr(ii)-1.96*se_qr(ii); ci2_qr(ii)=beta_qr(ii)+1.96*se_qr(ii); ci1_mr(ii)=beta_mr(ii)-1.96*se_mr(ii); ci2_mr(ii)=beta_mr(ii)+1.96*se_mr(ii); end; file print; put @10 'Parameter estimates from correlated error measurement raw model'/; put @35 'Independent' ; put @10 'Parameter'/ @10 'Type' @35 'Variable' @25 'Parameter' @57 'Estimate' @70 'SE' @77 '95% CI' @91 'p-value'/; put @10 'Intercept' @25 'mu t1' mri1 55-62 .4 seri1 65-72 .4 /; * @75 '(' ci1xbar3 76-81 .4 ',' ci2xbar3 83-88 .4 ')' /; put @25 'mu t2' mri2 55-62 .4 seri2 65-72 .4 /; * @75 '(' ci1xbar4 76-81 .4 ',' ci2xbar4 83-88 .4 ')' /; put @25 'alpha q1' alpha_q1 55-62 .4 se_alpha_q1 65-72 .4 /; * @75 '(' ci1_alphaq1 76-81 .4 ',' ci2_alphaq1 83-88 .4 ')' /; put @25 'alpha q2' alpha_q2 55-62 .4 se_alpha_q2 65-72 .4 /; * @75 '(' ci1_alphaq2 76-81 .4 ',' ci2_alphaq2 83-88 .4 ')' /; put @25 'alpha m1' alpha_m1 55-62 .4 se_alpha_m1 65-72 .4 /; * @75 '(' ci1_alpham1 76-81 .4 ',' ci2_alpham1 83-88 .4 ')' /; put @25 'alpha m2' alpha_m2 55-62 .4 se_alpha_m2 65-72 .4 /; * @75 '(' ci1_alpham2 76-81 .4 ',' ci2_alpham2 83-88 .4 ')' /; put @10 'Regression' @25 'beta q ' beta_q 55-62 .4 se_beta_q 65-72 .4 @75 '(' ci1_beta_q 76-81 .4 ',' ci2_beta_q 83-88 .4 ')' pval_beta_q 92-96 .3 /; put @25 'beta m ' beta_m 55-62 .4 se_beta_m 65-72 .4 @75 '(' ci1_beta_m 76-81 .4 ',' ci2_beta_m 83-88 .4 ')' pval_beta_m 92-96 .3 /; put @25 'gamma q' ; %do %while(&kk <= &nvar); put @35 "&&lbl&kk" beta_qr%eval(&kk+1) 55-62 .4 se_qr%eval(&kk+1) 65-72 .4 @75 '(' ci1_qr%eval(&kk+1) 76-81 .3 ',' ci2_qr%eval(&kk+1) 83-88 .3 ')' pval_qr(%eval(&kk+1)) 92-96 .3 /; %let kk=%eval(&kk+1); %end; put @35 ' ' /; put @25 'gamma m'; %do %while(&jj <= &nvar); put @35 "&&lbl&jj" beta_mr%eval(&jj+1) 55-62 .4 se_mr%eval(&jj+1) 65-72 .4 @75 '(' ci1_mr%eval(&jj+1) 76-81 .3 ',' ci2_mr%eval(&jj+1) 83-88 .3 ')' pval_mr(%eval(&jj+1)) 92-96 .3 /; %let jj=%eval(&jj+1); %end; put @35 ' ' /; put @10 'variance' @25 'sigma s**2' sigmasq_s 55-62 .4 se_sigmasq_s 65-72 .4 @75 '(' ci1_sigmasq_s 76-81 .4 ',' ci2_sigmasq_s 83-88 .4 ')' / @10 'component' ; put @25 'sigma d**2' sigmasq_d 55-62 .4 se_sigmasq_d 65-72 .4 @75 '(' ci1_sigmasq_d 76-81 .4 ',' ci2_sigmasq_d 83-88 .4 ')' /; put @25 'sigma t**2' sigmasq_t 55-62 .4 se_sigmasq_t 65-72 .4 @75 '(' ci1_sigmasq_t 76-81 .4 ',' ci2_sigmasq_t 83-88 .4 ')' /; put @25 'sigma q**2' sigmasq_q 55-62 .4 se_sigmasq_q 65-72 .4 @75 '(' ci1_sigmasq_q 76-81 .4 ',' ci2_sigmasq_q 83-88 .4 ')' /; put @25 'sigma r**2' sigmasq_r 55-62 .4 se_sigmasq_r 65-72 .4 @75 '(' ci1_sigmasq_r 76-81 .4 ',' ci2_sigmasq_r 83-88 .4 ')' /; put @25 'sigma m**2' sigmasq_m 55-62 .4 se_sigmasqm 65-72 .4 @75 '(' ci1_sigmasqm 76-81 .4 ',' ci2_sigmasqm 83-88 .4 ')' /; put @25 'sigma eq**2' sigmasq_eq 55-62 .4 se_sigmasq_eq 65-72 .4 @75 '(' ci1_beta_m 76-81 .4 ',' ci2_beta_m 83-88 .4 ')' /; put @25 'sigma er**2' sigmasq_er 55-62 .4 se_sigmasq_er 65-72 .4 @75 '(' ci1_sigmasq_er 76-81 .4 ',' ci2_sigmasq_er 83-88 .4 ')' /; put @25 'sigma em**2' sigmasq_em 55-62 .4 se_sigmasq_em 65-72 .4 @75 '(' ci1_sigmasq_em 76-81 .4 ',' ci2_sigmasq_em 83-88 .4 ')' /; put @10 'correlation' @25 'rho qr' rhoqr 55-62 .4 se_rhoqr 65-72 .4 @75 '(' ci1_rhoqr 76-81 .4 ',' ci2_rhoqr 83-88 .4 ')' /; put @25 'rho t' rho_t 55-62 .4 se_rho_t 65-72 .4 @75 '(' ci1_rho_t 76-81 .4 ',' ci2_rho_t 83-88 .4 ')' /; put @10 'deattenuation'@25 'lambda rq' lambda_rq 55-62 .3 se_lambda_rq 65-72 .4 @75 '(' ci1_lambda_rq 76-81 .4 ',' ci2_lambda_rq 83-88 .4 ')' / @10 'factor'; put @25 'lambda tq' lambda_tq 55-62 .3 se_lambdatq 65-72 .4 @75 '(' ci1 76-81 .4 ',' ci2 83-88 .4 ')' /; run; %mend; %macro numvar(list); %global nvar; %let k=1; %let nvar=1; %let cvar = %scan(&list, &k); %do %while(&cvar ne); %let nvar=%eval(&nvar+1); %let k=%eval(&k+1); %let cvar=%scan(&list, &k); %end; %let nvar=%eval(&nvar-1); %mend numvar;