******************************************* * GEE2loc.fit * * release 1 * ******************************************* * * * Program to fit the one locus model and * * the two locus model by solving GEEs * * * * Joanna Biernacka (08/04) * * * * Department of Public Health Sciences * * University of Toronto * * * * Based on the GENEFINDER program * * written by Chiu et al (Liang et al 2001)* ******************************************* ******************************************************************** * Note that in this program the 2-locus mean function * is parametrized in terms of delta2 = (C1,C2,tau1,tau2) * where C1 = E[S(tau1)|ASP] and C2 = E[S(tau2)|ASP]. * An alternative parameterization in terms of two "effect size" * parameters, C1* and C2*, is possible (see Biernacka et al. 2005). ******************************************************************** ********************************************************** * The main program to fit both models * ********************************************************** program fit parameter (mnloci=64,npair=15,mnfam=1000) c mnloci = maximum number of loci c npair = maximum number of sibpairs per family c mnfam = maximum number of families double precision delta1(2),info1(2,2),U1(2) double precision delta2(4),info2(4,4),U2(4) double precision inv1(2,2),tdelta1(2),tmp1,tmp2 double precision inv2(4,4),tdelta2(4),tmp3,tmp4 double precision covR1(2,2),sdC,sdT double precision covR2(4,4),sdC1,sdC2,sdTau1,sdTau2 double precision tinit1(2),tinit2(4),eps double precision st,posit double precision tmp integer count,msize,totitr1,totitr2 integer nloci,nfam,maxint character*64 file1,file2 double precision sumx,sumxsq,np,varS double precision sumsc,minsc,maxT,minT,maxC,minC double precision last(4,2),lastd(6,4),step double precision bestC,bestT,newC,newT double precision maxT1,maxT2,maxC1,maxC2 double precision minT1,minT2,minC1,minC2 double precision bestC1,bestT1,newC1,newT1 double precision bestC2,bestT2,newC2,newT2 integer conv1,conv2,new common st(mnfam,npair,mnloci-1),posit(mnloci-1) common nfam,nloci,msize(mnfam),eps,varS(mnloci-1) ************************************************************************** c parameter definitions: c delta1 = (C,tau) : the two parameters in the 1-locus model c delta2 = (C1,C2,tau1,tau2) : the four parameters in the 2-locus model c U1,U2 : score vectors for the 1-lcous and 2-locus models respectively c info1,info2 : information matrices for the 1-locus and 2-locus models c inv1,inv2 : inverses of info1 and info2, respectively c tinit1,tinit2 : initial values for delta1 and delta2, respectively c st(l,k,i) : IBD sharing at ith marker for the kth ASP in the lth family c posit(i) : position, in cM, of the ith marker c tdelta : change in the elements of delta at the last iteration c covR1 : robust covariance matrix for the 1-locus model parameter estimates c sdC, sdT: robust standard error estimates of C and T in the 1-locus model c covR2 : robust covariance matrix for the 2-locus model parameter estimates c sdC1, sdC2, sdT1, sdT2: robust standard error estimates c parameters in the 2-locus model c************************************************************************** c get parameters from standard input read(*, '(A)')file1 read(*, '(A)')file2 read(*, '(I)')nloci read(*, '(I)')nfam read(*,'(F)')eps read(*, '(I)')maxint read(*, '(F)')delta1(2) read(*, '(F)')delta1(1) read(*, '(F)')delta2(3) read(*, '(F)')delta2(4) read(*, '(F)')delta2(1) read(*, '(F)')delta2(2) do 100 k=1,(nloci-1) read(*, '(F)')posit(k) 100 continue c read in IBD sharing data from input files open(unit=11,file=file1,status='old') open(unit=10,file=file2,status='old') do 490 k=1, nfam read(11,*) ifam, msize(k) 490 continue do 499 l=1,nfam do 499 k=1,msize(l) do 499 i=1,nloci-1 read(10,*) l,k,i,st(l,k,i) 499 continue c calculate the empirical variance estimates c of IBD sharing at each marker do 550 i=1,nloci-1 sumx=0 sumxsq=0 np=0 do 549 j=1,nfam do 549 k=1,msize(j) sumx=sumx+st(j,k,i) sumxsq=sumxsq+st(j,k,i)**2 np=np+1 549 continue varS(i)=(sumxsq-sumx**2/np)/(np-1) 550 continue **************************************************** *********************** part 1 ******************* *************** fit the one-locus model ************ **************************************************** **************************************************** tinit1(1)=delta1(1) tinit1(2)=delta1(2) step=1.0 minsc=1000.0 conv1=0 count=0 600 count=count+1 call scinf1(delta1,info1,U1) call infinv1(info1,inv1) c if current values of parameters are in-bounds, c check if "better" than previous values if ((delta1(1).gt.-1.0).and.(delta1(1).lt.1.0).and. & (delta1(2).gt.posit(1)).and.(delta1(2).lt.posit(nloci-1))) then sumsc=abs(U1(1))+abs(U1(2)) if (sumsc .lt. minsc) then bestC=delta1(1) bestT=delta1(2) minsc=sumsc end if endif do 630 i=1,2 tdelta1(i)=delta1(i) 630 continue do 640 i=1,2 delta1(i)=delta1(i)+inv1(i,1)*U1(1)/step+inv1(i,2)*U1(2)/step 640 continue tmp1=abs(tdelta1(1)-delta1(1)) tmp2=abs(tdelta1(2)-delta1(2)) c keep parameter values from last 4 iterations of non-convergent cases if (count .gt. maxint-4) then last(count+4-maxint,1)=delta1(1) last(count+4-maxint,2)=delta1(2) end if c ******************************************************* c here start checking for convergence/max iterations, etc c ******************************************************* if (count .ge. maxint) then c try reducing step size and continue iterations if (step .lt. 2.0) then step=step+1.0 count=0 go to 600 end if c establish range of parameters in last 4 iterations maxT=max(last(1,2),last(2,2),last(3,2), & last(4,2)) minT=min(last(1,2),last(2,2),last(3,2), & last(4,2)) maxC=max(last(1,1),last(2,1),last(3,1), & last(4,1)) minC=min(last(1,1),last(2,1),last(3,1), & last(4,1)) c first check if T is out of bounds if ((delta1(2) .lt. posit(1)) .or. & (delta1(2) .gt. posit(nloci-1))) then conv1=4 write(6,4998) "1-locus model fit:" write(6,4998) "Tau: out of bounds" go to 9999 c check if "periodic" convergence else if ((last(4,2) .eq. last(2,2)) .and. & (last(3,2) .eq. last(1,2))) then conv1=2 write(6,4998) "1-locus model fit:" write(6,4998) "Periodic non-convergence." write(6,4998) "Grid-search for solution." c if not periodic... else if (maxT-minT .lt. 5) then conv1=3 write(6,4998) "1-locus model fit:" write(6,4998) "Non-convergence." write(6,4998) "Grid-search for solution." c if truely not convergent yet in bounds else conv1=5 write(6,4998) "1-locus model fit:" write(6,4998) "Non-convergence." write(6,4998) "No solution." go to 9999 end if c check if any parameters out of bounds and if so, put inside border c and continue iterations else if ((delta1(1) .gt. 1) .or. (delta1(1) .lt. -1) & .or. (delta1(2) .lt. posit(1)) .or. & (delta1(2) .gt. posit(nloci-1))) then if (delta1(1) .gt. 1.0) then delta1(1)=1.0 else if (delta1(1) .lt. -1.0) then delta1(1)=-1.0 end if if (delta1(2) .gt. posit(nloci-1)) then delta1(2)=posit(nloci-1) else if (delta1(2) .lt. posit(1)) then delta1(2)=posit(1) end if go to 600 c check convergence criteria else if ((tmp1 .gt. (1.0e-3)/step) .or. & (tmp2 .gt. (5.0e-3)/step)) then go to 600 else conv1=1 go to 999 end if c ************************************************************** c 2-dim grid search for (C,tau): c performed only if regular convergence to a solution c is not attained. c ************************************************************** 800 delta1(1)=minC-.025 802 delta1(1)=delta1(1)+.005 delta1(2)=minT-1.1 if (delta1(2).lt.posit(1)) then delta1(2)=posit(1)+.1 end if if (delta1(2).gt.posit(nloci-1)) then delta1(2)=posit(1)-.5 end if 805 delta1(2)=delta1(2)+.1 call scor1(delta1,U1) sumsc=abs(U1(1))+abs(U1(2)) if (sumsc .lt. minsc) then bestC=delta1(1) bestT=delta1(2) minsc=sumsc end if if (delta1(2) .gt. (posit(nloci-1)-.1)) goto 850 if (delta1(2) .lt. maxT+1) goto 805 if (delta1(1) .lt. maxC+.02) goto 802 850 new=0 delta1(1)=bestC-.0105 900 delta1(1)=delta1(1)+.0005 delta1(2)=bestT-1.0 905 delta1(2)=delta1(2)+.005 call scor1(delta1,U1) sumsc=abs(U1(1))+abs(U1(2)) if (sumsc .lt. minsc) then newC=delta1(1) newT=delta1(2) minsc=sumsc new=1 end if if (delta1(2) .lt. bestT+1.0) goto 905 if (delta1(1) .lt. bestC+.01) goto 900 if (new .eq. 0) then delta1(1)=bestC delta1(2)=bestT else delta1(1)=newC delta1(2)=newT end if if (minsc.ge.10.0) then conv1=6 write(6,4998) "Grid-search did not find a solution." goto 9999 end if 999 call covmat1(delta1,info1,U1,inv1,covR1) sdC=covR1(1,1)**0.5 sdT=covR1(2,2)**0.5 totitr1=(step-1.0)*maxint+count *********************************************** * Print results of fitting the 1-locus model * *********************************************** write(6,4998) "" write(6,4998) "Fitting the one locus model..." write(6,4998) "" write(6,4994) "intial value for C:",tinit1(1) write(6,4994) "intial value for tau:",tinit1(2) write(6,4995) "convergence code:",conv1 write(6,4995) "number of iterations:",totitr1 write(6,4998) "" write(6,4998) "Parameter estimates:" write(6,4994) "Estimate of C:",delta1(1) write(6,4994) "Estimate of Tau:",delta1(2) write(6,4998) "" write(6,4998) "Robust standard error estimates:" write(6,4994) "SD of estimate of C:",sdC write(6,4994) "SD of estimate of Tau:",sdT write(6,4998) "" c write(6,4998) "Score vector at solution" c write(6,4994) "First component",U1(1) c write(6,4994) "Second component",U1(2) **************************************************** *********************** part 2 ******************* *************** fit the two-locus model ************ **************************************************** **************************************************** tinit2(1)=delta2(1) tinit2(2)=delta2(2) tinit2(3)=delta2(3) tinit2(4)=delta2(4) step=1.0 minsc=1000.0 conv2=0 count=0 1600 count=count+1 call scinf2(delta2,info2,U2) call infinv2(info2,inv2) sumsc=abs(U2(1))+abs(U2(2))+abs(U2(3))+abs(U2(4)) if (sumsc .lt. minsc) then bestC1=delta2(1) bestT1=delta2(3) bestC2=delta2(2) bestT2=delta2(4) minsc=sumsc end if do 1630 i=1,4 tdelta2(i)=delta2(i) 1630 continue do 1640 i=1,4 delta2(i)=delta2(i)+(inv2(i,1)*U2(1)+inv2(i,2)*U2(2) & +inv2(i,3)*U2(3)+inv2(i,4)*U2(4))/step 1640 continue tmp1=abs(tdelta2(1)-delta2(1)) tmp2=abs(tdelta2(2)-delta2(2)) tmp3=abs(tdelta2(3)-delta2(3)) tmp4=abs(tdelta2(4)-delta2(4)) if ((count .lt. maxint) .and. (delta2(3) .gt. delta2(4))) then tmp=delta2(4) delta2(4)=delta2(3) delta2(3)=tmp tmp=delta2(2) delta2(2)=delta2(1) delta2(1)=tmp end if c keep last 6 iterations of non-convergent cases if ((step .gt. 2.0) .and. (count .gt. maxint-6)) then lastd(count+6-maxint,1)=delta2(1) lastd(count+6-maxint,2)=delta2(2) lastd(count+6-maxint,3)=delta2(3) lastd(count+6-maxint,4)=delta2(4) end if if (count .ge. maxint) then c try reducing step size and continue iterations if (step .lt. 3.0) then step=step+1.0 count=0 go to 1600 end if maxT1=max(lastd(1,3),lastd(2,3),lastd(3,3), & lastd(4,3),lastd(5,3),lastd(6,3)) minT1=min(lastd(1,3),lastd(2,3),lastd(3,3), & lastd(4,3),lastd(5,3),lastd(6,3)) maxC1=max(lastd(1,1),lastd(2,1),lastd(3,1), & lastd(4,1),lastd(5,1),lastd(6,1)) minC1=min(lastd(1,1),lastd(2,1),lastd(3,1), & lastd(4,1),lastd(5,1),lastd(6,1)) maxT2=max(lastd(1,4),lastd(2,4),lastd(3,4), & lastd(4,4),lastd(5,4),lastd(6,4)) minT2=min(lastd(1,4),lastd(2,4),lastd(3,4), & lastd(4,4),lastd(5,4),lastd(6,4)) maxC2=max(lastd(1,2),lastd(2,2),lastd(3,2), & lastd(4,2),lastd(5,2),lastd(6,2)) minC2=min(lastd(1,2),lastd(2,2),lastd(3,2), & lastd(4,2),lastd(5,2),lastd(6,2)) c First check if either T is out of bounds after max iterations if ((delta2(3) .lt. posit(1)) .or. & (delta2(3) .gt. posit(nloci-1)) .or. & (delta2(4) .lt. posit(1)) .or. & (delta2(4) .gt. posit(nloci-1))) then c also should still be able to do a score test <<<<<<<<<<<<<<< conv2=4 write(6,4998) "2-locus model fit:" write(6,4998) "At least one tau out of bounds." go to 9999 c then check if "periodic" convergence else if ((lastd(4,3) .eq. lastd(2,3)) .and. & (lastd(3,3) .eq. lastd(1,3)) .and. & (lastd(4,4) .eq. lastd(2,4)) .and. & (lastd(3,4) .eq. lastd(1,4))) then conv2=2 c if not periodic... else if ((maxT1-minT1 .lt. 5) .and. & (maxT2-minT2 .lt. 5)) then conv2=3 c if truely not convergent yet in bounds else c also should still be able to do a score test... conv2=5 write(6,4998) "2-locus model fit:" write(6,4998) "Non-convergence." write(6,4998) "No solution." go to 9999 end if else if ((delta2(1) .gt. 1) .or. (delta2(1) .lt. -1) .or. & (delta2(2) .gt. 1) .or. (delta2(2) .lt. -1) .or. & (delta2(3) .gt. posit(nloci-1)) .or. & (delta2(3) .lt. posit(1)) .or. & (delta2(4) .gt. posit(nloci-1)) .or. & (delta2(4) .lt. posit(1))) then if (delta2(1) .gt. 1.0) then delta2(1)=1.0 else if (delta2(1) .lt. -1.0) then delta2(1)=-1.0 end if if (delta2(2) .gt. 1.0) then delta2(2)=1.0 else if (delta2(2) .lt. -1.0) then delta2(2)=-1.0 end if if (delta2(3) .gt. posit(nloci-1)) then delta2(3)=posit(nloci-1) else if (delta2(3) .lt. posit(1)) then delta2(3)=posit(1) end if if (delta2(4) .gt. posit(nloci-1)) then delta2(4)=posit(nloci-1) else if (delta2(4) .lt. posit(1)) then delta2(4)=posit(1) end if go to 1600 else if ((tmp1 .gt. (1.0e-3)/step) .or. & (tmp2 .gt. (1.0e-3)/step) .or. (tmp3 .gt. (5.0e-3)/step) & .or. (tmp4 .gt. (5.0e-3)/step)) then go to 1600 else conv2=1 go to 1999 end if c ************************************************************** c 4-dim grid search for (C1,C2,tau1,tau2): c Performed only if regular convergence to a solution c is not attained. c ************************************************************** 1800 delta2(1)=minC1-.01 1801 delta2(1)=delta2(1)+.005 1802 delta2(2)=minC2-.01 1803 delta2(2)=delta2(2)+.005 delta2(3)=minT1-0.3 if (delta2(3).lt.posit(1)) then delta2(3)=posit(1) end if if (delta2(3).gt.posit(nloci-1)) then delta2(3)=posit(nloci-1)-.5 end if 1805 delta2(3)=delta2(3)+.1 delta2(4)=minT2-0.3 if (delta2(4).lt.posit(1)) then delta2(4)=posit(1) end if if (delta2(4).gt.posit(nloci-1)) then delta2(4)=posit(nloci-1)-.5 end if 1807 delta2(4)=delta2(4)+.1 if (delta2(3) .gt. posit(nloci-1)) goto 1850 if (delta2(4) .gt. posit(nloci-1)) goto 1850 call scor2(delta2,U2) sumsc=abs(U2(1))+abs(U2(2))+abs(U2(3))+abs(U2(4)) if (sumsc .lt. minsc) then bestC1=delta2(1) bestT1=delta2(3) bestC2=delta2(2) bestT2=delta2(4) minsc=sumsc end if if (delta2(4) .lt. maxT2+0.2) goto 1807 if (delta2(3) .lt. maxT1+0.2) goto 1805 if (delta2(2) .lt. maxC2+.005) goto 1803 if (delta2(1) .lt. maxC1+.005) goto 1801 1850 new=0 delta2(1)=bestC1-.004 1900 delta2(1)=delta2(1)+.001 delta2(2)=bestC2-.004 1902 delta2(2)=delta2(2)+.001 delta2(3)=bestT1-0.025 1905 delta2(3)=delta2(3)+.005 delta2(4)=bestT2-0.025 1907 delta2(4)=delta2(4)+.005 call scor2(delta2,U2) sumsc=abs(U2(1))+abs(U2(2))+abs(U2(3))+abs(U2(4)) if (sumsc .lt. minsc) then newC1=delta2(1) newT1=delta2(3) newC2=delta2(2) newT2=delta2(4) minsc=sumsc new=1 end if if (delta2(4) .lt. bestT2+0.02) goto 1907 if (delta2(3) .lt. bestT1+0.02) goto 1905 if (delta2(2) .lt. bestC2+.003) goto 1902 if (delta2(1) .lt. bestC1+.003) goto 1900 if (new .eq. 0) then delta2(1)=bestC1 delta2(2)=bestC2 delta2(3)=bestT1 delta2(4)=bestT2 else delta2(1)=newC1 delta2(2)=newC2 delta2(3)=newT1 delta2(4)=newT2 end if if (minsc.ge.10.0) then write(6,4998) "Grid-search did not find a solution." goto 9999 end if 1999 call covmat2(delta2,info2,U2,inv2,covR2) sdC1=covR2(1,1)**0.5 sdC2=covR2(2,2)**0.5 sdTau1=covR2(3,3)**0.5 sdTau2=covR2(4,4)**0.5 totitr2=(step-1.0)*maxint+count *********************************************** * Print results of fitting the 2-locus model * *********************************************** write(6,4998) "" write(6,4998) "Fitting the two locus model..." write(6,4998) "" write(6,4994) "intial value for C1:",tinit2(1) write(6,4994) "intial value for C2:",tinit2(2) write(6,4994) "intial value for tau1:",tinit2(3) write(6,4994) "intial value for tau2:",tinit2(4) write(6,4995) "convergence code:",conv2 write(6,4995) "number of iterations:",totitr2 write(6,4998) "" write(6,4998) "Parameter estimates:" write(6,4994) "Estimate of C1:",delta2(1) write(6,4994) "Estimate of C2:",delta2(2) write(6,4994) "Estimate of Tau1:",delta2(3) write(6,4994) "Estimate of Tau2:",delta2(4) write(6,4998) "" write(6,4998) "Robust standard error estimates:" write(6,4994) "SD of estimate of C1:",sdC1 write(6,4994) "SD of estimate of C2:",sdC2 write(6,4994) "SD of estimate of Tau1:",sdTau1 write(6,4994) "SD of estimate of Tau2:",sdTau2 write(6,4998) "" c write(6,4998) "Score vector at solution:" c write(6,4994) "First component",U2(1) c write(6,4994) "Second component",U2(2) c write(6,4994) "Third component",U2(3) c write(6,4994) "Fourth component",U2(4) *************end of 2-locus model fit******************************* *************************************************************** ********* print formats ************ *************************************************************** 4994 format(a25,f12.6) 4995 format(a25,i3) 4998 format(a35) 9999 continue end c*********************************************************************** c******************** subroutines for 1-locus model ******************** c*********************************************************************** subroutine scinf1(delta,info,score) parameter (mnloci=64,npair=15,mnfam=1000) double precision delta(2),info(2,2),score(2) double precision eps,varS,tmp double precision d,dd integer msize,nloci,nfam,npair double precision st,posit double precision score1(mnfam),score2(mnfam) common st(mnfam,npair,mnloci-1),posit(mnloci-1) common nfam,nloci,msize(mnfam),eps,varS(mnloci-1) do 1700 i=1,2 score(i)=0.d0 do 1699 j=1,2 info(i,j)=0.d0 1699 continue 1700 continue do 1820 i=1,nfam score1(i)=0.0 score2(i)=0.0 do 1810 k=1,msize(i) do 1810 j=1,nloci-1 if (posit(j) .gt. (delta(2)+eps)) then dd=-1 else if (posit(j) .lt. (delta(2)-eps)) then dd=1 else dd=-1*(posit(j)-delta(2))/eps end if if (abs(posit(j)-delta(2)) .gt. eps) then d=abs(posit(j)-delta(2)) else d=1/(2*eps)*(posit(j)-delta(2))**2+0.5*eps end if tmp=dexp(-.08*d)*(delta(1)-delta(1)**2-.5)+.5 info(1,1)=info(1,1)+dexp(-.04*d*2)/tmp info(1,2)=info(1,2)+dexp(-0.04*d)*delta(1) & *dexp(-.04*d)*(-.04)*dd/tmp info(2,1)=info(1,2) info(2,2)=info(2,2)+(delta(1)*dexp(-.04*d)* & .04*dd)**2/tmp score1(i)=score1(i)+dexp(-.04*d)/tmp* & (st(i,k,j)-1-delta(1)*dexp(-.04*d)) score2(i)=score2(i)+ & delta(1)*dexp(-.04*d)*(-.04)*dd/ & tmp*(st(i,k,j)-1-delta(1)*dexp(-.04*d)) 1810 continue score(1)=score1(i)+score(1) score(2)=score2(i)+score(2) 1820 continue end c******************************************************************* subroutine infinv1(info,inv) double precision info(2,2),inv(2,2),tmp tmp=info(1,1)*info(2,2)-info(1,2)**2 if (tmp .lt. 0.0000000001) then info(1,1)=info(1,1)+.00000001 info(2,2)=info(2,2)+.00000001 tmp=info(1,1)*info(2,2)-info(1,2)**2 end if inv(1,1)=1/tmp*info(2,2) inv(1,2)=1/tmp*info(1,2)*(-1.0) inv(2,1)=inv(1,2) inv(2,2)=1/tmp*info(1,1) return end c********************************************************************** subroutine scor1(delta,score) parameter (mnloci=64,npair=15,mnfam=1000) double precision delta(2),score(2) double precision eps,varS,tmp double precision st,posit double precision d,dd double precision score1(mnfam), score2(mnfam) integer msize,nloci,nfam,npair common st(mnfam,npair,mnloci-1),posit(mnloci-1) common nfam,nloci,msize(mnfam),eps,varS(mnloci-1) do 2700 i=1,2 score(i)=0.d0 2700 continue do 2820 i=1,nfam score1(i)=0.0 score2(i)=0.0 do 2810 k=1,msize(i) do 2810 j=1,nloci-1 if (posit(j) .gt. (delta(2)+eps)) then dd=-1.0 else if (posit(j) .lt. (delta(2)-eps)) then dd=1.0 else dd=-1*(posit(j)-delta(2))/eps end if if (abs(posit(j)-delta(2)) .gt. eps) then d=abs(posit(j)-delta(2)) else d=1/(2*eps)*(posit(j)-delta(2))**2+0.5*eps end if tmp=dexp(-.08*d)*(delta(1)-delta(1)**2-.5)+.5 score1(i)=score1(i)+dexp(-.04*d)/tmp* & (st(i,k,j)-1-delta(1)*dexp(-.04*d)) score2(i)=score2(i)+ & delta(1)*dexp(-.04*d)*(-.04)*dd/ & tmp*(st(i,k,j)-1-delta(1)*dexp(-.04*d)) 2810 continue score(1)=score1(i)+score(1) score(2)=score2(i)+score(2) 2820 continue end c****************************************************************************** subroutine covmat1(delta,info,score,inv,cov) parameter (mnloci=64,npair=15,mnfam=1000) double precision delta(2),info(2,2),score(2) double precision eps,varS,tmp double precision d,dd integer msize,nloci,nfam double precision st,posit double precision score1(mnfam),score2(mnfam) double precision cov(2,2),inv(2,2),cov2(2,2) common st(mnfam,npair,mnloci-1),posit(mnloci-1) common nfam,nloci,msize(mnfam),eps,varS(mnloci-1) c first find score and info do 1700 i=1,2 score(i)=0.d0 do 1699 j=1,2 info(i,j)=0.d0 1699 continue 1700 continue do 1820 i=1,nfam score1(i)=0.0 score2(i)=0.0 do 1810 k=1,msize(i) do 1810 j=1,nloci-1 if (posit(j) .gt. (delta(2)+eps)) then dd=-1 else if (posit(j) .lt. (delta(2)-eps)) then dd=1 else dd=-1*(posit(j)-delta(2))/eps end if if (abs(posit(j)-delta(2)) .gt. eps) then d=abs(posit(j)-delta(2)) else d=1/(2*eps)*(posit(j)-delta(2))**2+0.5*eps end if tmp=dexp(-.08*d)*(delta(1)-delta(1)**2-.5)+.5 info(1,1)=info(1,1)+dexp(-.04*d*2)/tmp info(1,2)=info(1,2)+dexp(-0.04*d)*delta(1) & *dexp(-.04*d)*(-.04)*dd/tmp info(2,1)=info(1,2) info(2,2)=info(2,2)+(delta(1)*dexp(-.04*d)* & .04*dd)**2/tmp score1(i)=score1(i)+dexp(-.04*d)/tmp* & (st(i,k,j)-1-delta(1)*dexp(-.04*d)) score2(i)=score2(i)+ & delta(1)*dexp(-.04*d)*(-.04)*dd/ & tmp*(st(i,k,j)-1-delta(1)*dexp(-.04*d)) 1810 continue score(1)=score1(i)+score(1) score(2)=score2(i)+score(2) 1820 continue c invert the info matrix call infinv1(info,inv) c The covariance matrix for the estimators do 3877 i=1,2 do 3877 j=1,2 cov(i,j)=0.d0 cov2(i,j)=0.d0 3877 continue do 3880 l=1,nfam cov2(1,1)=cov2(1,1)+score1(l)*score1(l) cov2(1,2)=cov2(1,2)+score1(l)*score2(l) cov2(2,1)=cov2(2,1)+score2(l)*score1(l) cov2(2,2)=cov2(2,2)+score2(l)*score2(l) 3880 continue do 3870 i1=1,2 do 3870 j1=1,2 do 3870 i2=1,2 do 3870 j2=1,2 tmp=inv(i1,j2)*cov2(j2,i2)*inv(i2,j1) cov(i1,j1)=cov(i1,j1)+tmp 3870 continue end c*********************************************************************** c******************** subroutines for 2-locus model ******************** c*********************************************************************** subroutine scinf2(delta,info,score) parameter (mnloci=64,npair=15,mnfam=1000) double precision delta(4),info(4,4),score(4) double precision mu,dmudC1,dmudC2,dmudT1,dmudT2 double precision eps,varS integer msize,nloci,nfam double precision st,posit double precision score1(mnfam),score2(mnfam) double precision score3(mnfam),score4(mnfam) common st(mnfam,npair,mnloci-1),posit(mnloci-1) common nfam,nloci,msize(mnfam),eps,varS(mnloci-1) do 1700 i=1,4 score(i)=0.d0 do 1699 j=1,4 info(i,j)=0.d0 1699 continue 1700 continue do 1820 i=1,nfam score1(i)=0.0 score2(i)=0.0 score3(i)=0.0 score4(i)=0.0 do 1810 k=1,msize(i) do 1810 j=1,nloci-1 if (posit(j) .lt. delta(3)) then mu=1+dexp(.04*(posit(j)-delta(3)))*delta(1) dmudC1=dexp(.04*posit(j)-.04*delta(3)) dmudC2=0.0 dmudT1=(-.04)*dexp(.04*posit(j)-.04*delta(3))*delta(1) dmudT2=0.0 else if (posit(j) .gt. delta(4)) then mu=1+dexp(-.04*(posit(j)-delta(4)))*delta(2) dmudC1=0.0 dmudC2=dexp(-.04*posit(j)+.04*delta(4)) dmudT1=0.0 dmudT2=.04*dexp(-.04*posit(j)+.04*delta(4))*delta(2) else mu=1+(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))*(delta(1)+1)+ & (dexp(.04*posit(j)-.04*delta(4))- & dexp(-.04*posit(j)-.04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))*(delta(2)+1)- & (dexp(-.04*posit(j)+.04*delta(3))+ & dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)-.04*delta(4))) dmudC1=(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)- & .08*delta(4)))/(1-dexp(.08*delta(3)-.08*delta(4))) dmudC2=(dexp(.04*posit(j)-.04*delta(4))- & dexp(-.04*posit(j)-.04*delta(4)+ & .08*delta(3)))/(1-dexp(.08*delta(3)-.08*delta(4))) dmudT1=(.04*dexp(-.04*posit(j)+.04*delta(3))- & .04*dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(1)+1)+.08*(dexp(-.04*posit(j)+ & .04*delta(3))-dexp(.04*posit(j)+ & .04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(1)+1)*dexp(.08*delta(3)-.08*delta(4))- & .08*dexp(-.04*posit(j)-.04*delta(4)+ & .08*delta(3))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(2)+1)+.08*(dexp(.04*posit(j)- & .04*delta(4))-dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(2)+1)*dexp(.08*delta(3)- & .08*delta(4))-.04*dexp(-.04*posit(j)+ & .04*delta(3))/ & (1+dexp(.04*delta(3)-.04*delta(4)))+ & .04*(dexp(-.04*posit(j)+ & .04*delta(3))+dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)- & .04*delta(4)))**2*dexp(.04*delta(3)-.04*delta(4)) dmudT2=.08*dexp(.04*posit(j)+.04*delta(3)- & .08*delta(4))/(1-dexp(.08*delta(3)- & .08*delta(4)))*(delta(1)+1)- & .08*(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(1)+1)*dexp(.08*delta(3)-.08*delta(4))+ & (-.04*dexp(.04*posit(j)-.04*delta(4))+ & .04*dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(2)+1)-.08*(dexp(.04*posit(j)- & .04*delta(4))-dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(2)+1)*dexp(.08*delta(3)-.08*delta(4))+ & .04*dexp(.04*posit(j)-.04*delta(4))/ & (1+dexp(.04*delta(3)-.04*delta(4)))- & .04*(dexp(-.04*posit(j)+.04*delta(3))+ & dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)-.04*delta(4)))**2* & dexp(.04*delta(3)-.04*delta(4)) end if info(1,1)=info(1,1)+dmudC1*dmudC1/varS(j) info(1,2)=info(1,2)+dmudC1*dmudC2/varS(j) info(1,3)=info(1,3)+dmudC1*dmudT1/varS(j) info(1,4)=info(1,4)+dmudC1*dmudT2/varS(j) info(2,2)=info(2,2)+dmudC2*dmudC2/varS(j) info(2,3)=info(2,3)+dmudC2*dmudT1/varS(j) info(2,4)=info(2,4)+dmudC2*dmudT2/varS(j) info(3,3)=info(3,3)+dmudT1*dmudT1/varS(j) info(3,4)=info(3,4)+dmudT1*dmudT2/varS(j) info(4,4)=info(4,4)+dmudT2*dmudT2/varS(j) score1(i)=score1(i)+dmudC1*(st(i,k,j)-mu)/varS(j) score2(i)=score2(i)+dmudC2*(st(i,k,j)-mu)/varS(j) score3(i)=score3(i)+dmudT1*(st(i,k,j)-mu)/varS(j) score4(i)=score4(i)+dmudT2*(st(i,k,j)-mu)/varS(j) 1810 continue score(1)=score1(i)+score(1) score(2)=score2(i)+score(2) score(3)=score3(i)+score(3) score(4)=score4(i)+score(4) 1820 continue end c****************************************************************************** subroutine infinv2(info,inv) double precision info(4,4),inv(4,4),tmp double precision a11,a12,a13,a14,a22,a23,a24,a33,a34,a44 double precision b11,b12,b13,b14,b22,b23,b24,b33,b34,b44 a11=info(1,1) a12=info(1,2) a13=info(1,3) a14=info(1,4) a22=info(2,2) a23=info(2,3) a24=info(2,4) a33=info(3,3) a34=info(3,4) a44=info(4,4) tmp=(a11*a22*a33*a44-a11*a22*a34**2-a11*a23**2*a44+ & 2*a11*a23*a24*a34-a11*a24**2*a33-a12**2*a33*a44+ & a12**2*a34**2+2*a12*a23*a13*a44- & 2*a12*a23*a14*a34-2*a12*a24*a13*a34+2*a12*a24*a14*a33- & a22*a13**2*a44+2*a13*a22*a14*a34+a13**2*a24**2- & 2*a13*a24*a14*a23-a22*a14**2*a33+a14**2*a23**2) 1825 if (tmp .lt. 0.00000001) then a11=a11+.000001 a22=a22+.000001 a33=a33+.000001 a44=a44+.000001 tmp=(a11*a22*a33*a44-a11*a22*a34**2-a11*a23**2*a44+ & 2*a11*a23*a24*a34-a11*a24**2*a33-a12**2*a33*a44+ & a12**2*a34**2+2*a12*a23*a13*a44- & 2*a12*a23*a14*a34-2*a12*a24*a13*a34+2*a12*a24*a14*a33- & a22*a13**2*a44+2*a13*a22*a14*a34+a13**2*a24**2- & 2*a13*a24*a14*a23-a22*a14**2*a33+a14**2*a23**2) end if if (tmp .lt. 0.00000001) goto 1825 b11 = (a22*a33*a44-a22*a34**2-a23**2*a44+ & 2*a23*a24*a34-a24**2*a33)/tmp b12 = -(a12*a33*a44-a12*a34**2-a23*a13*a44+ & a23*a14*a34+a24*a13*a34-a24*a14*a33)/tmp b13 = (a12*a23*a44-a12*a24*a34-a22*a13*a44+ & a22*a14*a34+a13*a24**2-a24*a14*a23)/tmp b14 = -(a12*a23*a34-a12*a24*a33-a22*a13*a34+ & a22*a14*a33+a23*a13*a24-a14*a23**2)/tmp b22 = (a11*a33*a44-a11*a34**2-a13**2*a44+ & 2*a13*a14*a34-a14**2*a33)/tmp b23 = -(a11*a23*a44-a11*a24*a34-a13*a12*a44+ & a13*a14*a24+a14*a12*a34-a14**2*a23)/tmp b24 = (a11*a23*a34-a11*a24*a33-a13*a12*a34+ & a13**2*a24+a14*a12*a33-a14*a13*a23)/tmp b33 = (a11*a22*a44-a11*a24**2-a12**2*a44+ & 2*a12*a14*a24-a14**2*a22)/tmp b34 = -(a11*a22*a34-a11*a24*a23-a12**2*a34+ & a12*a13*a24+a14*a12*a23-a14*a13*a22)/tmp b44 = (a11*a22*a33-a11*a23**2-a12**2*a33+ & 2*a12*a13*a23-a13**2*a22)/tmp inv(1,1)=b11 inv(1,2)=b12 inv(1,3)=b13 inv(1,4)=b14 inv(2,1)=b12 inv(2,2)=b22 inv(2,3)=b23 inv(2,4)=b24 inv(3,1)=b13 inv(3,2)=b23 inv(3,3)=b33 inv(3,4)=b34 inv(4,1)=b14 inv(4,2)=b24 inv(4,3)=b34 inv(4,4)=b44 end c*************************************************************************** subroutine scor2(delta,score) parameter (mnloci=64,npair=15,mnfam=1000) double precision delta(4),score(4) double precision mu,dmudC1,dmudC2,dmudT1,dmudT2 double precision eps,varS integer msize,nloci,nfam double precision st,posit double precision score1(mnfam),score2(mnfam) double precision score3(mnfam),score4(mnfam) common st(mnfam,npair,mnloci-1),posit(mnloci-1) common nfam,nloci,msize(mnfam),eps,varS(mnloci-1) do 2700 i=1,4 score(i)=0.d0 2700 continue do 2820 i=1,nfam score1(i)=0.0 score2(i)=0.0 score3(i)=0.0 score4(i)=0.0 do 2810 k=1,msize(i) do 2810 j=1,nloci-1 if (posit(j) .lt. delta(3)) then mu=1+dexp(.04*(posit(j)-delta(3)))*delta(1) dmudC1=dexp(.04*posit(j)-.04*delta(3)) dmudC2=0.0 dmudT1=(-.04)*dexp(.04*posit(j)-.04*delta(3))*delta(1) dmudT2=0.0 else if (posit(j) .gt. delta(4)) then mu=1+dexp(-.04*(posit(j)-delta(4)))*delta(2) dmudC1=0.0 dmudC2=dexp(-.04*posit(j)+.04*delta(4)) dmudT1=0.0 dmudT2=.04*dexp(-.04*posit(j)+.04*delta(4))*delta(2) else mu=1+(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))*(delta(1)+1)+ & (dexp(.04*posit(j)-.04*delta(4))- & dexp(-.04*posit(j)-.04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))*(delta(2)+1)- & (dexp(-.04*posit(j)+.04*delta(3))+ & dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)-.04*delta(4))) dmudC1=(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)- & .08*delta(4)))/(1-dexp(.08*delta(3)-.08*delta(4))) dmudC2=(dexp(.04*posit(j)-.04*delta(4))- & dexp(-.04*posit(j)-.04*delta(4)+ & .08*delta(3)))/(1-dexp(.08*delta(3)-.08*delta(4))) dmudT1=(.04*dexp(-.04*posit(j)+.04*delta(3))- & .04*dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(1)+1)+.08*(dexp(-.04*posit(j)+ & .04*delta(3))-dexp(.04*posit(j)+ & .04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(1)+1)*dexp(.08*delta(3)-.08*delta(4))- & .08*dexp(-.04*posit(j)-.04*delta(4)+ & .08*delta(3))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(2)+1)+.08*(dexp(.04*posit(j)- & .04*delta(4))-dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(2)+1)*dexp(.08*delta(3)- & .08*delta(4))-.04*dexp(-.04*posit(j)+ & .04*delta(3))/ & (1+dexp(.04*delta(3)-.04*delta(4)))+ & .04*(dexp(-.04*posit(j)+ & .04*delta(3))+dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)- & .04*delta(4)))**2*dexp(.04*delta(3)-.04*delta(4)) dmudT2=.08*dexp(.04*posit(j)+.04*delta(3)- & .08*delta(4))/(1-dexp(.08*delta(3)- & .08*delta(4)))*(delta(1)+1)- & .08*(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(1)+1)*dexp(.08*delta(3)-.08*delta(4))+ & (-.04*dexp(.04*posit(j)-.04*delta(4))+ & .04*dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(2)+1)-.08*(dexp(.04*posit(j)- & .04*delta(4))-dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(2)+1)*dexp(.08*delta(3)-.08*delta(4))+ & .04*dexp(.04*posit(j)-.04*delta(4))/ & (1+dexp(.04*delta(3)-.04*delta(4)))- & .04*(dexp(-.04*posit(j)+.04*delta(3))+ & dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)-.04*delta(4)))**2* & dexp(.04*delta(3)-.04*delta(4)) end if score1(i)=score1(i)+dmudC1*(st(i,k,j)-mu)/varS(j) score2(i)=score2(i)+dmudC2*(st(i,k,j)-mu)/varS(j) score3(i)=score3(i)+dmudT1*(st(i,k,j)-mu)/varS(j) score4(i)=score4(i)+dmudT2*(st(i,k,j)-mu)/varS(j) 2810 continue score(1)=score1(i)+score(1) score(2)=score2(i)+score(2) score(3)=score3(i)+score(3) score(4)=score4(i)+score(4) 2820 continue end c****************************************************************************** subroutine covmat2(delta,info,score,inv,cov) parameter (mnloci=64,npair=15,mnfam=1000) double precision delta(4),info(4,4),score(4) double precision mu,dmudC1,dmudC2,dmudT1,dmudT2 double precision eps,varS integer msize,nloci,nfam double precision st,posit double precision score1(mnfam),score2(mnfam) double precision score3(mnfam),score4(mnfam) double precision cov(4,4),tmp,inv(4,4),cov2(4,4) common st(mnfam,npair,mnloci-1),posit(mnloci-1) common nfam,nloci,msize(mnfam),eps,varS(mnloci-1) c first find score and info do 1700 i=1,4 score(i)=0.d0 do 1699 j=1,4 info(i,j)=0.d0 1699 continue 1700 continue do 1820 i=1,nfam score1(i)=0.0 score2(i)=0.0 score3(i)=0.0 score4(i)=0.0 do 1810 k=1,msize(i) do 1810 j=1,nloci-1 if (posit(j) .lt. delta(3)) then mu=1+dexp(.04*(posit(j)-delta(3)))*delta(1) dmudC1=dexp(.04*posit(j)-.04*delta(3)) dmudC2=0.0 dmudT1=(-.04)*dexp(.04*posit(j)-.04*delta(3))*delta(1) dmudT2=0.0 else if (posit(j) .gt. delta(4)) then mu=1+dexp(-.04*(posit(j)-delta(4)))*delta(2) dmudC1=0.0 dmudC2=dexp(-.04*posit(j)+.04*delta(4)) dmudT1=0.0 dmudT2=.04*dexp(-.04*posit(j)+.04*delta(4))*delta(2) else mu=1+(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))*(delta(1)+1)+ & (dexp(.04*posit(j)-.04*delta(4))- & dexp(-.04*posit(j)-.04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))*(delta(2)+1)- & (dexp(-.04*posit(j)+.04*delta(3))+ & dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)-.04*delta(4))) dmudC1=(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)- & .08*delta(4)))/(1-dexp(.08*delta(3)-.08*delta(4))) dmudC2=(dexp(.04*posit(j)-.04*delta(4))- & dexp(-.04*posit(j)-.04*delta(4)+ & .08*delta(3)))/(1-dexp(.08*delta(3)-.08*delta(4))) dmudT1=(.04*dexp(-.04*posit(j)+.04*delta(3))- & .04*dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(1)+1)+.08*(dexp(-.04*posit(j)+ & .04*delta(3))-dexp(.04*posit(j)+ & .04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(1)+1)*dexp(.08*delta(3)-.08*delta(4))- & .08*dexp(-.04*posit(j)-.04*delta(4)+ & .08*delta(3))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(2)+1)+.08*(dexp(.04*posit(j)- & .04*delta(4))-dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(2)+1)*dexp(.08*delta(3)- & .08*delta(4))-.04*dexp(-.04*posit(j)+ & .04*delta(3))/ & (1+dexp(.04*delta(3)-.04*delta(4)))+ & .04*(dexp(-.04*posit(j)+ & .04*delta(3))+dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)- & .04*delta(4)))**2*dexp(.04*delta(3)-.04*delta(4)) dmudT2=.08*dexp(.04*posit(j)+.04*delta(3)- & .08*delta(4))/(1-dexp(.08*delta(3)- & .08*delta(4)))*(delta(1)+1)- & .08*(dexp(-.04*posit(j)+.04*delta(3))- & dexp(.04*posit(j)+.04*delta(3)-.08*delta(4)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(1)+1)*dexp(.08*delta(3)-.08*delta(4))+ & (-.04*dexp(.04*posit(j)-.04*delta(4))+ & .04*dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))* & (delta(2)+1)-.08*(dexp(.04*posit(j)- & .04*delta(4))-dexp(-.04*posit(j)- & .04*delta(4)+.08*delta(3)))/ & (1-dexp(.08*delta(3)-.08*delta(4)))**2* & (delta(2)+1)*dexp(.08*delta(3)-.08*delta(4))+ & .04*dexp(.04*posit(j)-.04*delta(4))/ & (1+dexp(.04*delta(3)-.04*delta(4)))- & .04*(dexp(-.04*posit(j)+.04*delta(3))+ & dexp(.04*posit(j)-.04*delta(4)))/ & (1+dexp(.04*delta(3)-.04*delta(4)))**2* & dexp(.04*delta(3)-.04*delta(4)) end if info(1,1)=info(1,1)+dmudC1*dmudC1/varS(j) info(1,2)=info(1,2)+dmudC1*dmudC2/varS(j) info(1,3)=info(1,3)+dmudC1*dmudT1/varS(j) info(1,4)=info(1,4)+dmudC1*dmudT2/varS(j) info(2,2)=info(2,2)+dmudC2*dmudC2/varS(j) info(2,3)=info(2,3)+dmudC2*dmudT1/varS(j) info(2,4)=info(2,4)+dmudC2*dmudT2/varS(j) info(3,3)=info(3,3)+dmudT1*dmudT1/varS(j) info(3,4)=info(3,4)+dmudT1*dmudT2/varS(j) info(4,4)=info(4,4)+dmudT2*dmudT2/varS(j) score1(i)=score1(i)+dmudC1*(st(i,k,j)-mu)/varS(j) score2(i)=score2(i)+dmudC2*(st(i,k,j)-mu)/varS(j) score3(i)=score3(i)+dmudT1*(st(i,k,j)-mu)/varS(j) score4(i)=score4(i)+dmudT2*(st(i,k,j)-mu)/varS(j) 1810 continue score(1)=score1(i)+score(1) score(2)=score2(i)+score(2) score(3)=score3(i)+score(3) score(4)=score4(i)+score(4) 1820 continue c invert the info matrix call infinv2(info,inv) c The covariance matrix for the estimators do 3877 i=1,4 do 3877 j=1,4 cov(i,j)=0.d0 cov2(i,j)=0.d0 3877 continue do 3880 l=1,nfam cov2(1,1)=cov2(1,1)+score1(l)*score1(l) cov2(1,2)=cov2(1,2)+score1(l)*score2(l) cov2(1,3)=cov2(1,3)+score1(l)*score3(l) cov2(1,4)=cov2(1,4)+score1(l)*score4(l) cov2(2,1)=cov2(2,1)+score2(l)*score1(l) cov2(2,2)=cov2(2,2)+score2(l)*score2(l) cov2(2,3)=cov2(2,3)+score2(l)*score3(l) cov2(2,4)=cov2(2,4)+score2(l)*score4(l) cov2(3,1)=cov2(3,1)+score3(l)*score1(l) cov2(3,2)=cov2(3,2)+score3(l)*score2(l) cov2(3,3)=cov2(3,3)+score3(l)*score3(l) cov2(3,4)=cov2(3,4)+score3(l)*score4(l) cov2(4,1)=cov2(4,1)+score4(l)*score1(l) cov2(4,2)=cov2(4,2)+score4(l)*score2(l) cov2(4,3)=cov2(4,3)+score4(l)*score3(l) cov2(4,4)=cov2(4,4)+score4(l)*score4(l) 3880 continue do 3870 i1=1,4 do 3870 j1=1,4 do 3870 i2=1,4 do 3870 j2=1,4 tmp=inv(i1,j2)*cov2(j2,i2)*inv(i2,j1) cov(i1,j1)=cov(i1,j1)+tmp 3870 continue end c******************** extra subroutines******************************* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c routine to invert a matrix by Gaussian elimination c Ainv=inverse(A) A(n,n) Ainv(n,n) c Note: sizes maxed at 64 x 64 - re dimension for larger matrices; c D matrix is extended. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine invert(A,Ainv,n) double precision A(64,64),Ainv(64,64) double precision D(64,128) c initialize the reduction matrix n2 = 2*n do 1 i = 1,n do 2 j = 1,n D(i,j) = A(i,j) d(i,n+j) = 0. 2 continue D(i,n+i) = 1. 1 continue c do the reduction do 3 i = 1,n alpha = D(i,i) if(alpha .eq. 0.) go to 300 do 4 j = 1,n2 D(i,j) = D(i,j)/alpha 4 continue do 5 k = 1,n if((k-i).eq.0) go to 5 beta = D(k,i) do 6 j = 1,n2 D(k,j) = D(k,j) - beta*D(i,j) 6 continue 5 continue 3 continue c copy result into output matrix do 7 i = 1,n do 8 j = 1,n Ainv(i,j) = D(i,j+n) 8 continue 7 continue return 300 print *,'*** ERROR: Singular matrix ***' return end c****************************************************************************** c***************************************************************************