m_qtlmap_analyse_modlin_cox
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_analyse_modlin_cox
DESCRIPTION
NOTES
SEE ALSO
current_chr
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
current_chr
DESCRIPTION
The current chromosome while the likelihood calculs
current_ic
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
current_ic
DESCRIPTION
The current trait to analyse
f0
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
f0
DESCRIPTION
value of the likelihood under H0
fm0
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
fm0
DESCRIPTION
value of the likelihood by full-sib family under H0
DIMENSIONS
np
fm1
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
fm1
DESCRIPTION
value of the likelihood by full-sib family under H1
DIMENSIONS
np
fp0
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
fp0
DESCRIPTION
value of the likelihood by sire family under H0
DIMENSIONS
np
fp1
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
fp1
DESCRIPTION
value of the likelihood by sire family under H1
DIMENSIONS
np
Hy
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
Hy
DESCRIPTION
DIMENSIONS
ntot
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
ntot
DESCRIPTION
posi_h1
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
posi_h1
DESCRIPTION
rangy
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
rangy
DESCRIPTION
DIMENSIONS
sig1
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
sig1
DESCRIPTION
The standart deviation under H0
DIMENSIONS
np
Ssy
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
Ssy
DESCRIPTION
DIMENSIONS
Sy
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
Sy
DESCRIPTION
DIMENSIONS
xmu1m
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
xmu1m
DESCRIPTION
The polygenic mean for each full sib family
DIMENSIONS
np
xmu1p
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
xmu1p
DESCRIPTION
The polygenic mean for each sire family under H0
DIMENSIONS
np
yr
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
fm1
DESCRIPTION
DIMENSIONS
yt
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Variables ]
NAME
yt
DESCRIPTION
DIMENSIONS
calcul_rang
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
calcul_rang
DESCRIPTION
Calcul du rang de y
INPUTS
ic : index of the trait
SOURCE
282 subroutine calcul_rang(ic) 283 integer , dimension(:),allocatable :: nryt 284 integer :: nm1,nm2,nd1,nd2, ifail,ic,kd,jm,ip, nt 285 integer :: stat 286 ! 287 allocate (nryt(nd),STAT=stat) 288 call check_allocate(stat,'nryt [m_qtlmap_analyse_modlin_cox]') 289 290 ntot=0 291 do ip=1,np 292 nm1=nmp(ip)+1; nm2=nmp(ip+1) 293 do jm=nm1,nm2 294 nd1=ndm(jm)+1; nd2=ndm(jm+1) 295 do kd=nd1,nd2 296 if(presentc(ic,kd)) then 297 ntot=ntot+1 298 yt(ntot)=y(ic,kd) 299 yr(ntot)=0.d0 300 nryt(ntot)=1 301 endif 302 enddo 303 enddo 304 enddo 305 306 ifail=0 307 CALL MATH_QTLMAP_M01DAF(yt,1,Ntot,'Descending',nryt,ifail) 308 ! write(nficout,*) (i,yt(i),nryt(i),i=1,nt) 309 nt=0 310 do ip=1,np 311 nm1=nmp(ip)+1; nm2=nmp(ip+1) 312 do jm=nm1,nm2 313 nd1=ndm(jm)+1; nd2=ndm(jm+1) 314 do kd=nd1,nd2 315 if(presentc(ic,kd)) then 316 nt=nt+1 317 rangy(kd)=nryt(nt) 318 yr(rangy(kd))=y(ic,kd) 319 ! write(nficout,*) ic,kd, nt, y(ic,kd), rangy(kd) 320 endif 321 enddo 322 enddo 323 enddo 324 325 deallocate (nryt) 326 327 end subroutine calcul_rang
end_analyse_modlin_cox
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
end_analyse_modlin_cox
DESCRIPTION
deallocation of solution/buffer arrays
NOTES
SOURCE
256 subroutine end_analyse_modlin_cox 257 deallocate (sig1) 258 deallocate (xmu1p) 259 deallocate (xmu1m) 260 deallocate (fp0) 261 !deallocate (fp1) 262 deallocate (fm0) 263 !deallocate (fm1) 264 deallocate (yt) 265 deallocate (yr) 266 deallocate (rangy) 267 ! deallocate (Hy) 268 ! deallocate (Sy) 269 ! deallocate (Ssy) 270 end subroutine end_analyse_modlin_cox
funct_0qtl_modlin_cox
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
funct_0qtl_modlin_cox
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H0
INPUTS
ic : index of the trait est_moy : general mean
NOTES
------------------- MODELE DE COX ET METHODE SANS DERIVES-----------------
Calcul du max de vraisemblance sous H1 et H0
Sous programme appele par e04jyf, optimiseur utilise par cherche
LE DENOMINATEUR DE l(kd,i) est approxime:
LE CALCUL DE LA VRAISEMBLANCE EST REALISE EN 3 phases:
Dans la 1ere boucle, on calcul pour chaque animal de Hy(kd) utilise dans la
3eme boucle et de Sy(rang(ic,kkd)) utilise dans la 2eme boucle, pour calculer
les contributions de chaque animal.
la 3eme boucle calcule la vraisemblance proprement dit.
SOURCE
495 subroutine funct_0qtl_modlin_cox(n,x,f,iuser,user) 496 use m_qtlmap_analyse_lin_gen, only : nbfem,nbef,ntnifix,nbco,covdir,nivdir,corniv,vecsol 497 498 integer , intent(in) :: n 499 ! 500 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 501 ! 502 real (kind=dp) ,dimension(n), intent(in) :: x 503 real (kind=dp) , intent(inout) :: f 504 integer , dimension(1), intent(in) :: iuser 505 real (kind=dp) ,dimension(1), intent(in) :: user 506 507 ! Tableaux dimensionnes selon nm, le nombre de meres 508 real (kind=dp) ,dimension(nm) :: effm 509 510 511 integer ::jnit,ip,nm1,nm2,jm,nd1,nd2,kkd,ilev,ief,ico,icar,iid,i 512 real (kind=dp) :: vmere,vpf,v,l,dv 513 ! 514 ! 515 ! INITIALISATION 516 icar = current_ic 517 if (np.lt.2.and.(count(estime(icar,:))).lt.2) jnit=0 518 if (np.ge.2.and.(count(estime(icar,:))).lt.2) jnit=1 519 if (np.ge.2.and.(count(estime(icar,:))).ge.2) jnit=2 520 521 f=0.d0 522 523 ! CALCULS pour chaque descendant du terme de participation � la contribution 524 Sy=0.d0 525 do ip=1,np 526 nm1=nmp(ip)+1 527 nm2=nmp(ip+1) 528 do jm=nm1,nm2 529 ! on ne consid�re que les m�res 530 nd1=ndm(jm)+1 531 nd2=ndm(jm+1) 532 do kkd=nd1,nd2 533 IF(presentc(icar,kkd)) then 534 v=1.d0 535 ! effet du pere 536 if (np.ge.2) then 537 if(vecsol(nivdir(kkd,1))) then 538 ilev=corniv(nivdir(kkd,1)) 539 v=v*x(ilev) 540 endif 541 end if 542 !effet de la m�re 543 if (nbfem.ge.2) then 544 if(estime(icar,jm)) then 545 if(vecsol(nivdir(kkd,2))) then 546 ilev=corniv(nivdir(kkd,2)) 547 v=v*x(ilev) 548 end if 549 end if 550 endif 551 !effet fixes 552 do ief=jnit+1,jnit+nbef 553 if(vecsol(nivdir(kkd,ief))) then 554 ilev=corniv(nivdir(kkd,ief)) 555 v=v*x(ilev) 556 end if 557 end do 558 ! effet des covariables 559 do ico=1,nbco 560 if(vecsol(ntnifix+ico)) then 561 ilev=corniv(ntnifix+ico) 562 v=v*(x(ilev)**covdir(kkd,ico)) 563 end if 564 end do 565 566 dv = (v) 567 ! if ( dv /= dv ) 568 if ( dv > huge(dv)) then 569 call log_mess('algorithm goes out of bound, iteration was left!!!!!!',DEBUG_DEF) 570 fm1=INIFINY_REAL_VALUE 571 fp1=INIFINY_REAL_VALUE 572 f=INIFINY_REAL_VALUE 573 return 574 end if 575 Sy(rangy(kkd))=dv 576 ENDIF !fin de la condition sur les perf 577 enddo !fin boucle sur kd 578 enddo !fin de la boucle sur jm 579 enddo !fin de la boucle sur ip 580 581 ! CALCUL DES CONTRIBUTIONS 582 do ip=1,np 583 fp0(ip)=0.d0 584 nm1=nmp(ip)+1 585 nm2=nmp(ip+1) 586 do jm=nm1,nm2 587 fm0(jm)=0.d0 588 nd1=ndm(jm)+1 589 nd2=ndm(jm+1) 590 do kkd=nd1,nd2 591 IF(presentc(icar,kkd)) then 592 Ssy(kkd)=0.d0 593 ind_sort : do iid=1,ntot 594 if (yr(rangy(kkd)).GT.yr(iid)) exit ind_sort 595 Ssy(kkd)=Ssy(kkd)+Sy(iid) 596 enddo ind_sort 597 IF ((ndelta(icar,kkd).EQ.1)) then 598 if (Sy(rangy(kkd)) == 0.d0.or.SSy(kkd)==0.d0) then 599 fm0(jm)=INIFINY_REAL_VALUE 600 fp0=INIFINY_REAL_VALUE 601 f=INIFINY_REAL_VALUE 602 call log_mess('algorithm goes out of bound, iteration was left!!!!!!',DEBUG_DEF) 603 return 604 else 605 ! print *,kkd, 'ssy:',SSy(kkd), 'sy:', sy(rangy(kkd)) 606 fm0(jm)= fm0(jm)-(dlog(Sy(rangy(kkd)))-dlog(SSy(kkd))) 607 !print *,'SOUS H0', ndelta(icar,kkd),kkd,Sy(rangy(kkd)),SSy(kkd), 'fm0', fm0(jm) 608 endif 609 ENDIF 610 ENDIF !fin de la condition sur les perf 611 ! write(nficout,*) icar,kkd, rangy(kkd),ndelta(icar,kkd), l 612 enddo !fin boucle sur kkd 613 fp0(ip)=fp0(ip)+fm0(jm) 614 f=f+fm0(jm) 615 !print *, 'FM0', fm0(jm) 616 enddo !fin boucle sur jm 617 !print *, 'FP0', fp0(ip) 618 enddo !fin boucle sur ip 619 620 ! print *,f 621 return 622 end subroutine funct_0qtl_modlin_cox
funct_1qtl_modlin_cox
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
funct_1qtl_modlin_cox
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H1
INPUTS
n : number of parameter
x : the parameter's values
iuser :
user :
OUTPUTS
f : the value of the likelihood
NOTES
------------------- MODELE DE COX ET METHODE SANS DERIVES-----------------
Calcul du max de vraisemblance sous H1
Sous programme appele par e04jyf, o
LE DENOMINATEUR DE l(kd,i) est approxim�:
LE CALCUL DE LA VRAISEMBLANCE EST REALISE EN 3 phases:
Dans la 1�boucle, on calcul pour chaque animal de Hy(kd) utilis� dans la
3� boucle et de Sy(rang(ic,kkd)) utilis� dans la 2�boucle, pour calculer
les contributions de chaque animal.
la 3�boucle calcule la vraisemblance proprement dit.
SOURCE
924 subroutine funct_1qtl_modlin_cox(n,x,f,iuser,user) 925 use m_qtlmap_analyse_lin_gen, only : nbfem,nbef,ntnifix,nbco,covdir,nivdir,corniv,vecsol 926 927 integer , intent(in) :: n 928 ! 929 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 930 ! 931 real (kind=dp) ,dimension(n), intent(in) :: x 932 real (kind=dp) , intent(inout) :: f 933 integer , dimension(1), intent(in) :: iuser 934 real (kind=dp) ,dimension(1), intent(in) :: user 935 936 real (kind=dp) :: my, gy,dvq 937 938 939 integer ::init,jnit,ip,nm1,nm2,jm,nd1,nd2,kd,kkd,ilev,ief,ico,icar,nt,iid,i,ig, ngeno1, ngeno2, nbdv 940 real (kind=dp) :: vmere,vpf,v, vq,l 941 ! 942 ! 943 fp1=0.d0 944 fm1=0.d0 945 f=0.d0 946 icar = current_ic 947 ! jnit: numero du 1° effet fixe, init: numero du 1° effet polygénique père 948 if (np==1.and.(count(estime(icar,:))).lt.2) jnit=2 949 if (np==1.and.(count(estime(icar,:)))==0) jnit=1 950 if (np==1.and.(count(estime(icar,:))).ge.2) jnit=3 951 init=3 952 if (np.ge.2.and.(count(estime(icar,:))).lt.2) jnit=3 953 if (np.ge.2.and.(count(estime(icar,:)))==0 ) then;jnit=2; init=2; endif 954 if (np.ge.2.and.(count(estime(icar,:))).ge.2) jnit=4 955 956 f=0.d0 957 Sy = 0.d0 958 ! CALCULS pour chaque descendant du terme de participation � la contribution 959 b_pe: do ip=1,np 960 nm1=nmp(ip)+1 961 nm2=nmp(ip+1) 962 b_me: do jm=nm1,nm2 963 ngeno1=ngenom(current_chr,jm)+1 964 ngeno2=ngenom(current_chr,jm+1) 965 b_gme: do ig=ngeno1,ngeno2 966 nd1=ngend(current_chr,ig)+1 967 nd2=ngend(current_chr,ig+1) 968 vpf=0.d0 969 b_indgm: do kd=nd1,nd2 970 Hy(kd)=0.d0 971 kkd=ndesc(current_chr,kd) 972 if(presentc(icar,kkd)) then 973 v=1.d0 974 ! effet du pere 975 if (np.ge.2) then 976 if(vecsol(nivdir(kkd,init))) then 977 ilev=corniv(nivdir(kkd,init)) 978 v=v*x(ilev) 979 end if 980 endif 981 !effet de la m�re 982 983 !print *,'nivdir(kkd,4):',nivdir(kkd,4),' s1:',size(vecsol) 984 985 if (nfem.ge.2) then 986 if(estime(icar,jm)) then 987 if(vecsol(nivdir(kkd,4))) then 988 ilev=corniv(nivdir(kkd,4)) 989 v=v*x(ilev) 990 end if 991 end if 992 endif 993 !effet fixes 994 do ief=jnit+1,jnit+nbef 995 if(vecsol(nivdir(kkd,ief))) then 996 ilev=corniv(nivdir(kkd,ief)) 997 v=v*x(ilev) 998 end if 999 end do 1000 ! effet des covariables 1001 do ico=1,nbco 1002 if(vecsol(ntnifix+ico)) then 1003 ilev=corniv(ntnifix+ico) 1004 v=v*(x(ilev))**(covdir(kkd,ico)) 1005 end if 1006 end do 1007 ! calcul des contribution en fonction des effets QTL pere et mere 1008 b_pdd : do i=1,4 1009 ! effet QTL du pere 1010 vq=v 1011 if(vecsol(nivdir(kkd,1))) then 1012 ilev=corniv(nivdir(kkd,1)) 1013 if (i==1.or.i==2) then 1014 vq=vq 1015 else 1016 vq=vq*x(ilev) 1017 endif 1018 end if 1019 !effet QTL de la m�re 1020 if(estime(icar,jm)) then 1021 if (vecsol(nivdir(kkd,2))) then 1022 ilev=corniv(nivdir(kkd,2)) 1023 if (i==1.or.i==3) then 1024 vq=vq 1025 else 1026 vq=vq*x(ilev) 1027 endif 1028 end if 1029 end if 1030 ! if (vq.GT.1000.d0 .or. vq.LT.-1000.d0) then 1031 ! print *, 'vq=',vq, 'Vq est trop GRAND!' 1032 ! stop 1033 ! endif 1034 dvq = (vq) 1035 !print *, 'kd=', kd, 'dvq=', dvq 1036 if ( dvq > huge(dvq)) then 1037 ! print *, 'OUT kd=', kd, 'dvq=', dvq 1038 call log_mess('algorithm goes out of bound, iteration was left!!!!!!',DEBUG_DEF) 1039 print *, 'gros dvq=', dvq 1040 fm1=INIFINY_REAL_VALUE 1041 fp1=INIFINY_REAL_VALUE 1042 f=INIFINY_REAL_VALUE 1043 return 1044 endif 1045 1046 Hy(kd)=(pdd(current_chr,kd,i,posi_h1)*(dvq))+Hy(kd) 1047 enddo b_pdd 1048 Sy(rangy(kkd))=Sy(rangy(kkd))+(probg(current_chr,ig)*Hy(kd)) 1049 1050 ENDIF 1051 enddo b_indgm 1052 enddo b_gme 1053 enddo b_me 1054 enddo b_pe 1055 ! CALCUL DES DENOMINATEUR POUR CHAQUE KKD 1056 do ip=1,np 1057 fp0(ip)=0.d0 1058 nm1=nmp(ip)+1 1059 nm2=nmp(ip+1) 1060 do jm=nm1,nm2 1061 nd1=ndm(jm)+1 1062 nd2=ndm(jm+1) 1063 fm0(jm)=0.d0 1064 do kkd=nd1,nd2 1065 IF(presentc(icar,kkd)) then 1066 Ssy(kkd)=0.d0 1067 ind_sort : do iid=1,ntot 1068 if (yr(rangy(kkd)).GT.yr(iid)) exit ind_sort 1069 Ssy(kkd)=Ssy(kkd)+Sy(iid) 1070 enddo ind_sort 1071 ENDIF !fin de la condition sur les perf 1072 enddo !fin boucle sur kkd 1073 ! calcul de la vraisemblance 1074 vmere=0.d0 1075 nbdv=0 1076 ngeno1=ngenom(current_chr,jm)+1 1077 ngeno2=ngenom(current_chr,jm+1) 1078 do ig=ngeno1,ngeno2 1079 gy=1.d0 1080 nd1=ngend(current_chr,ig)+1 1081 nd2=ngend(current_chr,ig+1) 1082 vpf=0.d0 1083 do kd=nd1,nd2 1084 kkd=ndesc(current_chr,kd) 1085 if(presentc(icar,kkd).and.ndelta(icar,kkd).EQ.1) then 1086 nbdv=nbdv+1 1087 ! condition pour éviter des erreurs arithmetiques du genre division par 0 huge= plus grande valeur possible et tiny=plus petite valeur possible 1088 ! if (Ssy(kkd).GE.huge(Ssy(kkd))) Ssy(kkd)=huge(Ssy(kkd)) 1089 ! if (Ssy(kkd).LE.tiny(Ssy(kkd))) Ssy(kkd)=tiny(Ssy(kkd)) 1090 gy=gy*(Hy(kd)*ntot)/(Ssy(kkd)) 1091 ! print *,ntot,'sous H1', kd,Hy(kd),Ssy(kkd), 'gy',gy, 'probm', probg(current_chr,ig) 1092 endif 1093 enddo !boucle kkd 1094 vmere= vmere+(probg(current_chr,ig)*gy) 1095 !print *, 'vmere=', vmere 1096 1097 enddo !boucle ig 1098 if (vmere == 0) then 1099 call log_mess('algorithm goes out of bound, iteration was left!!!!!!',DEBUG_DEF) 1100 fm1=INIFINY_REAL_VALUE 1101 fp1=INIFINY_REAL_VALUE 1102 f=INIFINY_REAL_VALUE 1103 return 1104 else 1105 fm1(jm)=-dlog(vmere)+(nbdv*(dlog(ntot*1.d0))) 1106 ! print *, 'fm', fm1(jm), vmere 1107 endif 1108 fp1(ip)=fp1(ip)+fm1(jm) 1109 f=f+fm1(jm) 1110 enddo !fin boucle sur jm 1111 enddo !fin boucle sur np 1112 !print *,'f:',f,'effet QTL', x(1) 1113 return 1114 end subroutine funct_1qtl_modlin_cox
init_analyse_modlin_cox
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
init_analyse_modlin_cox
DESCRIPTION
Initialisation/allocation of solution/buffer arrays
NOTES
SOURCE
209 subroutine init_analyse_modlin_cox 210 integer :: stat 211 allocate (sig1(np),STAT=stat) 212 call check_allocate(stat,'sig1 [m_qtlmap_analyse_modlin]') 213 allocate (xmu1p(np),STAT=stat) 214 call check_allocate(stat,'xmu1p [m_qtlmap_analyse_modlin]') 215 allocate (xmu1m(nm),STAT=stat) 216 call check_allocate(stat,'xmu1m [m_qtlmap_analyse_modlin]') 217 allocate (fp0(np),STAT=stat) 218 call check_allocate(stat,'fp0 [m_qtlmap_analyse_modlin_cox]') 219 ! allocate (fp1(np),STAT=stat) 220 ! call check_allocate(stat,'fp1 [m_qtlmap_analyse_modlin_cox]') 221 allocate (fm0(nm),STAT=stat) 222 call check_allocate(stat,'fm0 [m_qtlmap_analyse_modlin_cox]') 223 ! allocate (fm1(nm),STAT=stat) 224 ! call check_allocate(stat,'fm1 [m_qtlmap_analyse_modlin_cox]') 225 ! allocations used in calcul_rang 226 allocate (yr(nd),STAT=stat) 227 call check_allocate(stat,'yr [m_qtlmap_analyse_modlin_cox]') 228 allocate (yt(nd),STAT=stat) 229 call check_allocate(stat,'yt [m_qtlmap_analyse_modlin_cox]') 230 allocate (rangy(nd),STAT=stat) 231 call check_allocate(stat,'rangy [m_qtlmap_analyse_modlin_cox]') 232 233 ! allocation used in opti_0qtl_modlin_cox 234 ! allocate (Hy(nd),STAT=stat) 235 ! call check_allocate(stat,'Hy [m_qtlmap_analyse_modlin_cox]') 236 ! allocate (Sy(nd),STAT=stat) 237 ! call check_allocate(stat,'Sy [m_qtlmap_analyse_modlin_cox]') 238 ! allocate (Ssy(nd),STAT=stat) 239 ! call check_allocate(stat,'Ssy [m_qtlmap_analyse_modlin_cox]') 240 241 end subroutine init_analyse_modlin_cox
opti_0qtl_modlin_cox
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
opti_0qtl_modlin_cox
DESCRIPTION
Calcul de la vraisemblance 0 QTL, 1 caractere , effets parasites inclus
INPUTS
ic : index of the trait est_moy : general mean
SOURCE
340 subroutine opti_0qtl_modlin_cox(ic,est_moy) 341 use m_qtlmap_analyse_lin_gen, only :contingence_cox,precision,vecsol,precis, & 342 nbnivest,ntniv,par0,precis0,vecsol0,xx 343 344 integer , intent(in) :: ic 345 logical , intent(in) :: est_moy 346 ! 347 ! Divers 348 349 integer :: iuser(1) 350 integer ,dimension(:),allocatable :: iw 351 real (kind=dp) ,dimension(:),allocatable :: par,borni,borns,w 352 real (kind=dp) :: user(1) 353 real (kind=dp) :: f 354 integer :: npar,ibound,ip,ix,ifail,i,indest,ntotp 355 integer :: km,jm 356 logical itest 357 ! 358 !****************************************************************************** 359 ! 360 ! recherche des elements estimables � partir de la d�composition de choleski 361 ! comptage et recodification de ces elements (on commence par mettre � 0 les 362 ! compteurs destin�s aux tests des effes fix�s 363 ! 364 ! allocate (Hy(nd),STAT=stat) 365 ! call check_allocate(stat,'Hy [m_qtlmap_analyse_modlin_cox]') 366 allocate (Sy(nd)) 367 allocate (Ssy(nd)) 368 ! call check_allocate(stat,'Ssy [m_qtlmap_analyse_modlin_cox]') 369 370 itest=.false. 371 !print *, 'est_moy',est_moy 372 call contingence_cox(ic,0,itest,est_moy) 373 call precision(xx,precis) 374 ! Parametres de maximisation 375 !npar=np+nbnivest 376 !on n'a pas � prendre en compte un ecart type par pere avec cox 377 npar=nbnivest 378 379 allocate (borni(npar)) 380 allocate (borns(npar)) 381 allocate (par(npar)) 382 383 ibound=0 384 ! pas besoin d'initialiser les ecart types par pere comme il n' y en a pas 385 !do ip=1,np 386 ! borni(ip)=1.d-6 387 ! borns(ip)=1.d6 388 !end do 389 !do ix=np+1,npar 390 ! on a passé les bornes de 1.d4 à 50 car sinon on risque un plantage avec exp() 391 ! un risque de exp(100)=5.d21 s'est largement suffisant /ref 392 do ix=1,npar 393 borni(ix)=1.do-6 394 borns(ix)=1.d4 395 end do 396 ! 397 ! Point de depart 398 ! sous cox il y a ni moyenne ni ecart type intra pere 399 !par(np+1)=1.d0 400 !do ip=1,np 401 !par(ip)=sig0(ip) 402 !par(np+1)=par(np+1)+xmu0p(ip) 403 !end do 404 !par(np+1)=par(np+1)/dble(np) 405 !do ix=np+2,npar 406 do ix= 1, npar 407 par(ix)=1.d0 408 end do 409 ! 410 ! Optimisation de la vraisemblance 411 ifail=1 412 current_ic = ic ! pour le mode lineaire.... 413 !print *, 'npar sous H0',npar 414 if (npar==0) then 415 call funct_0qtl_modlin_cox(npar,par,f,iuser,user) 416 else 417 call minimizing_funct(npar,ibound,funct_0qtl_modlin_cox,borni,borns,par,f,iuser,user,ifail) 418 endif 419 f0=f 420 ! print *, 'Likelihood sous HO=', f0 421 do i=1,npar 422 par0(i)=par(i) 423 ! print *, 'sous H0: par(i) =', par(i), 'npar=',npar 424 end do 425 do i=1,ntniv 426 vecsol0(i)=vecsol(i) 427 precis0(i)=precis(i) 428 end do 429 430 do ip = 1,np 431 ! sig1(ip)=par(ip) 432 sig1(ip)=0.d0 433 end do 434 435 ! indest=1 436 indest=0 437 if (np.ge.2) then 438 do ip = 1,np 439 ! if (vecsol(1+ip)) then 440 if (vecsol(ip)) then 441 indest=indest+1 442 ! xmu1p(ip)=par(np+indest) 443 xmu1p(ip)=par(indest) 444 end if 445 end do 446 endif 447 448 ! ntot=np+1 449 ntotp=0 450 if (np.ge.2) ntotp=np 451 km=0 452 if (nm.ge.2) then 453 do jm=1,nm 454 if (estime(ic,jm).and.opt_sib.eq.2)then 455 km=km+1 456 if(vecsol(ntotp+km)) then 457 indest=indest+1 458 ! xmu1m(jm)=par(np+indest) 459 xmu1m(jm)=par(indest) 460 end if 461 end if 462 end do 463 endif 464 deallocate (borni) 465 deallocate (borns) 466 deallocate (par) 467 deallocate (Sy) 468 deallocate (Ssy) 469 470 return 471 end subroutine opti_0qtl_modlin_cox
opti_1qtl_modlin_cox
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
opti_1qtl_modlin_cox
DESCRIPTION
Calcul de la vraisemblance 1 QTL, 1 caractere , effets parasites inclus
INPUTS
ic : index of the trait
OUTPUTS
lrtsol : LRT : curves and maximum under H(nqtl) fmax : maximum value of the likelihood supnbnivest : number of parameter estimated at the maximum value
SOURCE
640 subroutine opti_1qtl_modlin_cox(ic,lrtsol,fmax,supnbnivest) 641 642 integer , intent(in) :: ic 643 type(TYPE_LRT_SOLUTION) , intent(out) :: lrtsol 644 real (kind=dp) , intent(out) :: fmax 645 integer , intent(out) :: supnbnivest 646 ! 647 ! 648 ! Divers 649 650 integer :: iuser(1) 651 integer ,dimension(:),allocatable :: iw 652 real (kind=dp) ,dimension(:),allocatable :: val,par,borni,borns,w 653 real (kind=dp) :: user(1) 654 real (kind=dp) :: f 655 integer :: npar,indest 656 integer :: km,jm,chr,ifem 657 real(kind=dp) ,dimension(:,:) ,pointer :: listef 658 integer ,dimension(:,:) ,pointer :: listenbnivest 659 real(kind=dp) ,dimension(:,:,:) ,pointer :: listepar 660 661 real (kind=dp) :: xlrt_t,f1 662 663 logical itest,stvecsol(size(vecsol1)) 664 integer :: ip,i,n,ilong,ibound,j,ifail,ii,iip,indexchr(nchr),ntotal,nprime 665 integer :: save_nteffmax,save_ntnivmax 666 667 call new(1,lrtsol) 668 !****************************************************************************** 669 ! 670 ! Calcul de la vraisemblance sous H1 671 672 ! 673 ! initialisation 674 ! on utilisera : 675 ! PAR (vecteur des param�tres � optimiser), 676 ! CORNIV (vecteur des positions, parmi NTNIV, des effets extimables), 677 ! SOLVEC (vecteur disant l'estimabilit� des effets) 678 ! VAL le vecteur complet (ntniv positions) des valeurs des niveaux des effets 679 ! STVAL,STVECSOL et STCORNIV les copies de VAL, VECSOL et CORNIV � l'it�ration 680 ! n-1 quand on analyse la position n 681 682 ibound=0 683 ! pas besoin d'initialiser les ecart types par pere comme il n' y en a pas 684 !do ip=1,np 685 ! borni(ip)=1.d-6 686 ! borns(ip)=1.d6 687 !end do 688 !do ix=np+1,npar 689 690 do i=1,size(vecsol) 691 vecsol(i)=.true. 692 end do 693 694 lrtsol%lrtmax=-1.d75 695 lrtsol%chrmax=0 696 lrtsol%nxmax=0 697 698 current_ic = ic 699 700 ntotal=0 701 do chr=1,nchr 702 ntotal=ntotal+get_npo(chr) 703 indexchr(chr)=ntotal 704 end do 705 706 allocate (listef(nchr,get_maxnpo())) 707 allocate (listenbnivest(nchr,get_maxnpo())) 708 allocate (listepar(nchr,get_maxnpo(),np + ntnivmax)) 709 710 call end_contingence 711 save_ntnivmax = ntnivmax 712 save_nteffmax = nteffmax 713 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(par,val) & 714 !$OMP PRIVATE(i,ifem,iip,ii,stvecsol,npar,j,ifail,f1,n,chr,borni,borns) 715 ntnivmax = save_ntnivmax 716 nteffmax = save_nteffmax 717 allocate ( val( ntnivmax ) ) 718 allocate ( par( ntnivmax) ) 719 allocate ( borni( ntnivmax) ) 720 allocate ( borns( ntnivmax) ) 721 allocate (fp1(np)) 722 allocate (fm1(nm)) 723 allocate (Hy(nd)) 724 allocate (Sy(nd)) 725 allocate (Ssy(nd)) 726 727 ! 728 ! Point de depart 729 ! 730 do i=1,size(par) 731 par(i)=1.d0 732 end do 733 734 do i=1,size(val) 735 val(i)=par(i) 736 end do 737 738 739 call init_contingence 740 ! 741 ! Marche le long du chromosome 742 !$OMP DO 743 do nprime=1,ntotal 744 n=nprime 745 do chr=nchr,1,-1 746 if ( indexchr(chr) >= nprime) then 747 current_chr = chr 748 else 749 n=n-get_npo(chr) 750 end if 751 end do 752 chr=current_chr 753 ! 754 ! on stocke les conditions en n 755 do i=1,ntniv 756 stvecsol(i)=vecsol(i) 757 end do 758 ! 759 ! preparation de la matrice d'incidence 760 ! 761 call prepinc(current_chr,n,ic,"LA ") 762 ! 763 ! recherche des elements estimables � partir de la d�composition de choleski 764 ! comptage et recodification de ces elements (on commence par mettre � 0 les 765 ! compteurs destin�s aux tests des effes fix�s 766 ! 767 768 itest=.false. 769 call contingence_cox(ic,1,itest,.false.) 770 ! Parametres de maximisation 771 npar=nbnivest 772 ! on a passé les bornes de 1.d4 à 50 car sinon on risque un plantage avec exp() 773 ! un risque de exp(100)=5.d21 s'est largement suffisant /ref 774 do i=1,npar 775 borni(i)=1.do-5 776 borns(i)=1.d4 777 end do 778 779 780 781 ! 782 ! Point de depart 783 ! sous cox il y a ni moyenne ni ecart type intra pere 784 785 do i= 1, npar 786 par(i)=1.d0 787 end do 788 789 do i=1,ntniv 790 if(vecsol(i))then 791 par(i)=1.d0 792 if(stvecsol(i)) par(i)=val(i) 793 end if 794 stvecsol(i)=vecsol(i) 795 end do 796 ! 797 ! Optimisation de la vraisemblance 798 ifail=1 799 posi_h1=n 800 call minimizing_funct(npar,ibound,funct_1qtl_modlin_cox,borni,borns,par,f1,iuser,user,ifail) 801 print *,n,f1 802 803 do i=1,ntniv 804 if(vecsol(i)) then 805 val(i)=par(i) 806 else 807 val(i)=9999.d0 808 end if 809 end do 810 ! 811 ! on garde les valeurs du LRT pour pouvopir dessiner la courbe de vraisemblance 812 ! 813 if ( f1 < INIFINY_REAL_VALUE ) then 814 lrtsol%lrt1(chr,n)=-2.d0*(f1-f0) 815 else 816 lrtsol%lrt1(chr,n)=0.d0 817 end if 818 819 listef(chr,n)=f1 820 listenbnivest(chr,n)=nbnivest 821 listepar(chr,n,:npar)=par(:npar) 822 ! 823 ! on met les profil et effets QTL / prog�niteur dnas les ficheir ad hoc 824 ! 825 iip=0 826 do ii=1,np 827 lrtsol%xlrp(chr,ii,n)=-2.d0*(fp1(ii)-fp0(ii)) 828 if ( vecsol(ii)) then 829 iip=iip+1 830 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 831 lrtsol%pater_eff(chr,ii,n)=par(iip) 832 833 end if 834 end do 835 836 ifem=0 837 do ii=1,nm 838 lrtsol%xlrm(chr,ii,n)=-2.d0*(fm1(ii)-fm0(ii)) 839 if ( estime(ic,ii)) then 840 ifem=ifem+1 841 if ( vecsol(np+ifem) ) then 842 iip=iip+1 843 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 844 lrtsol%mater_eff(chr,ifem,n)=par(iip) 845 end if 846 end if 847 end do 848 849 end do 850 !$OMP END DO 851 call end_contingence 852 deallocate ( val ) 853 deallocate ( par ) 854 deallocate (fp1) 855 deallocate (fm1) 856 deallocate (Hy) 857 deallocate (Sy) 858 deallocate (Ssy) 859 deallocate ( borni ) 860 deallocate ( borns ) 861 !$OMP END PARALLEL 862 863 !recherche du maximum 864 do chr=1,nchr 865 do n=1,get_npo(chr) 866 if(lrtsol%lrtmax(0) < lrtsol%lrt1(chr,n)) then 867 lrtsol%nxmax(0)=n 868 lrtsol%chrmax(0)=chr 869 lrtsol%lrtmax(0)=lrtsol%lrt1(chr,n) 870 fmax=listef(chr,n) 871 supnbnivest=listenbnivest(chr,n) 872 par1(:supnbnivest+np)=listepar(chr,n,:supnbnivest+np) 873 end if 874 end do 875 end do 876 deallocate (listef) 877 deallocate (listenbnivest) 878 deallocate (listepar) 879 880 ! 881 ! on calcule la pr�cision des estimation au point correspondant 882 ! au LRT maximum 883 ! 884 call init_contingence 885 call prepinc(lrtsol%chrmax(0),lrtsol%nxmax(0),ic,"LA ") 886 call contingence_cox(ic,1,itest,.false.) 887 call precision(xx,precis) 888 889 do i=1,ntniv 890 precis1(i)=precis(i) 891 vecsol1(i)=vecsol(i) 892 end do 893 894 return 895 end subroutine opti_1qtl_modlin_cox
set_solution_hypothesis0
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
set_solution_hypothesis0
DESCRIPTION
INPUTS
OUTPUTS
SOURCE
1300 subroutine set_solution_hypothesis0(ic,incsol) 1301 integer ,intent(in) :: ic 1302 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1303 1304 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 1305 integer :: indest,isol,ipar,nlevel,ief,i,ife 1306 1307 real(kind=dp) :: par_ref 1308 1309 incsol%hypothesis=0 1310 ! Polygenic family 1311 nteff=0 1312 if (np.ge.2) nteff = 1 1313 if ( count(estime(ic,:)) > 1 ) nteff = 2 1314 1315 !Fixed effect and covariate 1316 if ( modele(ic,1) > 0 ) nteff = nteff+1 1317 if ( modele(ic,2) > 0 ) nteff = nteff+1 1318 1319 maxNbPar = max(np,count(estime(ic,:))) 1320 1321 !max numbre de covariable ? 1322 maxNbPar = max(maxNbPar,modele(ic,2)) 1323 1324 nlevel=0 1325 !max nombre de niveau pour un effet fixe ? 1326 do i=1,modele(ic,1) 1327 nlevel=nlev(modele(ic,3+i))+nlevel 1328 end do 1329 maxNbPar = max(maxNbPar,nlevel) 1330 1331 allocate (incsol%groupeName(nteff)) 1332 allocate (incsol%eqtl_print(nteff)) 1333 allocate (incsol%nbParameterGroup(nteff)) 1334 allocate (incsol%parameterName(nteff,maxNbPar)) 1335 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1336 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1337 1338 incsol%parameterVecsol=.true. 1339 incsol%eqtl_print=.true. 1340 1341 ieff=0 1342 indest = 0 1343 isol=0 1344 1345 if (np.ge.2) then 1346 ieff=ieff+1 1347 par_ref=par0(indest+1) 1348 incsol%groupeName(ieff) = 'Sire polygenic effects' 1349 incsol%nbParameterGroup(ieff)=np 1350 do ip=1,np 1351 isol=isol+1 1352 if ( vecsol0(isol) ) then 1353 indest=indest+1 1354 incsol%paramaterValue(ieff,ip) = (par0(indest)/par_ref) 1355 else 1356 incsol%paramaterValue(ieff,ip) = 1.d0 1357 end if 1358 1359 incsol%parameterName(ieff,ip) ='Sire '//trim(pere(ip)) 1360 incsol%parameterVecsol(ieff,ip) = vecsol0(isol) 1361 end do 1362 endif 1363 if ( count(estime(ic,:)) > 1 ) then 1364 par_ref=par0(indest+1) 1365 ieff = ieff +1 1366 incsol%groupeName(ieff) = 'Dam polygenic effects' 1367 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 1368 ifem=0 1369 do ip=1,np 1370 do jm=nmp(ip)+1,nmp(ip+1) 1371 if ( estime(ic,jm)) then 1372 isol=isol+1 1373 ifem=ifem+1 1374 if ( vecsol0(isol) ) then 1375 indest=indest+1 1376 incsol%paramaterValue(ieff,ifem) = (par0(indest)/par_ref) 1377 else 1378 incsol%paramaterValue(ieff,ifem) = 1.d0 1379 end if 1380 1381 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1382 incsol%parameterVecsol(ieff,ifem) = vecsol0(isol) 1383 1384 end if 1385 end do 1386 end do 1387 end if 1388 1389 !Fixed effect 1390 if ( modele(ic,1) > 0 ) then 1391 ieff=ieff+1 1392 incsol%eqtl_print(ieff)=.false. 1393 incsol%groupeName(ieff) = 'fixed effects' 1394 incsol%nbParameterGroup(ieff)=nlevel 1395 ife=0 1396 do ief=1,nbef 1397 par_ref=par0(indest+1) 1398 do i=1,nlev(modele(ic,3+ief)) 1399 ife=ife+1 1400 isol=isol+1 1401 1402 if ( vecsol0(isol) ) then 1403 indest=indest+1 1404 incsol%paramaterValue(ieff,ife) = (par0(indest)/par_ref) 1405 else 1406 incsol%paramaterValue(ieff,ife) = 1.d0 1407 end if 1408 incsol%parameterName(ieff,ife) = trim(namefix(modele(ic,3+ief)))//' level '//trim(str(i)) 1409 incsol%parameterVecsol(ieff,ife) = vecsol0(isol) 1410 1411 end do 1412 end do 1413 end if 1414 1415 !Covariate 1416 if ( modele(ic,2) > 0 ) then 1417 ieff=ieff+1 1418 incsol%eqtl_print(ieff)=.false. 1419 incsol%groupeName(ieff) = 'Covariates' 1420 incsol%nbParameterGroup(ieff)=modele(ic,2) 1421 1422 do ief=1,modele(ic,2) 1423 isol=isol+1 1424 if ( vecsol0(isol) ) then 1425 indest=indest+1 1426 incsol%paramaterValue(ieff,ief) = (par0(indest)) 1427 else 1428 incsol%paramaterValue(ieff,ief) = 1.d0 1429 end if 1430 incsol%parameterName(ieff,ief) = trim(namecov(modele(ic,3+modele(ic,1)+ief))) 1431 incsol%parameterVecsol(ieff,ief) = vecsol0(isol) 1432 end do 1433 end if 1434 1435 end subroutine set_solution_hypothesis0
set_solution_hypothesis1
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
set_solution_hypothesis1
DESCRIPTION
INPUTS
OUTPUTS
SOURCE
1450 subroutine set_solution_hypothesis1(ic,incsol) 1451 integer ,intent(in) :: ic 1452 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1453 1454 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 1455 integer :: ntlev,nbtp,jef,lp,indest 1456 integer :: isol,ipar,nlevel,ief,i,ife 1457 real(kind=dp) :: par_ref 1458 1459 incsol%hypothesis=1 1460 ! Polygenic family, QTL effect 1461 nteff=1 1462 if (np.ge.2) nteff = 2 1463 if ((np==1).and. count(estime(ic,:)) ==1 ) nteff = 2 1464 if ((np==1).and. count(estime(ic,:)) > 1 ) nteff = 3 1465 if ((np.ge.2).and. count(estime(ic,:))==1 ) nteff = 3 1466 if ((np.ge.2).and. count(estime(ic,:))> 1 ) nteff = 4 1467 !Fixed effect and covariate 1468 if ( modele(ic,1) > 0 ) nteff = nteff+1 1469 if ( modele(ic,2) > 0 ) nteff = nteff+1 1470 1471 maxNbPar = max(np,count(estime(ic,:))) 1472 !max numbre de covariable ? 1473 maxNbPar = max(maxNbPar,modele(ic,2)) 1474 1475 nlevel=0 1476 !max nombre de niveau pour un effet fixe ? 1477 do i=1,modele(ic,1) 1478 nlevel=nlev(modele(ic,3+i))+nlevel 1479 end do 1480 maxNbPar = max(maxNbPar,nlevel) 1481 1482 ntlev=1 1483 nbtp = 3 + modele(ic,1)+modele(ic,2)!+modele(ic,3) 1484 do jef=1,modele(ic,3) 1485 ntlev=ntlev*nlev(modele(ic,nbtp+jef)) 1486 end do 1487 1488 !max nombre de niveau pour un effet fixe en interaction avec le qtl ? 1489 maxNbPar = max(maxNbPar,ntlev*np) 1490 maxNbPar = max(maxNbPar,ntlev*count(estime(ic,:))) 1491 1492 allocate (incsol%groupeName(nteff)) 1493 allocate (incsol%eqtl_print(nteff)) 1494 allocate (incsol%nbParameterGroup(nteff)) 1495 allocate (incsol%parameterName(nteff,maxNbPar)) 1496 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1497 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1498 allocate (incsol%qtl_groupeName(1,1)) 1499 1500 incsol%parameterVecsol=.true. 1501 incsol%eqtl_print=.true. 1502 1503 ieff=1 1504 incsol%qtl_groupeName(1,1)=ieff 1505 incsol%groupeName(ieff) = 'Sire QTL effects' 1506 incsol%nbParameterGroup(ieff)=np*ntlevp 1507 1508 indest = 0 1509 isol=0 1510 ife=0 1511 do ip=1,np 1512 do lp=1,ntlev 1513 ife=ife+1 1514 isol=isol+1 1515 ! print *, 'vecsol',vecsol1(1), 'par', par1(1) 1516 if ( vecsol1(isol) ) then 1517 indest=indest+1 1518 incsol%paramaterValue(ieff,ife) = (par1(indest)) 1519 else 1520 incsol%paramaterValue(ieff,ife) = 1.d0 1521 end if 1522 1523 incsol%parameterName(ieff,ife) ='Sire '//trim(pere(ip)) 1524 incsol%parameterVecsol(ieff,ife) = vecsol1(isol) 1525 end do 1526 end do 1527 1528 if ( count(estime(ic,:)) > 0 ) then 1529 ieff = ieff +1 1530 incsol%groupeName(ieff) = 'Dam QTL effects' 1531 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 1532 ifem=0 1533 do jm=1,nm 1534 if ( estime(ic,jm)) then 1535 do lp=1,ntlev 1536 isol=isol+1 1537 ifem=ifem+1 1538 if ( vecsol1(isol) ) then 1539 indest=indest+1 1540 incsol%paramaterValue(ieff,ifem) = (par1(indest)) 1541 else 1542 incsol%paramaterValue(ieff,ifem) = 1.d0 1543 end if 1544 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm)) 1545 incsol%parameterVecsol(ieff,ifem) = vecsol1(isol) 1546 end do 1547 end if 1548 end do 1549 end if 1550 1551 1552 1553 if (np.ge.2) then 1554 ieff=ieff+1 1555 incsol%groupeName(ieff) = 'Sire polygenic effects' 1556 incsol%nbParameterGroup(ieff)=np 1557 par_ref=par1(indest+1) 1558 do ip=1,np 1559 isol=isol+1 1560 if ( vecsol1(isol) ) then 1561 indest=indest+1 1562 incsol%paramaterValue(ieff,ip) = (par1(indest)/par_ref) 1563 else 1564 incsol%paramaterValue(ieff,ip) = 1.d0 1565 end if 1566 1567 incsol%parameterName(ieff,ip) ='Sire '//trim(pere(ip)) 1568 incsol%parameterVecsol(ieff,ip) = vecsol1(isol) 1569 end do 1570 endif 1571 1572 if ( count(estime(ic,:)) > 1 ) then 1573 ieff = ieff +1 1574 incsol%groupeName(ieff) = 'Dam polygenic effects' 1575 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 1576 ifem=0 1577 par_ref=par1(indest+1) 1578 do ip=1,np 1579 do jm=nmp(ip)+1,nmp(ip+1) 1580 if ( estime(ic,jm)) then 1581 isol=isol+1 1582 ifem=ifem+1 1583 if ( vecsol1(isol) ) then 1584 indest=indest+1 1585 incsol%paramaterValue(ieff,ifem) = (par1(indest)/par_ref) 1586 else 1587 incsol%paramaterValue(ieff,ifem) = 1.d0 1588 end if 1589 1590 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1591 incsol%parameterVecsol(ieff,ifem) = vecsol1(isol) 1592 end if 1593 end do 1594 end do 1595 end if 1596 1597 !Fixed effect 1598 if ( modele(ic,1) > 0 ) then 1599 ieff=ieff+1 1600 incsol%eqtl_print(ieff)=.false. 1601 incsol%groupeName(ieff) = 'fixed effects' 1602 incsol%nbParameterGroup(ieff)=nlevel 1603 ife=0 1604 do ief=1,nbef 1605 par_ref=par1(indest+1) 1606 do i=1,nlev(modele(ic,3+ief)) 1607 isol=isol+1 1608 ife=ife+1 1609 if ( vecsol1(isol) ) then 1610 indest=indest+1 1611 incsol%paramaterValue(ieff,ife) = (par1(indest)/par_ref) 1612 else 1613 incsol%paramaterValue(ieff,ife) = 1.d0 1614 end if 1615 incsol%parameterName(ieff,ife) = trim(namefix(modele(ic,3+ief)))//' level '//trim(str(i)) 1616 incsol%parameterVecsol(ieff,ife) = vecsol1(isol) 1617 end do 1618 end do 1619 end if 1620 1621 !Covariate 1622 if ( modele(ic,2) > 0 ) then 1623 ieff=ieff+1 1624 incsol%eqtl_print(ieff)=.false. 1625 incsol%groupeName(ieff) = 'Covariates' 1626 incsol%nbParameterGroup(ieff)=modele(ic,2) 1627 1628 do ief=1,modele(ic,2) 1629 isol=isol+1 1630 if ( vecsol1(isol) ) then 1631 indest=indest+1 1632 incsol%paramaterValue(ieff,ief) =(par1(indest)) 1633 else 1634 incsol%paramaterValue(ieff,ief) = 1.d0 1635 end if 1636 incsol%parameterName(ieff,ief) = trim(namecov(modele(ic,3+modele(ic,1)+ief))) 1637 incsol%parameterVecsol(ieff,ief) = vecsol1(isol) 1638 end do 1639 end if 1640 end subroutine set_solution_hypothesis1
test_lin_cox
[ Top ] [ m_qtlmap_analyse_modlin_cox ] [ Subroutines ]
NAME
test_lin_cox
DESCRIPTION
Test des differents effets de nuisance du mod�le par une LRT comp�r� � une chi2
INPUTS
OUTPUTS
SOURCE
1129 subroutine test_lin_cox(chr,ic,est_moy,supnbnivest,fmax,nposx) 1130 use m_qtlmap_analyse_lin_gen, only : prepinc,contingence_cox,nbef,nbco,meff,mcov,mint,nbnivest 1131 1132 logical , intent(in) :: est_moy 1133 integer , intent(in) :: chr,ic 1134 integer , intent(in) :: supnbnivest 1135 real (kind=dp) , intent(in) :: fmax 1136 integer , intent(in) :: nposx 1137 1138 ! 1139 ! Divers 1140 logical :: itest 1141 integer :: iuser(1) 1142 real (kind=dp) ,dimension(:),allocatable :: par,borni,borns 1143 1144 real (kind=dp) :: user(1),prob,xlrt_t,f1 1145 integer :: nbint,iecd,iecq,ief,ibound,npar,ip,i,ifail 1146 integer :: nbreduit 1147 1148 write(nficout,600) 1149 600 format(//,80('*')/'test of the effets of the model',//) 1150 1151 ! 1152 ! preparation de la matrice d'incidence 1153 ! 1154 call prepinc(chr,nposx,ic,"LA ") 1155 ! 1156 ! on recalcule la vraisemblance en retirant les effets du mod�les un � un 1157 ! 1158 ! 1159 ! on commence par les effets hierarchis�s dans les effets qtl 1160 ! 1161 1162 nbef=modele(ic,1) 1163 nbco=modele(ic,2) 1164 nbint=modele(ic,3) 1165 1166 iecd=0 1167 iecq=0 1168 1169 allocate ( par( np + ntnivmax ) ) 1170 allocate ( borni( np + ntnivmax ) ) 1171 allocate ( borns( np + ntnivmax ) ) 1172 1173 allocate (fp1(np)) 1174 allocate (fm1(nm)) 1175 allocate (Hy(nd)) 1176 allocate (Sy(nd)) 1177 allocate (Ssy(nd)) 1178 1179 do ief=1,nbint+nbef+nbco 1180 meff=0 1181 mcov=0 1182 mint=0 1183 if(ief.le.nbef)meff=ief 1184 if(ief.gt.nbef.and.ief.le.(nbef+nbco))mcov=ief-nbef 1185 if(ief.gt.(nbef+nbco))mint=ief-(nbef+nbco) 1186 ! 1187 ! recherche des elements estimables � partir de la d�composition de choleski 1188 ! comptage et recodification de ces elements 1189 ! 1190 1191 itest=.true. 1192 call contingence_cox(ic,1,itest,est_moy) 1193 1194 ! 1195 ! la r�dustion du nombre d'effet estim�e est calcul�e 1196 ! 1197 nbreduit=supnbnivest-nbnivest 1198 1199 ! 1200 ! Parametres de maximisation 1201 ibound=0 1202 npar=nbnivest 1203 print *,'npar max', npar 1204 do i=1,npar 1205 borni(i)=1.do-6 1206 borns(i)=1.d20 1207 par(i)=1.d0 1208 end do 1209 ! 1210 ! Point de depart (on reprend les point d'arriv�e pr�c�dents 1211 ! 1212 do i=1,npar 1213 par(i)=1.d0 1214 end do 1215 ! 1216 ! Optimisation de la vraisemblance a la position dx 1217 ifail=1 1218 1219 call minimizing_funct(npar,ibound,funct_1qtl_modlin_cox,borni,borns,par,f1,iuser,user,ifail) 1220 ! 1221 xlrt_t=-2.d0*(fmax-f1) 1222 if(xlrt_t.le.1.do-8)xlrt_t=0.d0 1223 ! 1224 ! impression des tests des effets du mod�le 1225 ! 1226 1227 nbef=modele(ic,1) 1228 nbco=modele(ic,2) 1229 nbint=modele(ic,3) 1230 1231 if(ief.le.(nbef+nbco).and.iecd.eq.0) then 1232 iecd=1 1233 write(nficout,610) 1234 610 format(' Direct effects'/) 1235 1236 write(nficout,601) 1237 601 format('Tested effect df. Likelihood p-value'/ & 1238 ' ratio '/) 1239 end if 1240 1241 if(ief.le.nbef) then 1242 ifail=0 1243 ! prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nlev(modele(ic,3+ief))),ifail) 1244 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nbreduit),ifail) 1245 1246 write(nficout,602) trim(namefix(modele(ic,3+ief))),nbreduit,xlrt_t,prob 1247 ! & nlev(modele(ic,3+ief)),xlrt,prob 1248 end if 1249 602 format(a15,3x,i2,6x,f8.2,7x,f5.3) 1250 1251 if(ief.gt.nbef.and.ief.le.(nbef+nbco)) then 1252 ifail=0 1253 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(1.0),ifail) 1254 1255 write(nficout,603) trim(namecov(modele(ic,3+ief))),xlrt_t,prob 1256 end if 1257 603 format(a15,4x,'1',6x,f8.2,7x,f5.3) 1258 1259 if(ief.gt.(nbef+nbco).and.iecq.eq.0)then 1260 iecq=1 1261 write(nficout,611) 1262 611 format(/' Intra qtl effects'/) 1263 write(nficout,601) 1264 end if 1265 1266 if(ief.gt.(nbef+nbco)) then 1267 ifail=0 1268 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nbreduit),ifail) 1269 write(nficout,602) trim(namefix(modele(ic,3+ief))),nbreduit,xlrt_t,prob 1270 1271 end if 1272 1273 end do 1274 1275 deallocate (fp1) 1276 deallocate (fm1) 1277 deallocate (Hy) 1278 deallocate (Sy) 1279 deallocate (Ssy) 1280 deallocate ( par ) 1281 deallocate ( borni ) 1282 deallocate ( borns ) 1283 1284 return 1285 end subroutine test_lin_cox