!------------------------------------------------------------------- ! Fit a joint conditional odds ratio model without missing data ! Pre-specify groups of variables to be modeled ! Allow for zero coefficients !-------------------------------------------------------------------- program main parameter(TINY=1e-7) character*500 DataSetName !hold the name of the data file integer n,nn,np,niter,qa,qb,maxcat real eps integer unita,unitb,misscode integer nfield,Yfield,Xfield,Gfield real,allocatable:: dat(:) !Temperoarily store a record of the data. integer,Allocatable:: Efield(:),MD(:) real,Allocatable:: y(:),x(:),G(:),E(:,:) real,Allocatable:: Wa(:,:),WB(:,:) open(unit=5,file='ModelOR',status='old') !modeling information read(5,*)DataSetName !file name that stores the data read(5,*)niter,eps,maxcat read(5,*)np,nfield read(5,*)Yfield,Xfield,Gfield Allocate(dat(nfield),Efield(np)) read(5,*)(Efield(i),i=1,np) !locations of Z in the data file read(5,*) Misscode ! Code for missing values 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(y(nn),x(nn),G(nn),E(nn,np)) !read and allocate data open(unit=3,file=DataSetName,status='old') do i=1,nn read(3,*)(dat(j),j=1,nfield) y(i)=dat(Yfield) x(i)=dat(Xfield) do j=1,np E(i,j)=dat(Efield(j)) enddo enddo close(unit=3) ! Reduce the number of categories for quantitative trait call Round(y,nn,maxcat) call Round(x,nn,maxcat) !do j=1,np ! call Round(E(:,j),nn,maxcat) !enddo qa=3 qb=np allocate(WA(nn,qa),WB(nn,qb),md(qa+1)) do k=1,qa+1 md(k)=k enddo unita=7 unitb=8 open(unit=unita, file='estimate',status='unknown') open(unit=unitb, file='stderr',status='unknown') do nsnp=1,nfield-Gfield+1 ! number of SNPs write(*,*)nsnp,nsnp,nsnp open(unit=3,file=DataSetName,status='old') n=0 do i=1,nn !excluding sample with missing genotype read(3,*)(dat(j),j=1,nfield) if(abs(dat(Gfield+nsnp-1)-misscode)>10e-6 & .or.abs(Y(i)-misscode)>10e-6.or.abs(x(i)-misscode)>10e-6) then !Assume SNPs are listed consecutively n=n+1 WA(n,1)=Y(i) WA(n,2)=X(i) WA(n,3)=dat(Gfield+nsnp-1) do j=1,np WB(n,j)=E(i,j) enddo endif enddo ! Odds Ratio Model !Call JorML(WA,WB,n,qa,qb,md,qa,unita,unitb) Call JorMLnoi(WA(1:n,:),WB(1:n,:),n,qa,qb,md,qa,unita,unitb) ! no interaction model close(unit=3) enddo close(unit=unita) close(unit=unitb) !call resultsummary() deallocate(WA,WB,md) end program Main !------------------------------------------------------------- ! Summary of results from the secondary trait analysis !------------------------------------------------------------- !----------------------------------------------------------------------------------------------------- ! Maximum likelihood approach ! for the semiparametric joint odds ratio model of p(y|x) ! in the form of ! eta(y1;(y2,\cdots,yg)|x)...\eta(y(g-1);yg|y10,...,y(g-2)0,x) ! eta(y1;x|y20,...,yg0)...eta(yg;x|y10,...,y(g-1)0) ! p(y1|y20,...,yg0,x0)...p(yg|y10,...,y(g-1)0,x0) ! divided by a normalized constant, where ! y1 has dimension md(1) ! ...... ! yg has dimension md(g) ! mg is the number of Y variable groups modeled and md are the dimensions of each variable group ! md(1)+...+md(mg)=p ! We model each odds ratio function to have interactions up to the second order. !---------------------------------------------------------------------------------------------------- Subroutine JorML(Y,X,n,p,q,md,mg,unita,unitb) integer niter double precision eps parameter(niter=50,eps=1e-5) integer,intent(in):: n,p,q,mg,unita,unitb integer,intent(in):: md(mg+1) real,intent(in):: Y(n,p),X(n,q) real dy(p),dx(q) !to store y0,x0 integer ngam(mg) !dimension of parameters in the odds ratio model for eta(yk|y(k+1),...yg,x) integer ncat(mg) !number of categories in each modeling group integer indx(n,mg),cat(n+1,mg),Rcat(n+1,mg),loc(mg) !indx indicates the order addresses !cat indicates the starting position of each category from 1 to ncat in ordered sequence !Rcat indicates the category that an ordered observation has. real ycat(n,p) !denote the categories of y double precision prob, cprob, lkh double precision,allocatable:: lam(:,:) ! lam for the based function integer,allocatable::parind(:) !model specification integer,allocatable:: map(:) ! subset covariate matrix for output double precision,allocatable:: par(:),delta(:) double precision,allocatable:: tder(:),tempder(:),vder(:,:) double precision,allocatable:: der(:),der2(:,:) !!!!create order for variable group Call orderFD(y,n,p,md,mg,ncat,indx,cat,Rcat) do k=1,mg do j=1,ncat(k) do kk=md(k),md(k+1)-1 ycat(j,kk)=y(indx(cat(j,k),k),kk) enddo !write(*,*)(k, ycat(j,kk),kk=md(k),md(k+1)-1) enddo enddo nt=1 do k=1,mg nt=nt*ncat(k) enddo write(*,*)'Combinatorial number is:',nt if(nt>10000) then write(*,*)'Warning: Combinatorial number is too large. To reduce computation time,' write(*,*) 'rounding maybe considered for the continuous variables' endif nm=0 ! For the total number of regression parameters do k=1,mg ngam(k)=(md(k+1)-md(k))*(p+1-md(k+1)+q) !linear terms do j=k+1,mg !interaction terms for yyy and yyx ngam(k)=ngam(k)+(md(k+1)-md(k))*(md(j+1)-md(j))*(p+1-md(j+1)+q) enddo ngam(k)=ngam(k)+(md(k+1)-md(k))*q*(q-1)/2 !interactions for yxx nm=nm+ngam(k)+ncat(k)-1 enddo !odds ratio parameters and nonparametric parameters write(*,*)'The number of parameters is:',nm !!!!======================================================== allocate(parind(nm),map(nm),par(nm),lam(n,mg),delta(nm)) allocate(tder(nm),tempder(nm),vder(nm,nm)) allocate(der(nm),der2(nm,nm)) !!!!======================================================== !!initialize the parameters do k=1,nm par(k)=0 parind(k)=1 !all (linear and interaction terms) included first, and then adjust enddo lam=0 !Readjust model specification itemp1=1 !staring location itemp2=0 ! ending location ntotal=1 ! location of variance matrix to be output do k=1,mg itemp2=itemp2+(md(k+1)-md(k))*(p+1-md(k+1)+q) !linear terms parind(itemp1:itemp2)=1 !retain linear terms !write(*,*)itemp1,itemp2 itemp1=itemp2+1 do j=k+1,mg !interaction terms for yyy and yyx itemp2=itemp2+(md(k+1)-md(k))*(md(j+1)-md(j))*(p+1-md(j+1)+q) enddo itemp2=itemp2+(md(k+1)-md(k))*q*(q-1)/2 !interactions for yxx !write(*,*)itemp1,itemp2 parind(itemp1:itemp2)=0 !exclude interaction terms itemp1=itemp2+1 !update the starting point itemp2=itemp2+ncat(k)-1 !write(*,*)itemp1,itemp2 parind(itemp1:itemp2)=1 !always retain the nonparametric parameters itemp1=itemp2+1 !update the starting point enddo !complete adjustment of model specification ! itemp=0 ! do k=1,mg ! if(k==1) then ! itemp=itemp+(md(k+1)-md(k))*(p+1-md(k+1)+q)+1 !the first interaction is included ! do j=(md(k+1)-md(k))*(p+1-md(k+1)+q)+2,ngam(k) ! itemp=itemp+1 ! parind(itemp)=0 ! no other interaction is included ! enddo ! else ! itemp=itemp+(md(k+1)-md(k))*(p+1-md(k+1)+q) ! do j=(md(k+1)-md(k))*(p+1-md(k+1)+q)+1,ngam(k) ! itemp=itemp+1 ! parind(itemp)=0 !no interaction is included ! enddo ! endif ! itemp=itemp+ncat(k)-1 ! enddo !assign y0,x0 !dy=y(1,:) !dx=x(1,:) !do k=1,p ! dy(k)=sum(y(:,k))/n !enddo !do k=1,q ! dx(k)=sum(x(:,k))/n !enddo dy=0 dx=0 NR:Do Irep=1,NITER !!compute the first and second derivatives der=0 der2=0 itemp=0 do k=1,mg ! one group at a time itemp=itemp+ngam(k) do j=1,ncat(k)-1 ! first lambda itemp=itemp+1 lam(j,k)=par(itemp) enddo enddo lkh=0.0 subj:DO i=1,n cprob=0 !denomator of the conditional expectation do k=1,nm tempder(k)=0.0 do l=k,nm vder(k,l)=0.0 !temp. collect cross-product of first der. vder(l,k)=vder(k,l) enddo enddo CPCE:DO jt=1,nt !find conditional probabilities and conditional expectations itemp=jt-1 do k=1,mg loc(k)=1+itemp-int(itemp/ncat(k))*ncat(k) itemp=int(itemp/ncat(k)) enddo prob=0 tder=0.0 itemp=0 do k=1,mg do kk=md(k),md(k+1)-1 !It depends on how the parameters are stored! do j=k+1,mg do jj=md(j),md(j+1)-1 !linear terms in y itemp=itemp+1 tder(itemp)=parind(itemp)*(ycat(loc(k),kk)-dy(kk))*(ycat(loc(j),jj)-dy(jj)) !(y(indx(cat(loc(k),k),k),kk)-dy(kk)) & ! *(y(indx(cat(loc(j),j),j),jj)-dy(jj)) prob=prob+par(itemp)*tder(itemp) enddo enddo do j=1,q !linear term in x itemp=itemp+1 tder(itemp)=parind(itemp)*(ycat(loc(k),kk)-dy(kk))*(x(i,j)-dx(j)) !(y(indx(cat(loc(k),k),k),kk)-dy(kk))*(x(i,j)-dx(j)) prob=prob+par(itemp)*tder(itemp) enddo do j=k+1,mg !interactions in yyy do jj=md(j),md(j+1)-1 do l=j+1,mg do ll=md(l),md(l+1)-1 itemp=itemp+1 tder(itemp)=parind(itemp)*(ycat(loc(k),kk)-dy(kk)) & *(ycat(loc(j),jj)*ycat(loc(l),ll)-dy(jj)*dy(ll)) !(y(indx(cat(loc(k),k),k),kk)-dy(kk)) & !*(y(indx(cat(loc(j),j),j),jj)*y(indx(cat(loc(l),l),l),ll)-dy(jj)*dy(ll)) prob=prob+par(itemp)*tder(itemp) enddo enddo do l=1,q ! iteraction yyx itemp=itemp+1 tder(itemp)=parind(itemp)*(ycat(loc(k),kk)-dy(kk)) & *(ycat(loc(j),jj)*x(i,l)-dy(jj)*dx(l)) !(y(indx(cat(loc(k),k),k),kk)-dy(kk)) & !*(y(indx(cat(loc(j),j),j),jj)*x(i,l)-dy(jj)*dx(l)) prob=prob+par(itemp)*tder(itemp) enddo enddo enddo do j=1,q !interaction yxx do jj=j+1,q itemp=itemp+1 tder(itemp)=parind(itemp)*(ycat(loc(k),kk)-dy(kk))*(x(i,j)*x(i,jj)-dx(j)*dx(jj)) !(y(indx(cat(loc(k),k),k),kk)-dy(kk))*(x(i,j)*x(i,jj)-dx(j)*dx(jj)) prob=prob+par(itemp)*tder(itemp) enddo enddo enddo do j=1,ncat(k)-1 !conditional densities at a fixed condition itemp=itemp+1 if(loc(k)==j) then tder(itemp)=1.0 prob=prob+lam(loc(k),k) endif enddo enddo prob=exp(prob) do j=1,nm tempder(j)=tempder(j)+tder(j)*prob vder(j,j)=vder(j,j)+tder(j)*tder(j)*prob do l=j+1,nm vder(j,l)=vder(j,l)+tder(j)*tder(l)*prob enddo enddo cprob=cprob+prob enddo CPCE lkh=lkh-log(cprob) der=der-tempder/cprob !level 3: minus sign because the expectation is subtracted from the score do j=1,nm der2(j,j)=der2(j,j)+vder(j,j)/cprob-tempder(j)*tempder(j)/cprob/cprob do l=j+1,nm der2(j,l)=der2(j,l)+vder(j,l)/cprob-tempder(j)*tempder(l)/cprob/cprob enddo enddo !compute contributions from the observed data itemp=0 do k=1,mg do kk=md(k),md(k+1)-1 !It depends on how the parameters are stored! do j=k+1,mg do jj=md(j),md(j+1)-1 !linear terms in y itemp=itemp+1 der(itemp)=der(itemp)+parind(itemp)*(y(i,kk)-dy(kk))*(y(i,jj)-dy(jj)) lkh=lkh+par(itemp)*parind(itemp)*(y(i,kk)-dy(kk))*(y(i,jj)-dy(jj)) enddo enddo do j=1,q !linear term in x itemp=itemp+1 der(itemp)=der(itemp)+parind(itemp)*(y(i,kk)-dy(kk))*(x(i,j)-dx(j)) lkh=lkh+par(itemp)*parind(itemp)*(y(i,kk)-dy(kk))*(x(i,j)-dx(j)) enddo do j=k+1,mg !interactions yyx do jj=md(j),md(j+1)-1 do l=j+1,mg do ll=md(l),md(l+1)-1 itemp=itemp+1 der(itemp)=der(itemp)+parind(itemp)*(y(i,kk)-dy(kk)) & *(y(i,jj)*y(i,ll)-dy(jj)*dy(ll)) lkh=lkh+par(itemp)*parind(itemp)*(y(i,kk)-dy(kk)) & *(y(i,jj)*y(i,ll)-dy(jj)*dy(ll)) enddo enddo do l=1,q itemp=itemp+1 der(itemp)=der(itemp)+parind(itemp)*(y(i,kk)-dy(kk)) & *(y(i,jj)*x(i,l)-dy(jj)*dx(l)) lkh=lkh+par(itemp)*parind(itemp)*(y(i,kk)-dy(kk)) & *(y(i,jj)*x(i,l)-dy(jj)*dx(l)) enddo enddo enddo do j=1,q ! interactions yxx do jj=j+1,q itemp=itemp+1 der(itemp)=der(itemp)+parind(itemp)*(y(i,kk)-dy(kk))*(x(i,j)*x(i,jj)-dx(j)*dx(jj)) lkh=lkh+par(itemp)*parind(itemp)*(y(i,kk)-dy(kk))*(x(i,j)*x(i,jj)-dx(j)*dx(jj)) enddo enddo enddo do j=1,ncat(k)-1 !conditional densities at a fixed condition itemp=itemp+1 if(Rcat(indx(i,k),k)==j) then der(itemp)=der(itemp)+1.0 lkh=lkh+lam(j,k) endif enddo enddo enddo subj !1.invert information matrix ! write(*,12)(der(j),J=1,NM) ! do j=1,nm ! do l=j+1,nm ! der2(l,j)=der2(j,l) ! enddo ! write(*,12)(der2(j,l),l=1,j) ! 12 Format(10(F10.3,2x)) ! enddo !pause Call chol(der2,parind,nm,'yes') !2.Newton-Raphson update do k=1,nm delta(k)=0 if(parind(k)==1) then do j=1,nm delta(k)=delta(k)+der2(k,j)*der(j)*parind(j) enddo endif enddo if(irep<5) delta=delta/4 par=par+delta write(*,*) 'Iteration=',irep write(*,*) 'Convergence=',sum(abs(delta)),'LKH=',lkh if(abs(sum(delta))<=eps) then ! write(*,*)'Parameter estimates: exact zero means that term was not included in the model.' ! write(*,61)par*parind ! do k=1,mg ! write(*,*)'linear parameter:', 'k=',k ! write(*,*) par( ! ngam(k)=(md(k+1)-md(k))*(p+1-md(k+1)+q) !linear terms ! do j=k+1,mg-1 !interaction terms among y. no interaction with x ! ngam(k)=ngam(k)+(md(k+1)-md(k))*(md(j+1)-md(j))*(p+1-md(j+1)) ! enddo ! nm=nm+ngam(k) ! nm=nm+ncat(k)-1 ! enddo !odds ratio parameters and nonparametric parameters 61 format(1x, 10(f10.5,2x)) ! write(*,*)'Standard error: exact zero means that term was not included in the model.' ! write(*,61)(sqrt(parind(j)*der2(j,j)),j=1,nm) ! write(*,*)'Variance-Covariance for parameter estimates:' ! do k=1,nm ! if(parind(k)==1) then ! write(*,61)(parind(k)*parind(j)*der2(k,j),j=1,nm) ! endif ! enddo exit NR endif enddo NR itemp1=1 !starting location itemp2=0 !ending location ntotal=0 do k=1,mg itemp2=itemp2+(md(k+1)-md(k))*(p+1-md(k+1)+q) !linear terms write(*,*)'Linear terms in model ',k write(*,61)(par(j)*parind(j),j=itemp1,itemp2) do i=itemp1,itemp2 if(parind(i)==1) then ntotal=ntotal+1 map(ntotal)=i endif enddo itemp1=itemp2+1 !updating ending location do j=k+1,mg !interaction terms for yyy and yyx itemp2=itemp2+(md(k+1)-md(k))*(md(j+1)-md(j))*(p+1-md(j+1)+q) do i=itemp1,itemp2 if(parind(i)==1) then ntotal=ntotal+1 map(ntotal)=i endif enddo itemp1=itemp2+1 !updating ending location enddo itemp2=itemp2+(md(k+1)-md(k))*q*(q-1)/2 !interactions for yxx do i=itemp1,itemp2 if(parind(i)==1) then ntotal=ntotal+1 map(ntotal)=i endif enddo itemp1=itemp2+1 !updating ending location itemp2=itemp2+ncat(k)-1 itemp1=itemp2+1 enddo !odds ratio parameters and nonparametric parameters !write(*,61)(par(map(j)),j=1,ntotal) write(unita,61)(par(map(j)),j=1,ntotal) do i=1,ntotal !write(*,61)(der2(map(i),map(j)),j=1,i) write(unitb,61)(der2(map(i),map(j)),j=1,i) enddo ! write(*,61)(par(j)*parind(j),j=1,5) ! write(*,61)(sqrt(parind(j)*der2(j,j)),j=1,5) ! write(unita,61)(par(j)*parind(j),j=1,5) ! write(unitb,61)(sqrt(parind(j)*der2(j,j)),j=1,5) end subroutine JorML !----------------------------------------------------------------------------------------------------- ! Maximum likelihood approach ! for the semiparametric joint odds ratio model of p(y|x) ! in the form of ! eta(y1;(y2,\cdots,yg|x)...\eta(y(g-1);yg|y10,...,y(g-2)0,x) ! eta(y1;x|y20,...,yg0)...eta(yg;x|y10,...,y(g-1)0) ! p(y1|y20,...,yg0,x0)...p(yg|y10,...,y(g-1)0,x0) ! divided by a normalized constant, where ! y1 has dimension md(1) ! ...... ! yg has dimension md(g) ! mg is the number of Y variable groups modeled and md are the dimensions of each variable group ! md(1)+...+md(mg)=p ! NO INTERACTION IS ALLOWED IN THIS MODEL, see JorML FOR INTERACTIONS. !---------------------------------------------------------------------------------------------------- Subroutine JorMLnoi(Y,X,n,p,q,md,mg,unita,unitb) integer niter double precision eps parameter(niter=50,eps=1e-5) integer,intent(in):: n,p,q,mg,unita,unitb integer,intent(in):: md(mg+1) real,intent(in):: Y(n,p),X(n,q) real dy(p),dx(q) !to store y0,x0 integer ngam(mg) !dimension of parameters in the odds ratio model for eta(yk|y(k+1),...yg,x) integer ncat(mg) !number of categories in each modeling group integer indx(n,mg),cat(n+1,mg),Rcat(n+1,mg),loc(mg) !indx indicates the order addresses !cat indicates the starting position of each category from 1 to ncat in ordered sequence !Rcat indicates the category that an ordered observation has. real ycat(n,p) !denote the categories of y double precision prob, cprob, lkh double precision,allocatable:: lam(:,:) ! lam for the based function double precision,allocatable:: par(:),delta(:) double precision,allocatable:: tder(:),tempder(:),vder(:,:) double precision,allocatable:: der(:),der2(:,:) integer,allocatable:: nmap(:) integer itemp !!!!create order for variable group Call orderFD(y,n,p,md,mg,ncat,indx,cat,Rcat) do k=1,mg do j=1,ncat(k) do kk=md(k),md(k+1)-1 ycat(j,kk)=y(indx(cat(j,k),k),kk) enddo !write(*,*)(k, ycat(j,kk),kk=md(k),md(k+1)-1) enddo enddo nt=1 do k=1,mg nt=nt*ncat(k) enddo write(*,*)'Combinatorial number is:',nt if(nt>10000) then write(*,*)'Warning: Combinatorial number is too large. To reduce computation time,' write(*,*) 'rounding maybe considered for the continuous variables' endif nm=0 ! For the total number of regression parameters do k=1,mg ngam(k)=(md(k+1)-md(k))*(p+1-md(k+1)+q) !linear terms nm=nm+ngam(k)+ncat(k)-1 !noparametric part enddo !odds ratio parameters and nonparametric parameters write(*,*)'The number of parameters is:',nm !!!!======================================================== allocate(par(nm),lam(n,mg),delta(nm),nmap(nm)) allocate(tder(nm),tempder(nm),vder(nm,nm)) allocate(der(nm),der2(nm,nm)) !!!!======================================================== !!initialize the parameters par=0 lam=0 dy=0 dx=0 NR:Do Irep=1,NITER !!compute the first and second derivatives der=0 der2=0 itemp=0 do k=1,mg ! one group at a time itemp=itemp+ngam(k) do j=1,ncat(k)-1 ! first lambda itemp=itemp+1 lam(j,k)=par(itemp) enddo enddo lkh=0.0 subj:DO i=1,n cprob=0 !denomator of the conditional expectation do k=1,nm tempder(k)=0.0 do l=k,nm vder(k,l)=0.0 !temp. collect cross-product of first der. vder(l,k)=vder(k,l) enddo enddo CPCE:DO jt=1,nt !find conditional probabilities and conditional expectations itemp=jt-1 do k=1,mg loc(k)=1+itemp-int(itemp/ncat(k))*ncat(k) itemp=int(itemp/ncat(k)) enddo prob=0 tder=0.0 itemp=0 do k=1,mg do kk=md(k),md(k+1)-1 !It depends on how the parameters are stored! do j=k+1,mg do jj=md(j),md(j+1)-1 !linear terms in y itemp=itemp+1 tder(itemp)=(ycat(loc(k),kk)-dy(kk))*(ycat(loc(j),jj)-dy(jj)) !(y(indx(cat(loc(k),k),k),kk)-dy(kk)) & ! *(y(indx(cat(loc(j),j),j),jj)-dy(jj)) prob=prob+par(itemp)*tder(itemp) enddo enddo do j=1,q !linear term in x itemp=itemp+1 tder(itemp)=(ycat(loc(k),kk)-dy(kk))*(x(i,j)-dx(j)) !(y(indx(cat(loc(k),k),k),kk)-dy(kk))*(x(i,j)-dx(j)) prob=prob+par(itemp)*tder(itemp) enddo enddo do j=1,ncat(k)-1 !conditional densities at a fixed condition itemp=itemp+1 if(loc(k)==j) then tder(itemp)=1.0 prob=prob+lam(loc(k),k) endif enddo enddo prob=exp(prob) do j=1,nm tempder(j)=tempder(j)+tder(j)*prob vder(j,j)=vder(j,j)+tder(j)*tder(j)*prob do l=j+1,nm vder(j,l)=vder(j,l)+tder(j)*tder(l)*prob enddo enddo cprob=cprob+prob enddo CPCE lkh=lkh-log(cprob) der=der-tempder/cprob !level 3: minus sign because the expectation is subtracted from the score do j=1,nm der2(j,j)=der2(j,j)+vder(j,j)/cprob-tempder(j)*tempder(j)/cprob/cprob do l=j+1,nm der2(j,l)=der2(j,l)+vder(j,l)/cprob-tempder(j)*tempder(l)/cprob/cprob enddo enddo !compute contributions from the observed data itemp=0 do k=1,mg do kk=md(k),md(k+1)-1 !It depends on how the parameters are stored! do j=k+1,mg do jj=md(j),md(j+1)-1 !linear terms in y itemp=itemp+1 der(itemp)=der(itemp)+(y(i,kk)-dy(kk))*(y(i,jj)-dy(jj)) lkh=lkh+par(itemp)*(y(i,kk)-dy(kk))*(y(i,jj)-dy(jj)) enddo enddo do j=1,q !linear term in x itemp=itemp+1 der(itemp)=der(itemp)+(y(i,kk)-dy(kk))*(x(i,j)-dx(j)) lkh=lkh+par(itemp)*(y(i,kk)-dy(kk))*(x(i,j)-dx(j)) enddo enddo do j=1,ncat(k)-1 !conditional densities at a fixed condition itemp=itemp+1 if(Rcat(indx(i,k),k)==j) then der(itemp)=der(itemp)+1.0 lkh=lkh+lam(j,k) endif enddo enddo enddo subj nmap=1 !temporarily used as indicator for submatrix inversion Call chol(der2,nmap,nm,'yes') !2.Newton-Raphson update do k=1,nm delta(k)=0 do j=1,nm delta(k)=delta(k)+der2(k,j)*der(j) enddo enddo if(irep<5) delta=delta/4 par=par+delta write(*,*) 'Iteration=',irep write(*,*) 'Convergence=',sum(abs(delta)),'LKH=',lkh if(abs(sum(delta))<=eps) then 61 format(1x, 10(f10.5,2x)) exit NR endif enddo NR ntotal=0 itemp=0 do k=1,mg write(*,*)'Parameter estimates for model: ',k write(*,61)(par(j),j=itemp+1,itemp+(md(k+1)-md(k))*(p+1-md(k+1)+q)) do j=1,(md(k+1)-md(k))*(p+1-md(k+1)+q) ntotal=ntotal+1 nmap(ntotal)=itemp+j enddo itemp=itemp+ncat(k)-1+(md(k+1)-md(k))*(p+1-md(k+1)+q) enddo write(unita,61)(par(nmap(j)),j=1,ntotal) do i=1,ntotal !write(*,61)(der2(map(i),map(j)),j=1,i) write(unitb,61)(der2(nmap(i),nmap(j)),j=1,i) enddo end subroutine JorMLnoi !---|---------------------| ! | Sorting subroutine | !---|---------------------|---------------------------------------------------------- ! x(n,p): input sequence containing missing datacoded by 'miscode' ! missing values are listed in the end ! indx(n,mg): output of address order of the sequence for each variable group ! mg: the number variable groups ! ncat(mg): number of distinctive values in each variable group ! md(mg): starting position of each variable group !-------------------------------------------------------------------------------------- subroutine orderFD(x,n,p,md,mg,ncat,indx,cat,Rcat) parameter(TINY=1e-6) integer n,p,mg,indtemp integer p1 integer indx(n,mg),cat(n+1,mg),Rcat(n+1,mg) integer md(mg+1),ncat(mg) real x(n,p) real temp do 1000 ii=1,mg p1=md(ii) do 10 i=1,n indx(i,ii)=i 10 continue !-------Sort the first variable in the group with the same missing pattern do 60 i=1,n-1 temp=x(indx(i,ii),p1) do 70 j=(i+1),n if(temp.GT.x(indx(j,ii),p1))then temp=x(indx(j,ii),p1) indtemp=indx(i,ii) indx(i,ii)=indx(j,ii) indx(j,ii)=indtemp else endif 70 continue 60 continue !-------Sort the remaining in order and keeping the previous order do while(p1.lt.(md(ii+1)-1)) temp=x(indx(1,ii),p1) itemp=1 !length of the segment itempd=0 !indicator for further ordering do 110 i=2,n if(abs(temp-x(indx(i,ii),p1)).LT.TINY) then itemp=itemp+1 if(i.eq.n) then do 120 j=(i-itemp+1),i temp=x(indx(j,ii),(p1+1)) do 130 l=(j+1),i if(temp.gt.x(indx(l,ii),(p1+1)))then temp=x(indx(l,ii),(p1+1)) indtemp=indx(l,ii) indx(l,ii)=indx(j,ii) indx(j,ii)=indtemp endif 130 continue 120 continue endif else if(itemp.gt.1)then itempd=itempd+1 do 160 j=(i-itemp),(i-1) temp=x(indx(j,ii),(p1+1)) do 170 l=(j+1),(i-1) if(temp.gt.x(indx(l,ii),(p1+1)))then temp=x(indx(l,ii),(p1+1)) indtemp=indx(l,ii) indx(l,ii)=indx(j,ii) indx(j,ii)=indtemp else endif 170 continue 160 continue endif itemp=1 temp=x(indx(i,ii),p1) endif 110 continue if(itempd.eq.0) goto 200 p1=p1+1 enddo ! collect number of categories and catuency from sequence x(,q1) ! missing does not count as a category, frequency the starting points of ! new values(is cumulative.freq+1) 200 temp=x(indx(1,ii),p1) !c itemp=1 ncat(ii)=1 cat(1,ii)=1 do 150 i=2,n if(abs(temp-x(indx(i,ii),p1)).GT.TINY)then !c itemp=itemp+1 ncat(ii)=ncat(ii)+1 cat(ncat(ii),ii)=i temp=x(indx(i,ii),p1) endif 150 continue cat(ncat(ii)+1,ii)=n+1 ! cumulative frequency do 220 j=1,n do 230 l=1,ncat(ii) if(cat(l,ii).le.j.and.j.lt.cat(l+1,ii)) then Rcat(indx(j,ii),ii)=l goto 220 endif 230 continue 220 continue 1000 CONTINUE end subroutine orderFD !------------------------------------------------------------- ! | inversion of a submatrix of a positive definite matrix | ! | by cholesky decomposition | !--|---------------------------------------------------------|------ ! | input a is the original matrix, output a is the decomposition ! | or the inversion of the original matrix, which is not kept. ! | b(n) indicates which components are in the submatrix to be inverted. ! | only the submatrix correspond to nonzero b(n) is decomposed or inverted. !------------------------------------------------------------------------ subroutine chol(aa,b,n,inverse) integer,intent(in):: n,b(n) character*3 inverse double precision det double precision aa(n,n) integer i,j,k,m,itemp double precision summ double precision, allocatable:: a(:,:) m=sum(b) !do k=1,n ! write(*,*)(aa(k,j),j=1,n) ! 2 format(10(F10.3,2x)) !enddo !write(*,*)m !pause allocate(a(m,m)) !extract the submatrix itemp=0 do k=1,n if(b(k)==1) then itemp=itemp+1 a(itemp,itemp)=aa(k,k) itemp2=itemp do j=k+1,n if(b(j)==1) then itemp2=itemp2+1 a(itemp,itemp2)=aa(k,j) a(itemp2,itemp)=a(itemp,itemp2) endif enddo endif enddo !matrix decomposition a(1,1)=sqrt(a(1,1)) det=a(1,1) do j=2,m a(j,1)=a(j,1)/a(1,1) enddo do i=2,m do, j=i,m summ=a(i,j) do, k=i-1,1,-1 summ=summ-a(i,k)*a(j,k) if(i.eq.j) then if(summ<=0) then write(*,*) i,j,k,summ pause 'choldc failed' endif a(i,i)=sqrt(summ) else a(j,i)=summ/a(i,i) endif enddo enddo det=det*a(i,i) ! write(*,*)i,det enddo IF(INVERSE=='yes')THEN do i=1,m a(i,i)=1.0/a(i,i) do j=i+1,m summ=0.0 do k=i,j-1 summ=summ-a(j,k)*a(k,i) enddo a(j,i)=summ/a(j,j) enddo enddo do j=1,m do i=j,m summ=0.0 do k=i,m summ=summ+a(k,i)*a(k,j) enddo if(i.eq.j) then a(i,i)=summ else a(i,j)=summ a(j,i)=a(i,j) endif enddo enddo ENDIF !Insert back the submatrix itemp=0 do k=1,n if(b(k)==1) then itemp=itemp+1 aa(k,k)=a(itemp,itemp) itemp2=itemp do j=k+1,n if(b(j)==1) then itemp2=itemp2+1 aa(k,j)=a(itemp,itemp2) aa(j,k)=aa(k,j) endif enddo endif enddo Deallocate(a) end subroutine chol !-------------------------------------------- ! Generate normal random variable !-------------------------------------------- double precision function rnorm() double precision x, y, z; z=2.0 do while(z>=1.0) x = random() * 2.0 - 1.0 y = random() * 2.0 - 1.0 z = x * x + y * y enddo z = x * sqrt(- 2.0 * log(z) / z ) rnorm=z end function rnorm !------------------------------------------------------------------------ ! Rounding Y to have ncat categories ! !---------------------------------------------------------------------- subroutine Round(y,n,maxcat) integer n,maxcat real y(n),yr(n,1) integer index(n),ncat real delta yr(:,1)=y call order(yr,n,1,index) rtemp=y(index(1)) ! collect the numbers of categories ncat=1 do i=2,n if(abs(y(index(i))-rtemp).GT.1.0e-6) then rtemp=y(index(i)) ncat=ncat+1 endif enddo if(ncat>maxcat) then delta=(maxval(y)-minval(y))/maxcat y=delta*int(y/delta) endif end subroutine round !---|---------------------| ! | Sorting subroutine | !---|---------------------|--------------------------------------------------------- ! x(n,p): input sequence ! index(n): output of address order of sequence ! sect(nsect+1): starting positions of blocks with the same vales after sorting ! m: the maximum dimension of sect with m>=nsect+1 !------------------------------------------------------------------------------------ subroutine order(x,n,p,index) parameter(TINY=10e-6) integer index(n),indtemp,p real x(n,p) !integer x(n,p) !real temp integer temp !number of ordered blocks and their ending positions. integer nsect(p),sect(p,n+1) Initial: do i=1,n index(i)=i enddo Initial !sort the first component FirstSort:do i=1,(n-1) temp=x(index(i),1) do j=(i+1),n if(temp.GT.x(index(j),1))then temp=x(index(j),1) indtemp=index(i) index(i)=index(j) index(j)=indtemp else endif enddo enddo FirstSort ! count the number of blocks and set the starting position of the blocks nsect(1)=1 sect(1,1)=1 temp=x(index(1),1) do i=2,n if(X(index(i),1).GT.temp) then temp=X(index(i),1) nsect(1)=nsect(1)+1 sect(1,nsect(1))=i endif enddo nsect(1)=nsect(1)+1 !the additional fictitious block and its starting position (for convenience). sect(1,nsect(1))=n+1 !write(*,*)'Starting position of each missing data pattern up to time t:' !write(9,*)'Starting position of each missing data pattern up to time t:' j=1 !write(*,*)'Time=',j,(sect(1,k),k=1,nsect(1)) !write(9,*)'Time=',j,(sect(1,k),k=1,nsect(1)) SortRest: do j=2,p nsect(j)=1 sect(j,1)=1 !initialize nsect: do k=1,nsect(j-1)-1 i:do i=sect(j-1,k),sect(j-1,k+1)-1 !within each block temp=x(index(i),j) l:do l=(i+1),sect(j-1,k+1)-1 if(temp.gt.x(index(l),j)) then temp=x(index(l),j) indtemp=index(i) index(i)=index(l) index(l)=indtemp endif enddo l enddo i !end of inner block ordering !count the number of subblocks and store the starting position of the subblocks in temp vector temp=x(index(sect(j-1,k)),j) do i=sect(j-1,k),sect(j-1,k+1)-1 if(X(index(i),j).GT.temp) then temp=X(index(i),j) nsect(j)=nsect(j)+1 sect(j,nsect(j))=i endif enddo nsect(j)=nsect(j)+1 sect(j,nsect(j))=sect(j-1,k+1) enddo nsect !write(*,*)'Time=',j,(sect(j,k),k=1,nsect(j)) !write(9,*)'Time=',j,(sect(j,k),k=1,nsect(j)) enddo SortRest end subroutine order