!------------------------------------------------------------- ! For Data Analysis ! Association between Y and X controlled Z !------------------------------------------------------ ! Y(n,np) ==> The first group of variables; ! X(n,nq) ==> The second group of variables; ! Z(n,nr) ==> The third group of variables. !------------------------------------------------------------- program main parameter(TINY=1e-7) character*4 ref !Last ordered observation. !ref='last' means last unordered observation character*300 DataSetName, outputName !hold the name of the data file and result integer np,nq,nr,niter,miss,missind real eps real dat(5000) !Temperoarily store a record of the data. integer unita integer nfield,ncaty,maxycat double precision result(5) ! five test result in p-values integer,Allocatable:: Yfield(:),Xfield(:),Zfield(:),yCfield(:) integer,Allocatable:: index(:),catX(:),catY(:) real,Allocatable:: y(:,:),x(:,:),z(:,:),stemp(:),yc(:,:),yr(:,:) real,Allocatable:: sdy(:),sdx(:) double precision,Allocatable:: gam(:),vgam(:,:),gam1(:),vgam1(:,:),gam2(:),vgam2(:,:) double precision,Allocatable:: pf(:),gf(:),pb(:),gb(:) !pf,gf---forward likelihood estimator !pb,gb---backward likelihood estimator double precision pvaluenorm real,allocatable:: score(:),sqvar(:),ce1(:,:),ce2(:,:),cov1(:,:,:),cov2(:,:,:),xc(:,:) real,allocatable:: beta(:),var(:,:) open(unit=5,file='ModelOR',status='old') !modeling information read(5,*)DataSetName !file name that stores the data read(5,*)OutputName ! file name for storing the result read(5,*)nfield,np,nq,nr,niter,eps,maxycat !write(*,*)nfield,np,nq,nr,niter,eps Allocate(Yfield(np),YCfield(np),Xfield(nq),Zfield(nr)) read(5,*)(Yfield(i),i=1,np) !locations of Y in the data file read(5,*)(YCfield(i),i=1,np)!locations of Y categorized in the data file read(5,*)(Xfield(i),i=1,nq) !locations of X in the data file read(5,*)(Zfield(i),i=1,nr) !locations of Z in the data file read(5,*) miss ! label of missing values unita=4 open(unit=unita,file=OutputName,status='unknown') write(unita,*)'This analysis computes the association of variables ' write(unita,*)(Yfield(k),k=1,np) write(unita,*)'with variables ' write(unita,*)(Xfield(k),k=1,nq) write(unita,*)'conditional on variables' write(unita,*)(zfield(k),k=1,nr) write(unita,*)'Label for missing values' write(unita,*) miss write(*,*)'This analysis computes the association of variables ' write(*,*)(Yfield(k),k=1,np) write(*,*)'with variables ' write(*,*)(Xfield(k),k=1,nq) write(*,*)'conditional on variables' write(*,*)(zfield(k),k=1,nr) write(*,*)'Label for missing values' write(*,*) miss close(unit=5) !determine the number of records in the data open(unit=3,file=DataSetName,status='old') nn=0 1020 read(3,*,end=1010)(dat(j),j=1,nfield) nn=nn+1 goto 1020 1010 continue close(unit=3) allocate(index(nn),catX(nn+1),catY(nn+1)) allocate(y(nn,np),yc(nn,np),x(nn,nq),z(nn,nr),stemp(nn)) allocate(sdy(np),sdx(nq)) allocate(gam(np*nq),vgam(np*nq,np*nq)) allocate(gam1(np*nr),vgam1(np*nr,np*nr)) allocate(gam2(nq*nr),vgam2(nq*nr,nq*nr)) allocate(pf(nn),gf(np*nr),pb(nn),gb(nq*nr)) allocate(ce1(nn,np),cov1(nn,np,np),ce2(nn,nq),cov2(nn,nq,nq),xc(nn,np+nr)) allocate(sqvar(nq),score(nq)) allocate(beta(nq+nr+1),var(nq+nr+1,nq+nr+1)) ref='lord' ! last on in ordered form n=nn write(unita,*)"sample size= ", n write(unita,*)" SNP "," missing% ", " Logistic ", " Linear ", " SPORT-r ", " SPORT-p ", " SPORT-d" ! count the number of categories in Y order(y,nn,np,index) rtemp=y(index(1),np) ! collect the numbers of categories ncat=1 do i=2,n if(abs(y(index(i),np)-rtemp).GT.TINY) then rtemp=y(index(i),np) ncat=ncat+1 endif enddo write(unita,*) "The number of categories in Y are:", ncat yr=y if(ncat>maxycat) then call roundv(yr,yc,n,p,maxycat) endif do nsnp=1,nfield-Xfield(1)+1 !number of snps fields write(*,*)"SNP index:", Nsnp !write(unita,*)"SNP index:", Nsnp n=0 open(unit=3,file=DataSetName,status='old') do i=1,nn read(3,*)(dat(j),j=1,nfield) missind=0 do j=1,nq ! No interaction if(abs(dat(Xfield(j)+(nsnp-1)*nq)-miss)<10e-6) then !Assume SNPs are listed consecutively missind=1 exit endif enddo if(missind==0) then n=n+1 do j=1,nq ! No interaction x(n,j)=dat(Xfield(j)+(nsnp-1)*nq) !Assume SNPs are listed consecutively enddo do j=1,np y(n,j)=dat(Yfield(j)) YC(n,j)=dat(YCfield(j)) enddo do j=1,nr z(n,j)=dat(Zfield(j)) enddo !do j=1,nq-1 ! with interaction ! x(n,j)=dat(Xfield(j)+(nsnp-1)*(nq-1)) !Assume SNPs are listed consecutively !enddo !x(n,nq)=x(n,1)*z(n,2) !interaction term else endif enddo close(unit=3) do j=1,nq ! rescale X to get standardized odds ratio parameter sum1=0.0 sum2=0.0 do i=1,n sum1=sum1+x(i,j) sum2=sum2+x(i,j)*x(i,j) enddo sum1=sum1/n sdx(j)=sqrt(sum2/n-sum1*sum1) !! sdx(j)=1.0 !prevent from rescaling !! write(*,*)'std for X(j) ',j,sdx(j) do i=1,n x(i,j)=x(i,j)/sdx(j) enddo enddo !! do 200 i=1,n !! round=50 ! rounded to 1/round !! do 170 j=1,nq !! x(i,j)=int(x(i,j)*round)/round !c170 continue !! do 180 j=1,np !! y(i,j)=int(y(i,j)*round)/round !c180 continue !c200 continue do j=1,np ! rescale Y to get standardized odds ratio parameter sum1=0.0 sum2=0.0 do i=1,nn sum1=sum1+y(i,j) sum2=sum2+y(i,j)*y(i,j) enddo sum1=sum1/nn sdy(j)=sqrt(sum2/nn-sum1*sum1) !! sdy(j)=1.0 ! prevent from rescaling !! write(*,*)'std for Y(j) ',j,sdy(j) do i=1,nn y(i,j)=y(i,j)/sdy(j) enddo enddo write(*,*)"Total N=",nn, ", Observed n=", n, ", Missing n=", nn-n !write(unita,*)"Total N=",nn, ", Observed n=", n, ", Missing n=", nn-n write(*,*)"case-control analysis" !write(unita,*)"case-control analysis" Call OR(yc,x,z,n,np,nq,nr,index,gam,vgam,gf,pf,niter,eps,caty,ncaty,ref,unita) !Call OR(x,yc,z,n,nq,np,nr,index,gam,vgam,gf,pf,niter,eps,caty,ncaty,ref,unita) ! write(*,*)"Parameter estimates" ! write(*,55)((gam((k-1)*nq+j),j=1,nq),k=1,np) ! write(*,*)'Standard errors:' ! write(*,55)((sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) write(*,*)'Wald tests' write(*,55)((gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) write(*,*)'Two-sided P-values' write(*,56)((2*pvaluenorm(gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j))),j=1,nq),k=1,np) result(1)=2*pvaluenorm(gam(1)/sqrt(vgam(1,1))) ! write(*,*)"Variance-Covariance estimates" ! do l=1,np ! do ll=1,nq ! lll=(l-1)*nq+ll ! write(*,55)((vgam((k-1)*nq+j,lll),j=1,nq),k=1,np) ! enddo ! enddo ! write(unita,*)"Odds ratio estimates" ! write(unita,55)((gam((k-1)*nq+j),j=1,nq),k=1,np) ! write(unita,*)'Standard errors:' ! write(unita,55)((sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) ! write(unita,*)'Wald tests' !write(unita,55)((gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) ! write(unita,*)'Two-sided P-values' !write(unita,56)((2*pvaluenorm(gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j))),j=1,nq),k=1,np) ! Linear model analysis write(*,*)"Linear regression analysis" !write(unita,*)"Linear regression analysis" xc(:,1:nq)=x xc(:,(nq+1):(nq+nr))=z Call Linear(y(:,1),xc,nn,nq+nr,beta,var,sigma1) write(*,55)(beta(k)/sqrt(var(k,k)),k=1,nq) write(*,56)(2*pvaluenorm(DBLE(beta(k)/sqrt(var(k,k)))),k=1,nq) result(2)=2*pvaluenorm(DBLE(beta(1)/sqrt(var(1,1)))) !write(unita,55)(beta(k)/sqrt(var(k,k)),k=1,nq) ! write(unita,*)'Two-sided P-values' !write(unita,56)(2*pvaluenorm(DBLE(beta(k)/sqrt(var(k,k)))),k=1,nq) write(*,*)"Retrospective dose-response analysis" !write(unita,*)"Retrospective dose-response analysis" !write(unita,*)"dose-response analysis" ! forward likelihood analysis Call OR(x,y,z,n,nq,np,nr,index,gam,vgam,gf,pf,niter,eps,caty,ncaty,ref,unita) !Call OR(x,y,z,n,nq,np,nr,index,gam,vgam,gf,pf,niter,eps,caty,ncaty,ref,unita) ! write(*,*)"Parameter estimates" ! write(*,55)((gam((k-1)*nq+j),j=1,nq),k=1,np) ! write(*,*)'Standard errors:' ! write(*,55)((sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) write(*,*)'Wald tests' write(*,55)((gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) write(*,*)'Two-sided P-values' write(*,56)((2*pvaluenorm(gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j))),j=1,nq),k=1,np) result(3)=2*pvaluenorm(gam(1)/sqrt(vgam(1,1))) ! write(*,*)"Variance-Covariance estimates" ! do l=1,np ! do ll=1,nq ! lll=(l-1)*nq+ll ! write(*,55)((vgam((k-1)*nq+j,lll),j=1,nq),k=1,np) ! enddo ! enddo ! write(unita,*)"Odds ratio estimates" ! write(unita,55)((gam((k-1)*nq+j),j=1,nq),k=1,np) ! write(unita,*)'Standard errors:' ! write(unita,55)((sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) ! write(unita,*)'Wald tests' ! write(unita,55)((gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) ! write(unita,*)'Two-sided P-values' ! write(unita,56)((2*pvaluenorm(gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j))),j=1,nq),k=1,np) ! write(unita,*)"Variance-Covariance estimates" ! do l=1,np ! do ll=1,nq ! lll=(l-1)*nq+ll ! write(unita,55)((vgam((k-1)*nq+j,lll),j=1,nq),k=1,np) ! enddo ! enddo ! write(*,*) "Complete the retrospective likelihood approach" ! write(unita,*) "Complete the retrospective likelihood approach" write(*,*)"Prospective dose-response analysis" !write(unita,*)"Prospective dose-response analysis" !write(unita,*)"dose-response analysis" ! forward likelihood analysis ! yr is y with possible rounding Call OR(yr,x,z,n,np,nq,nr,index,gam,vgam,gf,pf,niter,eps,caty,ncaty,ref,unita) !Call OR(x,y,z,n,nq,np,nr,index,gam,vgam,gf,pf,niter,eps,caty,ncaty,ref,unita) ! write(*,*)"Parameter estimates" ! write(*,55)((gam((k-1)*nq+j),j=1,nq),k=1,np) ! write(*,*)'Standard errors:' ! write(*,55)((sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) write(*,*)'Wald tests' write(*,55)((gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) write(*,*)'Two-sided P-values' write(*,56)((2*pvaluenorm(gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j))),j=1,nq),k=1,np) result(4)=2*pvaluenorm(gam(1)/sqrt(vgam(1,1))) ! write(*,*)"Variance-Covariance estimates" ! do l=1,np ! do ll=1,nq ! lll=(l-1)*nq+ll ! write(*,55)((vgam((k-1)*nq+j,lll),j=1,nq),k=1,np) ! enddo ! enddo ! write(unita,*)"Odds ratio estimates" ! write(unita,55)((gam((k-1)*nq+j),j=1,nq),k=1,np) ! write(unita,*)'Standard errors:' ! write(unita,55)((sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) ! write(unita,*)'Wald tests' !write(unita,55)((gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j)),j=1,nq),k=1,np) ! write(unita,*)'Two-sided P-values' ! write(unita,56)((2*pvaluenorm(gam((k-1)*nq+j)/sqrt(vgam((k-1)*nq+j,(k-1)*nq+j))),j=1,nq),k=1,np) ! write(unita,*)"Variance-Covariance estimates" ! do l=1,np ! do ll=1,nq ! lll=(l-1)*nq+ll ! write(unita,55)((vgam((k-1)*nq+j,lll),j=1,nq),k=1,np) ! enddo ! enddo ! write(*,*) "Complete the prospective likelihood approach" ! write(unita,*) "Complete the prospective likelihood approach" ! fit model and estimate the conditional means and variances for observed covariates. write(*,*)"Doubly robust test" !write(unita,*)"Doubly robust test" call ORsimple(yr,z,nn,np,nr,index,gam1,vgam1,pf,niter,eps,caty,ncaty,ref,unita,ce1,cov1) call ORsimple(x,z,nn,nq,nr,index,gam2,vgam2,pf,niter,eps,caty,ncaty,ref,unita,ce2,cov2) ! calculate the efficient test ! For np=1,nq=1 do k=1,nq score(k)=sum((Y(:,1)-ce1(:,1))*(X(:,k)-ce2(:,k))) sqvar(k)=sqrt(sum(cov1(:,1,1)*cov2(:,k,k))) enddo write(*,*)'Wald tests' write(*,55)(score(k)/sqvar(k),k=1,nq) write(*,*)'Two-sided P-values' write(*,56) (2*pvaluenorm(DBLE(score(k)/sqvar(k))),k=1,nq) result(5)=2*pvaluenorm(DBLE(score(1)/sqvar(1))) !write(unita,55)(score(k)/sqvar(k),k=1,nq) !write(unita,56)(2*pvaluenorm(DBLE(score(k)/sqvar(k))),k=1,nq) 55 format(1x,8(f10.4,2x)) 56 Format(10(ES10.1,2x)) write(unita,57)nsnp,(nn-n)*100.0/nn,(result(k),k=1,5) 57 Format(I6,2x,f6.1,2x,6(ES10.1,2x)) enddo close(unit=unita) end program main !--|------------------------------------------| ! | Direct Newton Raphson method for finding | ! | the maximizer in the linear OR model | !--|------------------------------------------| ! gam: odds ratio estimator ! vgam: variance estimate !--------------------------------------------------------- Subroutine OR(y,x,z,n,np,nq,nr,index,gam,vgam,gf,pf,niter,eps,cat,ncat,ref,unita) parameter(TINY=10e-6) integer n,np,nq,nr,niter,ncat,IRP ! integer reference point integer index(n),cat(n+1) character*4 ref integer unita real y(n,np),x(n,nq),z(n,nr),eps double precision gam(np*nq),vgam(np*nq,np*nq) double precision gf(np*nr),pf(n) ! gf---nuissance gamma; pf---nuissance probabilities integer itemp real rtemp double precision,allocatable:: para(:) !nonparametric parameter double precision,allocatable:: dint(:),cp(:,:)!conditional probability double precision,allocatable:: ce(:,:),cov(:,:,:),vtg(:,:) !conditional mean and variance double precision,allocatable:: der(:),der2(:,:),tder2(:,:),der22(:,:) !first and second derivatives double precision,allocatable:: tp(:),tder(:)!total parameter vector allocate(tder(np*nq)) allocate(vtg(np*nq,np*nq)) !!!!!Preparation call order(y,n,np,index) if(ref.eq.'last') then IRP=n elseif(ref.eq.'lord') then IRP=index(n) else write(*,*)'No reference point set' return endif rtemp=y(index(1),np) ! collect the numbers of categories itemp=1 cat(1)=1 do 60 i=2,n if(abs(y(index(i),np)-rtemp).GT.TINY) then rtemp=y(index(i),np) itemp=itemp+1 cat(itemp)=i ! The beginning index. endif 60 continue ncat=itemp cat(ncat+1)=n+1 !!!!!! Initialization na=np*nq nb=np*nr m=na+nb+ncat-1 allocate(para(ncat),dint(n),cp(n,ncat),ce(n,np),cov(n,np,np)) allocate(der(m),der2(m,m),tder2(m,m),der22(m,m),tp(m)) ! total parameter do 70 k=1,ncat-1 para(k)=log(1.0*(cat(k+1)-cat(k))/(cat(ncat+1)-cat(ncat))) !nonparametric parameter 70 continue para(ncat)=0 do 80 k=1,na gam(k)=0.0 !odds ratio parameter 80 continue do 85 k=1,nb gf(k)=0.0 !nuissance odds ratio parameter 85 continue DO 1000 IArep=1,niter ! overall iteration do 800 i=1,n ! initialize conditional probability, expectation, and covariance of Y. dint(i)=0.0 do 810 k=1,np ce(i,k)=0.0 do 820 l=1,np cov(i,k,l)=0.0 820 continue 810 continue do 830 j=1,ncat! ! compute denominator integral cp(i,j)=0 do 840 k=1,np do 845 l=1,nq cp(i,j)=cp(i,j)+gam((k-1)*nq+l)*(y(index(cat(j)),k)-y(IRP,k))*(x(i,l)-x(IRP,l)) 845 continue do 848 ll=1,nr cp(i,j)=cp(i,j)+gf((k-1)*nr+ll)*(y(index(cat(j)),k)-y(IRP,k))*(z(i,ll)-z(IRP,ll)) 848 continue 840 continue cp(i,j)=exp(cp(i,j)+para(j)) dint(i)=dint(i)+cp(i,j) !for denominator do 850 k=1,np ! for conditional probability, mean, and variance ce(i,k)=ce(i,k)+y(index(cat(j)),k)*cp(i,j) ! contrast y with the largest, butallow exponent to be subtracted by a constant do 860 l=1,np cov(i,k,l)=cov(i,k,l)+y(index(cat(j)),k)*y(index(cat(j)),l)*cp(i,j) ! contrast y with the largest, butallow exponent to be subtracted by a constant 860 continue 850 continue !write(*,*)IArep,j,ncat 830 continue do 870 j=1,ncat ! CONDITIONAL PROBABILITY cp(i,j)=cp(i,j)/dint(i) 870 continue do 880 k=1,np ! complete the conditional expectation calculation ce(i,k)=ce(i,k)/dint(i) do 890 l=1,k cov(i,k,l)=cov(i,k,l)/dint(i)-ce(i,k)*ce(i,l) cov(i,l,k)=cov(i,k,l) 890 continue 880 continue 800 continue ! ! Find the first and second derivatives for nonparametric part ! do 900 j=1,m if(j.lt.ncat) then der(j)=cat(j+1)-cat(j) else der(j)=0.0 endif do 910 k=1,j der2(j,k)=0.0 der2(k,j)=0.0 910 continue 900 continue do 90 i=1,n do 100 j=1,ncat-1 der(j)=der(j)-cp(i,j) do 110 l=1,j-1 der2(j,l)=der2(j,l)+cp(i,j)*cp(i,l) 110 continue der2(j,j)=der2(j,j)-cp(i,j)*(1-cp(i,j)) 100 continue ! ! Find the first and second deivatives for parametric part ! do 140 k=1,np do 150 j=1,nq kj=ncat-1+(k-1)*nq+j der(kj)=der(kj)+(y(i,k)-ce(i,k))*(x(i,j)-x(IRP,j)) do 160 kk=1,np do 170 jj=1,nq kkjj=ncat-1+(kk-1)*nq+jj der2(kj,kkjj)=der2(kj,kkjj)-cov(i,k,kk)*(x(i,j)-x(IRP,j))*(x(i,jj)-x(IRP,jj)) 170 continue 160 continue 150 continue do 180 j=1,nr kj=ncat-1+na+(k-1)*nr+j der(kj)=der(kj)+(y(i,k)-ce(i,k))*(z(i,j)-z(IRP,j)) do 190 kk=1,np do 200 jj=1,nq ! xz cross-product terms kkjj=ncat-1+(kk-1)*nq+jj der2(kj,kkjj)=der2(kj,kkjj)-cov(i,k,kk)*(z(i,j)-z(IRP,j))*(x(i,jj)-x(IRP,jj)) 200 continue 190 continue do 210 kk=1,np do 220 jj=1,nr ! zz terms kkjj=ncat-1+na+(kk-1)*nr+jj der2(kj,kkjj)=der2(kj,kkjj)-cov(i,k,kk)*(z(i,j)-z(IRP,j))*(z(i,jj)-z(IRP,jj)) 220 continue 210 continue 180 continue 140 continue ! ! Find the joint second derivative of nonparametric ! and parametric parts ! do 230 k=1,np do 240 j=1,nq kj=ncat-1+(k-1)*nq+j do 250 l=1,ncat-1 der2(kj,l)=der2(kj,l)-(y(index(cat(l)),k)-ce(i,k))*cp(i,l)*(x(i,j)-x(IRP,j)) 250 continue 240 continue do 260 j=1,nr kj=ncat-1+na+(k-1)*nr+j do 270 l=1,ncat-1 der2(kj,l)=der2(kj,l)-(y(index(cat(l)),k)-ce(i,k))*cp(i,l)*(z(i,j)-z(IRP,j)) 270 continue 260 continue 230 continue 90 continue ! sum over subjects ! reflecting the second derivatives do 120 j=1,m ! m=ncat-1+na+nb; na=np*nq; nb=np*nr do 130 l=1,j-1 der2(j,l)=-der2(j,l) der2(l,j)=der2(j,l) tder2(j,l)=der2(j,l) tder2(l,j)=der2(l,j) 130 continue der2(j,j)=-der2(j,j) tder2(j,j)=der2(j,j) !tder2 is reserved for computing influence scores 120 continue !!!! Find the inverse of the second derivative matrix call chol(der2,tp,m,m) !!!! parameter updating delta=0 do 280 j=1,m tp(j)=0 do 290 l=1,m tp(j)=tp(j)+der2(j,l)*der(l) 290 continue delta=delta+abs(tp(j)) 280 continue do 300 j=1,ncat-1 para(j)=para(j)+tp(j) 300 continue do 310 k=1,np do 320 j=1,nq kj=(k-1)*nq+j gam(kj)=gam(kj)+tp(ncat-1+kj) 320 continue 310 continue do 330 k=1,np do 340 j=1,nr kj=(k-1)*nr+j gf(kj)=gf(kj)+tp(ncat-1+na+kj) 340 continue 330 continue write(*,*)'Iteration=',IArep,' Convergence=',delta ! write(unita,*)'Iteration=',IArep,' Convergence=',delta ! write(*,*) (gam(j),j=1,na) ! write(*,*) (gf(j),j=1,nb) ! write(*,*) (para(j),j=1,ncat) ! pause !!!! check convergence if(delta.lt.eps) Then ! convergent ! (a). Find the variance for the gamma estimate do 350 k=1,np do 360 j=1,nq kj=(k-1)*nq+j do 370 kk=1,np do 380 jj=1,nq kkjj=(kk-1)*nq+jj vgam(kj,kkjj)=der2(kj+ncat-1,kkjj+ncat-1) 380 continue 370 continue 360 continue 350 continue ! do 336 k=1,na ! write(*,*)(vgam(k,j),j=1,na) !c336 continue ! (b). Find the attached probability for each observation do j=1,ncat pf(j)=para(j) enddo DEallocate(para,dint,cp,ce,cov,der,der2,tder2) DEallocate(tp,der22,vtg,tder) exit endif 1000 CONTINUE end subroutine OR !-------------------------------------------------- ! fit a linear regression model with fully data ! y---continuous variable !-------------------------------------------------- subroutine linear(y,x,n,p,beta,var,sigma) integer n,p real y(n),x(n,p) real beta(p+1),var(p+1,p+1) real sigma,temp double precision,allocatable::XX(:,:),XY(:),det(:),alpha(:) allocate(XX(p+1,p+1),XY(p+1),det(p+1),alpha(p+1)) det=0 XX=0 XY=0 do i=1,n do j=1,p XY(j)=XY(j)+x(i,j)*y(i) XX(j,j)=XX(j,j)+X(i,j)*X(i,j) do k=1,j-1 XX(j,k)=XX(j,k)+X(i,j)*X(i,k) enddo enddo XY(p+1)=XY(p+1)+y(i) XX(p+1,p+1)=XX(p+1,p+1)+1 do k=1,p XX(p+1,k)=XX(p+1,k)+X(i,k) enddo enddo do j=1,p+1 do k=1,j-1 XX(k,j)=XX(j,k) enddo enddo call chol(XX,det,p+1,p+1) alpha=matmul(XX,XY) sigma=0 do i=1,n temp=alpha(p+1) do j=1,p temp=temp+alpha(j)*x(i,j) enddo sigma=sigma+(y(i)-temp)*(y(i)-temp) enddo sigma=sigma/(n-p-1) ! write(*,*)'alpha=',alpha(1:p) ! write(*,*)'var=',(XX(j,j)*sigma,j=1,p) beta=alpha var=XX*sigma !! write(*,*)'sigma=',sigma ! beta=alpha(1:p)/sigma !odds ratio parameter ! do k=1,p ! var(k,k)=XX(k,k)/sigma+2*beta(k)*beta(k)/(n-p-1) ! do j=1,k-1 ! var(k,j)=XX(k,j)/sigma+2*beta(k)*beta(j)/(n-p-1) ! var(j,k)=var(k,j) ! enddo ! enddo !! write(*,*)'theta=',beta !! write(*,*)'std=',(sqrt(var(k,k)),k=1,p) DEallocate(XX,XY,det,alpha) end subroutine linear !--------------------------------------------------------------- ! Unconditional odds ratio model fit by Newton-Raphson !--------------------------------------------------------------- Subroutine ORsimple(y,x,n,np,nq,index,gam,vgam,pf,niter,eps,cat,ncat,ref,unita,ce,cov) parameter(TINY=10e-6) integer n,np,nq,niter,ncat,IRP ! integer reference point integer index(n),cat(n+1) character*4 ref integer unita real y(n,np),x(n,nq),eps double precision gam(np*nq),vgam(np*nq,np*nq) double precision pf(n) ! pf---nuissance probabilities real ce(n,np),cov(n,np,np) !conditional mean and variance of y given observed x integer itemp real rtemp double precision,allocatable:: para(:) !nonparametric parameter double precision,allocatable:: dint(:),cp(:,:)!conditional probability double precision,allocatable:: vtg(:,:) double precision,allocatable:: der(:),der2(:,:),tder2(:,:),der22(:,:) !first and second derivatives double precision,allocatable:: tp(:),tder(:)!total parameter vector allocate(tder(np*nq)) allocate(vtg(np*nq,np*nq)) !!!!!Preparation call order(y,n,np,index) if(ref.eq.'last') then IRP=n elseif(ref.eq.'lord') then IRP=index(n) else write(*,*)'No reference point set' return endif rtemp=y(index(1),np) ! collect the numbers of categories itemp=1 cat(1)=1 do 60 i=2,n if(abs(y(index(i),np)-rtemp).GT.TINY) then rtemp=y(index(i),np) itemp=itemp+1 cat(itemp)=i ! The beginning index. endif 60 continue ncat=itemp cat(ncat+1)=n+1 !!!!!! Initialization na=np*nq m=na+ncat-1 allocate(para(ncat),dint(n),cp(n,ncat)) allocate(der(m),der2(m,m),tder2(m,m),der22(m,m),tp(m)) ! total parameter do 70 k=1,ncat-1 para(k)=log(1.0*(cat(k+1)-cat(k))/(cat(ncat+1)-cat(ncat))) !nonparametric parameter 70 continue para(ncat)=0 do 80 k=1,na gam(k)=0.0 !odds ratio parameter 80 continue DO 1000 IArep=1,niter ! overall iteration do 800 i=1,n ! initialize conditional probability, expectation, and covariance of Y. dint(i)=0.0 do 810 k=1,np ce(i,k)=0.0 do 820 l=1,np cov(i,k,l)=0.0 820 continue 810 continue do 830 j=1,ncat! ! compute denominator integral cp(i,j)=0 do 840 k=1,np do 845 l=1,nq cp(i,j)=cp(i,j)+gam((k-1)*nq+l)*(y(index(cat(j)),k)-y(IRP,k))*(x(i,l)-x(IRP,l)) 845 continue 840 continue cp(i,j)=exp(cp(i,j)+para(j)) dint(i)=dint(i)+cp(i,j) !for denominator do 850 k=1,np ! for conditional probability, mean, and variance ce(i,k)=ce(i,k)+y(index(cat(j)),k)*cp(i,j) ! contrast y with the largest, butallow exponent to be subtracted by a constant do 860 l=1,np cov(i,k,l)=cov(i,k,l)+y(index(cat(j)),k)*y(index(cat(j)),l)*cp(i,j) ! contrast y with the largest, butallow exponent to be subtracted by a constant 860 continue 850 continue !write(*,*)IArep,j,ncat 830 continue do 870 j=1,ncat ! CONDITIONAL PROBABILITY cp(i,j)=cp(i,j)/dint(i) 870 continue do 880 k=1,np ! complete the conditional expectation calculation ce(i,k)=ce(i,k)/dint(i) do 890 l=1,k cov(i,k,l)=cov(i,k,l)/dint(i)-ce(i,k)*ce(i,l) cov(i,l,k)=cov(i,k,l) 890 continue 880 continue 800 continue ! ! Find the first and second derivatives for nonparametric part ! do 900 j=1,m if(j.lt.ncat) then der(j)=cat(j+1)-cat(j) else der(j)=0.0 endif do 910 k=1,j der2(j,k)=0.0 der2(k,j)=0.0 910 continue 900 continue do 90 i=1,n do 100 j=1,ncat-1 der(j)=der(j)-cp(i,j) do 110 l=1,j-1 der2(j,l)=der2(j,l)+cp(i,j)*cp(i,l) 110 continue der2(j,j)=der2(j,j)-cp(i,j)*(1-cp(i,j)) 100 continue ! ! Find the first and second deivatives for parametric part ! do 140 k=1,np do 150 j=1,nq kj=ncat-1+(k-1)*nq+j der(kj)=der(kj)+(y(i,k)-ce(i,k))*(x(i,j)-x(IRP,j)) do 160 kk=1,np do 170 jj=1,nq kkjj=ncat-1+(kk-1)*nq+jj der2(kj,kkjj)=der2(kj,kkjj)-cov(i,k,kk)*(x(i,j)-x(IRP,j))*(x(i,jj)-x(IRP,jj)) 170 continue 160 continue 150 continue 140 continue ! ! Find the joint second derivative of nonparametric ! and parametric parts ! do 230 k=1,np do 240 j=1,nq kj=ncat-1+(k-1)*nq+j do 250 l=1,ncat-1 der2(kj,l)=der2(kj,l)-(y(index(cat(l)),k)-ce(i,k))*cp(i,l)*(x(i,j)-x(IRP,j)) 250 continue 240 continue 230 continue 90 continue ! sum over subjects ! reflecting the second derivatives do 120 j=1,m ! m=ncat-1+na+nb; na=np*nq; nb=np*nr do 130 l=1,j-1 der2(j,l)=-der2(j,l) der2(l,j)=der2(j,l) tder2(j,l)=der2(j,l) tder2(l,j)=der2(l,j) 130 continue der2(j,j)=-der2(j,j) tder2(j,j)=der2(j,j) !tder2 is reserved for computing influence scores 120 continue !!!! Find the inverse of the second derivative matrix call chol(der2,tp,m,m) !!!! parameter updating delta=0 do 280 j=1,m tp(j)=0 do 290 l=1,m tp(j)=tp(j)+der2(j,l)*der(l) 290 continue delta=delta+abs(tp(j)) 280 continue do 300 j=1,ncat-1 para(j)=para(j)+tp(j) 300 continue do 310 k=1,np do 320 j=1,nq kj=(k-1)*nq+j gam(kj)=gam(kj)+tp(ncat-1+kj) 320 continue 310 continue !!!! check convergence if(delta.lt.eps) Then ! convergent ! (a). Find the variance for the gamma estimate do 350 k=1,np do 360 j=1,nq kj=(k-1)*nq+j do 370 kk=1,np do 380 jj=1,nq kkjj=(kk-1)*nq+jj vgam(kj,kkjj)=der2(kj+ncat-1,kkjj+ncat-1) 380 continue 370 continue 360 continue 350 continue ! do 336 k=1,na ! write(*,*)(vgam(k,j),j=1,na) !c336 continue ! (b). Find the attached probability for each observation do j=1,ncat pf(j)=para(j) enddo DEallocate(para,dint,cp,der,der2,tder2) DEallocate(tp,der22,vtg,tder) exit endif 1000 CONTINUE end subroutine ORsimple !--|----------------------------------------| ! | inversion of positive definite matrix | ! | by cholesky decomposition | !--|----------------------------------------|------------------ ! | input a is the original matrix, output a is the inversion ! | the originalmatrix is not kept. ! | a(nxn) only a submatrix of mxm is inverted. !--------------------------------------------------------------- subroutine chol(a,b,n,m) integer n,m double precision a(n,n),b(n) integer i,j,k double precision sum b(1)=sqrt(a(1,1)) do 5 j=2,m a(j,1)=a(j,1)/b(1) 5 continue do 10 i=2,m do 20 j=i,m sum=a(i,j) do 30 k=i-1,1,-1 sum=sum-a(i,k)*a(j,k) if(i.eq.j) then if(sum.le.0) then write(*,*)j,j pause 'choldc failed' endif b(i)=sqrt(sum) else a(j,i)=sum/b(i) endif 30 continue 20 continue 10 continue do 40 i=1,m a(i,i)=1.0/b(i) do 50 j=i+1,m sum=0.0 do 60 k=i,j-1 sum=sum-a(j,k)*a(k,i) 60 continue a(j,i)=sum/b(j) 50 continue 40 continue do 70 j=1,m do 80 i=j,m sum=0.0 do 90 k=i,m sum=sum+a(k,i)*a(k,j) 90 continue if(i.eq.j) then a(i,i)=sum else a(i,j)=sum a(j,i)=a(i,j) endif 80 continue 70 continue end subroutine chol !--|----------------------------------------| ! | matrix inversion by lower and upper | ! | triangular decomposition, n<=500 | !--|----------------------------------------| ! a(n,n): input matrix, a is changed ! b(n,n): output matrix, inversion of a ! n: the dimension of the original matrix !--------------------------------------------- subroutine solve(a,b,nn,n) integer n,nn,NMAX double precision a(nn,nn),b(nn,nn),TINY parameter(NMAX=600,TINY=1.0e-20) integer i,imax,j,k,ndum,indx(NMAX) double precision aamax,dum,sum,vv(NMAX) do 20 i=1,n indx(i)=i aamax=0.0 do 10 j=1,n if(abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 10 continue if(abs(aamax)<=TINY) pause 'singular matrix' vv(i)=1.0/aamax 20 continue do 90 j=1,n do 40 i=1,j-1 sum=a(i,j) do 30 k=1,i-1 sum=sum-a(i,k)*a(k,j) 30 continue a(i,j)=sum 40 continue aamax=0.0 do 60 i=j,n sum=a(i,j) do 50 k=1,j-1 sum=sum-a(i,k)*a(k,j) 50 continue a(i,j)=sum dum=vv(i)*abs(sum) if(dum.ge.aamax) then imax=i aamax=dum endif 60 continue if(j.ne.imax) then do 70 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 70 continue vv(imax)=vv(j) ndum=indx(imax) indx(imax)=indx(j) indx(j)=ndum endif ! indx(j)=imax if(abs(a(j,j))8) then y=(7*exp(-x*x/2)+16*exp(-x*x*(2-sqrt(2.0)))+(7+3.1415926*x*x/4)*exp(-x*x))/30 pvaluenorm=0.5*(y/2+y*y/8+y*y*y/16) endif end function pvaluenorm !------------------------------------------------------------------------ ! Rounding Y to have 2*ncat categories ! Both case and control defined by yc have equal number of ! categories ncat !---------------------------------------------------------------------- subroutine Roundv(y,yc,n,p,ncat) integer n,p,ncat real y(n,p),yc(n,p) real mid0,mid1,range0,rang1,delta0,delta1 mid1=sum(y*yc)/sum(yc) range1=maxval(y*yc+mid1*(1-yc))-minval(y*yc+mid1*(1-yc)) delta1=range1/ncat mid0=sum(y*(1-yc))/sum(1-yc) range0=maxval(y*(1-yc)+mid0*yc)-minval(y*(1-yc)+mid0*yc) delta0=range0/ncat y=int(y*yc/delta1)*delta1+int(y*(1-yc)/delta0)*delta0 end subroutine roundv