m_qtlmap_analyse_modlin
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_analyse_modlin
DESCRIPTION
NOTES
SEE ALSO
current_chr
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
current_chr
DESCRIPTION
the current chromosome to analyse while the likelihood calculus
current_ic
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
current_ic
DESCRIPTION
the current trait to analyse while the likelihood calculus
f0
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
f0
DESCRIPTION
value of the likelihood under H0
fm0
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
fm0
DESCRIPTION
value of the likelihood by full-sib family under H0
DIMENSIONS
nm
fm1
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
fm1
DESCRIPTION
value of the likelihood by full-sib family under H1
DIMENSIONS
nm
fp0
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
fp0
DESCRIPTION
value of the likelihood by sire family under H0
DIMENSIONS
np
fp1
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
fp1
DESCRIPTION
value of the likelihood by sire family under H1
DIMENSIONS
np
sig1
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
sig1
DESCRIPTION
The standart deviation under H0
DIMENSIONS
np
xmu1g
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
xmu1g
DESCRIPTION
The general mean
xmu1m
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
xmu1m
DESCRIPTION
The polygenic mean for each full sib family
DIMENSIONS
nm
xmu1p
[ Top ] [ m_qtlmap_analyse_modlin ] [ Variables ]
NAME
xmu1p
DESCRIPTION
The polygenic mean for each sire family under H0
DIMENSIONS
np
end_analyse_modlin
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
end_analyse_modlin
DESCRIPTION
deallocation of solution/buffer arrays
NOTES
SOURCE
182 subroutine end_analyse_modlin 183 deallocate (sig1) 184 deallocate (xmu1p) 185 deallocate (xmu1m) 186 deallocate (fp0) 187 deallocate (fp1) 188 deallocate (fm0) 189 deallocate (fm1) 190 end subroutine end_analyse_modlin
funct_0qtl_modlin
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
funct_0qtl_modlin
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H0
NOTES
SOURCE
322 subroutine funct_0qtl_modlin(n,x,f,iuser,user) 323 use m_qtlmap_analyse_lin_gen, only : nbfem,nbef,ntnifix,nbco,covdir,nivdir,corniv,vecsol 324 325 implicit none 326 integer , intent(in) :: n 327 ! 328 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 329 ! 330 real (kind=dp) ,dimension(n), intent(in) :: x 331 real (kind=dp) , intent(inout) :: f 332 integer , dimension(1), intent(in) :: iuser 333 real (kind=dp) ,dimension(1), intent(in) :: user 334 335 ! Tableaux dimensionnes selon nm, le nombre de meres 336 real (kind=dp) ,dimension(nm) :: effm 337 338 integer ::jnit,ip,nm1,nm2,jm,nd1,nd2,kkd,ilev,ief,ico,icar 339 real (kind=dp) :: sig,var,vmere,vpf,v 340 ! 341 342 !****************************************************************************** 343 ! 344 jnit=3 345 icar = current_ic 346 if (nbfem.eq.0)jnit=2 347 f=0.d0 348 do ip=1,np 349 sig=x(ip) 350 var=sig*sig 351 fp0(ip)=0.d0 352 nm1=nmp(ip)+1 353 nm2=nmp(ip+1) 354 do jm=nm1,nm2 355 effm(jm)=0.d0 356 vmere=0.d0 357 vpf=0.d0 358 ! 359 ! on ne consid�re que les m�res 360 ! 361 nd1=ndm(jm)+1 362 nd2=ndm(jm+1) 363 ! 364 do kkd=nd1,nd2 365 ! 366 if(presentc(icar,kkd)) then 367 effm(jm)=effm(jm)+1.d0 368 v=y(icar,kkd) 369 ! 370 if(vecsol(nivdir(kkd,1))) then 371 ilev=corniv(nivdir(kkd,1)) 372 v=v-x(np+ilev) 373 end if 374 375 ! 376 if(vecsol(nivdir(kkd,2))) then 377 ilev=corniv(nivdir(kkd,2)) 378 v=v-x(np+ilev) 379 end if 380 ! 381 if(estime(icar,jm)) then 382 if(vecsol(nivdir(kkd,3))) then 383 ilev=corniv(nivdir(kkd,3)) 384 v=v-x(np+ilev) 385 end if 386 end if 387 388 do ief=jnit+1,jnit+nbef 389 if(vecsol(nivdir(kkd,ief))) then 390 ilev=corniv(nivdir(kkd,ief)) 391 v=v-x(np+ilev) 392 end if 393 end do 394 395 do ico=1,nbco 396 if(vecsol(ntnifix+ico)) then 397 ilev=corniv(ntnifix+ico) 398 v=v-covdir(kkd,ico)*x(np+ilev) 399 end if 400 end do 401 402 vpf=vpf+v*v*cd(icar,kkd) 403 end if 404 405 end do 406 vmere=dexp(-0.5d0*vpf/var) 407 if (vmere == 0) then 408 fm0(jm)=INIFINY_REAL_VALUE 409 else 410 fm0(jm)=-dlog(vmere)+effm(jm)*dlog(sig) 411 end if 412 fp0(ip)=fp0(ip)+fm0(jm) 413 f=f+fm0(jm) 414 end do 415 416 end do 417 return 418 end subroutine funct_0qtl_modlin
funct_0qtl_modlin_family
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
funct_0qtl_modlin_family
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous H0
NOTES
SOURCE
429 subroutine funct_0qtl_modlin_family(ip,jm,n,x,f,iuser,user) 430 use m_qtlmap_analyse_lin_gen, only : nbfem,nbef,ntnifix,nbco,covdir,nivdir,corniv,vecsol 431 432 implicit none 433 integer , intent(in) :: ip,jm,n 434 ! 435 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 436 ! 437 real (kind=dp) ,dimension(n), intent(in) :: x 438 real (kind=dp) , intent(inout) :: f 439 integer , dimension(1), intent(in) :: iuser 440 real (kind=dp) ,dimension(1), intent(in) :: user 441 442 ! Tableaux dimensionnes selon nm, le nombre de meres 443 real (kind=dp) :: effm 444 445 integer ::jnit,nm1,nm2,nd1,nd2,kkd,ilev,ief,ico,icar 446 real (kind=dp) :: sig,var,vmere,vpf,v 447 ! 448 449 !****************************************************************************** 450 ! 451 jnit=3 452 icar = current_ic 453 if (nbfem.eq.0)jnit=2 454 f=0.d0 455 sig=x(ip) 456 var=sig*sig 457 effm=0.d0 458 vmere=0.d0 459 vpf=0.d0 460 ! 461 ! on ne consid�re que les m�res 462 ! 463 nd1=ndm(jm)+1 464 nd2=ndm(jm+1) 465 ! 466 do kkd=nd1,nd2 467 ! 468 if(presentc(icar,kkd)) then 469 effm=effm+1.d0 470 v=y(icar,kkd) 471 ! 472 if(vecsol(nivdir(kkd,1))) then 473 ilev=corniv(nivdir(kkd,1)) 474 v=v-x(np+ilev) 475 end if 476 477 ! 478 if(vecsol(nivdir(kkd,2))) then 479 ilev=corniv(nivdir(kkd,2)) 480 v=v-x(np+ilev) 481 end if 482 ! 483 if(estime(icar,jm)) then 484 if(vecsol(nivdir(kkd,3))) then 485 ilev=corniv(nivdir(kkd,3)) 486 v=v-x(np+ilev) 487 end if 488 end if 489 490 do ief=jnit+1,jnit+nbef 491 if(vecsol(nivdir(kkd,ief))) then 492 ilev=corniv(nivdir(kkd,ief)) 493 v=v-x(np+ilev) 494 end if 495 end do 496 497 do ico=1,nbco 498 if(vecsol(ntnifix+ico)) then 499 ilev=corniv(ntnifix+ico) 500 v=v-covdir(kkd,ico)*x(np+ilev) 501 end if 502 end do 503 504 vpf=vpf+v*v*cd(icar,kkd) 505 end if 506 507 end do 508 vmere=dexp(-0.5d0*vpf/var) 509 if (vmere == 0) then 510 f=INIFINY_REAL_VALUE 511 else 512 f=-dlog(vmere)+effm*dlog(sig) 513 end if 514 515 516 end subroutine funct_0qtl_modlin_family
funct_1qtl_modlin
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
funct_1qtl_modlin
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous Hypothese 1QTL 1carac
NOTES
SOURCE
789 subroutine funct_1qtl_modlin(n,x,f,iuser,user) 790 791 integer , intent(in) :: n 792 ! 793 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 794 ! 795 real (kind=dp) ,dimension(n), intent(in) :: x 796 real (kind=dp) , intent(inout) :: f 797 integer , dimension(1), intent(in) :: iuser 798 real (kind=dp) ,dimension(1), intent(in) :: user 799 800 ! Tableaux dimensionnes selon nm, le nombre de meres 801 real (kind=dp) ,dimension(nm) :: effm 802 803 ! 804 ! Divers 805 integer :: init,jnit,ip,nm1,nm2,jm,ngeno1,ngeno2,ig 806 integer :: nd1,nd2,kd,kkd,ilev,ief,ico,ic 807 real (kind=dp) :: sig,var,vmere,vpf,v 808 !****************************************************************************** 809 ! 810 811 ic = current_ic 812 init=4 813 jnit=5 814 if (nbfem.eq.0)then 815 init=3 816 jnit=3 817 end if 818 819 f=0.d0 820 do ip=1,np 821 sig=x(ip) 822 var=sig*sig 823 fp1(ip)=0.d0 824 nm1=nmp(ip)+1 825 nm2=nmp(ip+1) 826 do jm=nm1,nm2 827 vmere=0.d0 828 ngeno1=ngenom(current_chr,jm)+1 829 ngeno2=ngenom(current_chr,jm+1) 830 do ig=ngeno1,ngeno2 831 nd1=ngend(current_chr,ig)+1 832 nd2=ngend(current_chr,ig+1) 833 vpf=0.d0 834 effm(jm)=0.d0 835 do kd=nd1,nd2 836 kkd=ndesc(current_chr,kd) 837 838 if(presentc(ic,kkd)) then 839 effm(jm)=effm(jm)+1 840 v=y(ic,kkd) 841 842 ilev=corniv(nivdir(kkd,1)) 843 v=v-x(np+ilev) 844 if(vecsol(nivdir(kkd,2))) then 845 ilev=corniv(nivdir(kkd,2)) 846 v=v-x(np+ilev)*ppt(kd) 847 end if 848 849 if(estime(ic,jm)) then 850 if(vecsol(nivdir(kkd,3))) then 851 ilev=corniv(nivdir(kkd,3)) 852 v=v-x(np+ilev)*pmt(kd) 853 end if 854 end if 855 856 if(vecsol(nivdir(kkd,init))) then 857 ilev=corniv(nivdir(kkd,init)) 858 v=v-x(np+ilev) 859 end if 860 861 if(estime(ic,jm)) then 862 if(vecsol(nivdir(kkd,5))) then 863 ilev=corniv(nivdir(kkd,5)) 864 v=v-x(np+ilev) 865 end if 866 end if 867 868 do ief=jnit+1,jnit+nbef 869 870 if(vecsol(nivdir(kkd,ief))) then 871 ilev=corniv(nivdir(kkd,ief)) 872 v=v-x(np+ilev) 873 end if 874 end do 875 876 do ico=1,nbco 877 if(vecsol(ntnifix+ico)) then 878 ilev=corniv(ntnifix+ico) 879 v=v-covdir(kkd,ico)*x(np+ilev) 880 end if 881 end do 882 vpf=vpf+v*v*cd(ic,kkd) 883 end if 884 end do 885 vmere=vmere+probg(current_chr,ig)*dexp(-0.5d0*vpf/var) 886 end do 887 888 if (vmere == 0) then 889 fm1(jm)=INIFINY_REAL_VALUE 890 f=INIFINY_REAL_VALUE 891 return 892 else 893 fm1(jm)=-dlog(vmere)+dble(effm(jm))*dlog(sig) 894 end if 895 fp1(ip)=fp1(ip)+fm1(jm) 896 f=f+fm1(jm) 897 898 end do 899 end do 900 901 return 902 end subroutine funct_1qtl_modlin
funct_1qtl_modlin_family
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
funct_1qtl_modlin_family
DESCRIPTION
Calcul de la vraisemblance et de ses derivees partielles sous Hypothese 1QTL 1carac
NOTES
SOURCE
913 subroutine funct_1qtl_modlin_family(ip,jm,n,x,f,iuser,user) 914 915 integer , intent(in) :: ip,jm,n 916 ! 917 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 918 ! 919 real (kind=dp) ,dimension(n), intent(in) :: x 920 real (kind=dp) , intent(inout) :: f 921 integer , dimension(1), intent(in) :: iuser 922 real (kind=dp) ,dimension(1), intent(in) :: user 923 924 ! Tableaux dimensionnes selon nm, le nombre de meres 925 real (kind=dp) :: effm 926 927 ! 928 ! Divers 929 integer :: init,jnit,nm1,nm2,ngeno1,ngeno2,ig 930 integer :: nd1,nd2,kd,kkd,ilev,ief,ico,ic 931 real (kind=dp) :: sig,var,vmere,vpf,v 932 !****************************************************************************** 933 ! 934 935 ic = current_ic 936 init=4 937 jnit=5 938 if (nbfem.eq.0)then 939 init=3 940 jnit=3 941 end if 942 943 f=0.d0 944 sig=x(ip) 945 var=sig*sig 946 vmere=0.d0 947 ngeno1=ngenom(current_chr,jm)+1 948 ngeno2=ngenom(current_chr,jm+1) 949 do ig=ngeno1,ngeno2 950 nd1=ngend(current_chr,ig)+1 951 nd2=ngend(current_chr,ig+1) 952 vpf=0.d0 953 effm=0.d0 954 do kd=nd1,nd2 955 kkd=ndesc(current_chr,kd) 956 957 if(presentc(ic,kkd)) then 958 effm=effm+1 959 v=y(ic,kkd) 960 ilev=corniv(nivdir(kkd,1)) 961 v=v-x(np+ilev) 962 if(vecsol(nivdir(kkd,2))) then 963 ilev=corniv(nivdir(kkd,2)) 964 v=v-x(np+ilev)*ppt(kd) 965 end if 966 967 if(estime(ic,jm)) then 968 if(vecsol(nivdir(kkd,3))) then 969 ilev=corniv(nivdir(kkd,3)) 970 v=v-x(np+ilev)*pmt(kd) 971 end if 972 end if 973 974 if(vecsol(nivdir(kkd,init))) then 975 ilev=corniv(nivdir(kkd,init)) 976 v=v-x(np+ilev) 977 end if 978 979 if(estime(ic,jm)) then 980 if(vecsol(nivdir(kkd,5))) then 981 ilev=corniv(nivdir(kkd,5)) 982 v=v-x(np+ilev) 983 end if 984 end if 985 986 do ief=jnit+1,jnit+nbef 987 988 if(vecsol(nivdir(kkd,ief))) then 989 ilev=corniv(nivdir(kkd,ief)) 990 v=v-x(np+ilev) 991 end if 992 end do 993 994 do ico=1,nbco 995 if(vecsol(ntnifix+ico)) then 996 ilev=corniv(ntnifix+ico) 997 v=v-covdir(kkd,ico)*x(np+ilev) 998 end if 999 end do 1000 vpf=vpf+v*v*cd(ic,kkd) 1001 end if 1002 end do 1003 vmere=vmere+probg(current_chr,ig)*dexp(-0.5d0*vpf/var) 1004 end do 1005 1006 if (vmere == 0) then 1007 f=INIFINY_REAL_VALUE 1008 else 1009 f=-dlog(vmere)+dble(effm)*dlog(sig) 1010 end if 1011 1012 end subroutine funct_1qtl_modlin_family
init_analyse_modlin
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
init_analyse_modlin
DESCRIPTION
Initialisation/allocation of solution/buffer arrays
NOTES
SOURCE
153 subroutine init_analyse_modlin 154 integer :: stat 155 156 allocate (sig1(np),STAT=stat) 157 call check_allocate(stat,'sig1 [m_qtlmap_analyse_modlin]') 158 allocate (xmu1p(np),STAT=stat) 159 call check_allocate(stat,'xmu1p [m_qtlmap_analyse_modlin]') 160 allocate (xmu1m(nm),STAT=stat) 161 call check_allocate(stat,'xmu1m [m_qtlmap_analyse_modlin]') 162 allocate (fp0(np),STAT=stat) 163 call check_allocate(stat,'fp0 [m_qtlmap_analyse_modlin]') 164 allocate (fp1(np),STAT=stat) 165 call check_allocate(stat,'fp1 [m_qtlmap_analyse_modlin]') 166 allocate (fm0(nm),STAT=stat) 167 call check_allocate(stat,'fm0 [m_qtlmap_analyse_modlin]') 168 allocate (fm1(nm),STAT=stat) 169 call check_allocate(stat,'fm1 [m_qtlmap_analyse_modlin]') 170 171 end subroutine init_analyse_modlin
opti_0qtl_modlin
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
opti_0qtl_modlin
DESCRIPTION
Calcul de la vraisemblance 0 QTL, 1 caractere , effets parasites inclus
NOTES
SOURCE
201 subroutine opti_0qtl_modlin(ic) 202 use m_qtlmap_analyse_lin_gen, only : contingence,precision,vecsol,precis, & 203 nbnivest,ntniv,par0,precis0,vecsol0,xx 204 205 implicit none 206 207 integer , intent(in) ::ic 208 ! 209 ! Divers 210 211 integer :: iuser(1) 212 integer ,dimension(:),allocatable :: iw 213 real (kind=dp) ,dimension(:),allocatable :: par,borni,borns,w 214 real (kind=dp) :: user(1) 215 real (kind=dp) :: f 216 integer :: npar,ibound,ip,ix,ifail,i,indest,ntot 217 integer :: km,jm,kd,kd1,kd2,ii 218 logical itest 219 logical , dimension(:,:,:),pointer :: filter_inc 220 221 ! 222 !****************************************************************************** 223 ! 224 ! recherche des elements estimables � partir de la d�composition de choleski 225 ! comptage et recodification de ces elements (on commence par mettre � 0 les 226 ! compteurs destin�s aux tests des effes fix�s 227 ! 228 itest=.false. 229 call contingence(ic,0,itest,.true.) 230 call precision(xx,precis) 231 ! Parametres de maximisation 232 npar=np+nbnivest 233 234 allocate (borni(npar)) 235 allocate (borns(npar)) 236 allocate (par(npar)) 237 238 !MODIF - OPMIZATION 239 allocate (filter_inc(np,nm, npar)) 240 call set_filter_optim(ic,.true.,.false.,ntnivmax,ntniv,vecsol,xinc,filter_inc) 241 !FIN MODIF - OPMIZATION 242 243 ibound=0 244 do ip=1,np 245 borni(ip)=SIG_MIN 246 borns(ip)=SIG_MAX 247 end do 248 do ix=np+1,npar 249 borni(ix)=XMU_MIN 250 borns(ix)=XMU_MAX 251 end do 252 253 ! 254 ! Point de depart 255 par(np+1)=0.d0 256 do ip=1,np 257 par(ip)=sig0(ip) 258 par(np+1)=par(np+1)+xmu0p(ip) 259 end do 260 par(np+1)=par(np+1)/dble(np) 261 do ix=np+2,npar 262 par(ix)=0.d0 263 end do 264 ! 265 ! Optimisation de la vraisemblance 266 ifail=1 267 current_ic = ic ! pour le mode lineaire.... 268 !call minimizing_funct(npar,ibound,funct_0qtl_modlin,borni,borns,par,f,iuser,user,ifail) 269 call minimizing_funct_family(npar,ibound,funct_0qtl_modlin_family,filter_inc,fm0,fp0,borni,borns,par,f,iuser,user,ifail) 270 271 f0=f 272 do i=1,npar 273 par0(i)=par(i) 274 end do 275 do i=1,ntniv 276 vecsol0(i)=vecsol(i) 277 precis0(i)=precis(i) 278 end do 279 280 do ip = 1,np 281 sig1(ip)=par(ip) 282 end do 283 284 xmu1g=par(np+1) 285 indest=1 286 do ip = 1,np 287 if (vecsol(1+ip)) then 288 indest=indest+1 289 xmu1p(ip)=par(np+indest) 290 end if 291 end do 292 293 ntot=np+1 294 km=0 295 do jm=1,nm 296 if (estime(ic,jm).and.opt_sib.eq.2)then 297 km=km+1 298 if(vecsol(ntot+km)) then 299 indest=indest+1 300 xmu1m(jm)=par(np+indest) 301 end if 302 end if 303 end do 304 305 deallocate (borni) 306 deallocate (borns) 307 deallocate (par) 308 deallocate (filter_inc) 309 310 return 311 end subroutine opti_0qtl_modlin
opti_1qtl_modlin
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
opti_1qtl_modlin
DESCRIPTION
Calcul de la statistique de test le long du chromosome
NOTES
SOURCE
527 subroutine opti_1qtl_modlin(ic,lrtsol,fmax,supnbnivest) 528 529 integer , intent(in) :: ic 530 type(TYPE_LRT_SOLUTION) , intent(out) :: lrtsol 531 real (kind=dp) , intent(out) :: fmax 532 integer , intent(out) :: supnbnivest 533 ! 534 535 ! Divers 536 real (kind=dp) ,dimension(:),allocatable :: val,par,borni,borns 537 real (kind=dp) :: user(1) 538 real (kind=dp) :: fm(nm),fp(np),f1,save_f0 539 540 logical itest,stvecsol(size(vecsol1)) 541 integer :: ip,i,n,ilong,ibound,npar,ix,j,ifail,iuser(1),indexchr(nchr),ntotal,nprime 542 integer :: ii,chr,nkd,jm,ig,ifem,iip,kd,kd1,kd2,save_nteffmax,save_ntnivmax 543 logical , dimension(:,:,:),pointer :: filter_inc 544 real(kind=dp) ,dimension(:,:) ,pointer :: listef 545 integer ,dimension(:,:) ,pointer :: listenbnivest 546 real(kind=dp) ,dimension(:,:,:) ,pointer :: listepar 547 real(kind=dp) ,dimension(:) ,pointer :: save_fm0 548 real(kind=dp) ,dimension(:) ,pointer :: save_fp0 549 real(kind=dp) ,dimension(np) :: save_sig1 550 551 552 !****************************************************************************** 553 ! Calcul de la vraisemblance sous H1 554 555 save_fp0=>fp0 556 save_fm0=>fm0 557 save_ntnivmax = ntnivmax 558 save_nteffmax = nteffmax 559 save_f0 = f0 560 save_sig1 = sig1 561 562 call new(1,lrtsol) 563 564 ! initialisation 565 ! on utilisera : 566 ! PAR (vecteur des param�tres � optimiser), 567 ! CORNIV (vecteur des positions, parmi NTNIV, des effets extimables), 568 ! SOLVEC (vecteur disant l'estimabilit� des effets) 569 ! VAL le vecteur complet (ntniv positions) des valeurs des niveaux des effets 570 ! STVAL,STVECSOL et STCORNIV les copies de VAL, VECSOL et CORNIV � l'it�ration 571 ! n-1 quand on analyse la position n 572 ! 573 ! 574 lrtsol%lrtmax=-1.d75 575 lrtsol%chrmax=0 576 lrtsol%nxmax=0 577 578 ntotal=0 579 do chr=1,nchr 580 ntotal=ntotal+get_npo(chr) 581 indexchr(chr)=ntotal 582 end do 583 584 allocate (listef(nchr,get_maxnpo())) 585 allocate (listenbnivest(nchr,get_maxnpo())) 586 allocate (listepar(nchr,get_maxnpo(),np + ntnivmax)) 587 !on libere l espace pour le cas non openmp 588 call end_contingence 589 ! 590 ! Marche le long du chromosome 591 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(par,val) & 592 !$OMP PRIVATE(filter_inc,i,ifem,iip,ii,stvecsol,npar,j,ifail,fp,fm,f1,n,chr,borni,borns) 593 current_ic = ic 594 ntnivmax = save_ntnivmax 595 nteffmax = save_nteffmax 596 fp0=>save_fp0 597 fm0=>save_fm0 598 f0 = save_f0 599 600 allocate ( val( np+ntnivmax ) ) 601 allocate ( par( np + ntnivmax) ) 602 allocate ( borni( np + ntnivmax) ) 603 allocate ( borns( np + ntnivmax) ) 604 605 ! Point de depart 606 ! 607 do ip=1,np 608 par(ip)=save_sig1(ip) 609 end do 610 611 par(np+1)=xmu1g 612 613 do i=np+2,size(par) 614 par(i)=0.d0 615 end do 616 617 do i=np+1,size(val) 618 val(i-np)=par(i) 619 end do 620 621 call init_contingence 622 vecsol=.true. 623 allocate (filter_inc(np,nm,ntnivmax+np)) 624 !$OMP DO 625 do nprime=1,ntotal 626 n=nprime 627 do chr=nchr,1,-1 628 if ( indexchr(chr) >= nprime) then 629 current_chr = chr 630 else 631 n=n-get_npo(chr) 632 end if 633 end do 634 chr=current_chr 635 636 ! 637 ! on stocke les conditions en n 638 do i=1,ntniv 639 stvecsol(i)=vecsol(i) 640 end do 641 ! 642 ! preparation de la matrice d'incidence 643 ! 644 call prepinc(current_chr,n,ic,"LA ") 645 ! 646 ! recherche des elements estimables � partir de la d�composition de choleski 647 ! comptage et recodification de ces elements (on commence par mettre � 0 les 648 ! compteurs destin�s aux tests des effes fix�s 649 ! 650 itest=.false. 651 call contingence(ic,1,itest) 652 !MODIF - OPMIZATION 653 call set_filter_optim(ic,.true.,.false.,ntnivmax,ntniv,vecsol,xinc,filter_inc) 654 !FIN MODIF - OPMIZATION 655 656 ! Parametres de maximisation 657 ibound=0 658 npar=np+nbnivest 659 660 do ip=1,np 661 borni(ip)=SIG_MIN 662 borns(ip)=SIG_MAX 663 end do 664 do i=np+1,npar 665 borni(i)=XMU_MIN 666 borns(i)=XMU_MAX 667 end do 668 ! 669 ! Point de depart (on reprend les point d'arriv�e pr�c�dents 670 ! 671 do i=np+1,npar 672 par(i)=0.d0 673 end do 674 675 j=np 676 do i=1,ntniv 677 if(vecsol(i))then 678 j=j+1 679 par(j)=0.d0 680 if(stvecsol(i)) par(j)=val(i) 681 end if 682 stvecsol(i)=vecsol(i) 683 end do 684 685 ! Optimisation de la vraisemblance a la position dx 686 ifail=1 687 688 call minimizing_funct_family(npar,ibound,funct_1qtl_modlin_family,filter_inc,fm,fp,borni,borns,par,f1,iuser,user,ifail) 689 690 j=np 691 do i=1,ntniv 692 if(vecsol(i)) then 693 j=j+1 694 val(i)=par(j) 695 else 696 val(i)=9999.d0 697 end if 698 end do 699 ! 700 ! on garde les valeurs du LRT pour pouvopir dessiner la courbe de vraisemblance 701 ! 702 if ( f1 < INIFINY_REAL_VALUE ) then 703 lrtsol%lrt1(chr,n)=-2.d0*(f1-f0) 704 else 705 lrtsol%lrt1(chr,n)=0 706 end if 707 708 listef(chr,n)=f1 709 listenbnivest(chr,n)=nbnivest 710 listepar(chr,n,:npar)=par(:npar) 711 ! 712 ! on met les profil / prog�niteur dnas les ficheir ad hoc 713 ! 714 iip=np+1 715 do ii=1,np 716 lrtsol%xlrp(chr,ii,n)=-2.d0*(fp(ii)-fp0(ii)) 717 if ( vecsol(1+ii)) then 718 iip=iip+1 719 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 720 lrtsol%pater_eff(chr,ii,n)=par(iip) 721 end if 722 end do 723 724 ifem=0 725 do ii=1,nm 726 lrtsol%xlrm(chr,ii,n)=-2.d0*(fm(ii)-fm0(ii)) 727 if ( estime(ic,ii)) then 728 ifem=ifem+1 729 if ( vecsol(np+1+ifem) ) then 730 iip=iip+1 731 ! ATTENTION VALIDE DANS LE CAS SANS INTERACTION AVEC QTL 732 lrtsol%mater_eff(chr,ifem,n)=par(iip) 733 end if 734 end if 735 end do 736 end do ! nprime 737 !$OMP END DO 738 call end_contingence 739 deallocate (filter_inc) 740 deallocate ( val ) 741 deallocate ( par ) 742 deallocate ( borni ) 743 deallocate ( borns ) 744 !$OMP END PARALLEL 745 746 !recherche du maximum 747 do chr=1,nchr 748 do n=1,get_npo(chr) 749 if(lrtsol%lrtmax(0) < lrtsol%lrt1(chr,n)) then 750 lrtsol%nxmax(0)=n 751 lrtsol%chrmax(0)=chr 752 lrtsol%lrtmax(0)=lrtsol%lrt1(chr,n) 753 fmax=listef(chr,n) 754 supnbnivest=listenbnivest(chr,n) 755 par1(:supnbnivest+np)=listepar(chr,n,:supnbnivest+np) 756 end if 757 end do 758 end do 759 deallocate (listef) 760 deallocate (listenbnivest) 761 deallocate (listepar) 762 763 call init_contingence 764 765 ! on calcule la pr�cision des estimation au point correspondant 766 ! au LRT maximum 767 ! 768 call prepinc(lrtsol%chrmax(0),lrtsol%nxmax(0),ic,"LA ") 769 call contingence(ic,1,itest) 770 call precision(xx,precis) 771 772 do i=1,ntniv 773 precis1(i)=precis(i) 774 vecsol1(i)=vecsol(i) 775 end do 776 777 return 778 end subroutine opti_1qtl_modlin
set_solution_hypothesis0
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
set_solution_hypothesis0
DESCRIPTION
NOTES
SOURCE
1194 subroutine set_solution_hypothesis0(ic,incsol) 1195 integer ,intent(in) :: ic 1196 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1197 1198 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 1199 integer :: indest,isol,ipar,nlevel,ief,i,ife 1200 1201 real(kind=dp) ,dimension(size(par0)) :: par0_t 1202 1203 ! remise a l echelle des variables 1204 do ipar=1,size(par0) 1205 par0_t(ipar)=par0(ipar)*sigt(ic) 1206 end do 1207 1208 allocate (incsol%sig(1,np)) 1209 1210 incsol%hypothesis=0 1211 ! Mean, Polygenic family 1212 nteff = 2 1213 if ( count(estime(ic,:)) > 0 ) nteff = 3 1214 1215 !Fixed effect and covariate 1216 if ( modele(ic,1) > 0 ) nteff = nteff+1 1217 if ( modele(ic,2) > 0 ) nteff = nteff+1 1218 1219 maxNbPar = max(np,count(estime(ic,:))) 1220 1221 !max numbre de covariable ? 1222 maxNbPar = max(maxNbPar,modele(ic,2)) 1223 1224 nlevel=0 1225 !max nombre de niveau pour un effet fixe ? 1226 do i=1,modele(ic,1) 1227 nlevel=nlev(modele(ic,3+i))+nlevel 1228 end do 1229 maxNbPar = max(maxNbPar,nlevel) 1230 1231 allocate (incsol%groupeName(nteff)) 1232 allocate (incsol%eqtl_print(nteff)) 1233 allocate (incsol%nbParameterGroup(nteff)) 1234 allocate (incsol%parameterName(nteff,maxNbPar)) 1235 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1236 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1237 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 1238 incsol%eqtl_print=.true. 1239 incsol%parameterPrecis=0.d0 1240 incsol%parameterVecsol=.true. 1241 1242 do ip=1,np 1243 incsol%sig(1,ip) = par0_t(ip) 1244 end do 1245 1246 ieff=1 1247 incsol%groupeName(ieff) = 'General Mean' 1248 incsol%nbParameterGroup(ieff)=1 1249 incsol%parameterName(ieff,1) ='General Mean' 1250 incsol%paramaterValue(ieff,1) = par0_t(np+1)+xmut(ic) 1251 incsol%parameterVecsol(ieff,1) = vecsol0(1) 1252 incsol%parameterPrecis(ieff,1) = precis0(1) 1253 1254 1255 ieff=ieff+1 1256 incsol%groupeName(ieff) = 'Sire polygenic effects' 1257 incsol%nbParameterGroup(ieff)=np 1258 1259 indest = np+1 1260 isol=1 1261 do ip=1,np 1262 isol=isol+1 1263 if ( vecsol0(isol) ) then 1264 indest=indest+1 1265 incsol%paramaterValue(ieff,ip) = par0_t(indest) 1266 else 1267 incsol%paramaterValue(ieff,ip) = 0.d0 1268 end if 1269 1270 incsol%parameterName(ieff,ip) ='Sire '//trim(pere(ip)) 1271 incsol%parameterVecsol(ieff,ip) = vecsol0(isol) 1272 incsol%parameterPrecis(ieff,ip) = precis0(isol) 1273 end do 1274 1275 if ( count(estime(ic,:)) > 0 ) then 1276 ieff = ieff +1 1277 incsol%groupeName(ieff) = 'Dam polygenic effects' 1278 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 1279 ifem=0 1280 do ip=1,np 1281 do jm=nmp(ip)+1,nmp(ip+1) 1282 if ( estime(ic,jm)) then 1283 isol=isol+1 1284 ifem=ifem+1 1285 if ( vecsol0(isol) ) then 1286 indest=indest+1 1287 incsol%paramaterValue(ieff,ifem) = par0_t(indest) 1288 else 1289 incsol%paramaterValue(ieff,ifem) = 0.d0 1290 end if 1291 1292 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1293 incsol%parameterVecsol(ieff,ifem) = vecsol0(isol) 1294 incsol%parameterPrecis(ieff,ifem) = precis0(isol) 1295 end if 1296 end do 1297 end do 1298 end if 1299 1300 !Fixed effect 1301 if ( modele(ic,1) > 0 ) then 1302 ieff=ieff+1 1303 incsol%eqtl_print(ieff)=.false. 1304 incsol%groupeName(ieff) = 'Fixed effects' 1305 incsol%nbParameterGroup(ieff)=nlevel 1306 ife=0 1307 do ief=1,nbef 1308 do i=1,nlev(modele(ic,3+ief)) 1309 ife=ife+1 1310 isol=isol+1 1311 1312 if ( vecsol0(isol) ) then 1313 indest=indest+1 1314 incsol%paramaterValue(ieff,ife) = par0_t(indest) 1315 else 1316 incsol%paramaterValue(ieff,ife) = 0.d0 1317 end if 1318 incsol%parameterName(ieff,ife) = trim(namefix(modele(ic,3+ief)))//' Level '//trim(str(i)) 1319 incsol%parameterVecsol(ieff,ife) = vecsol0(isol) 1320 incsol%parameterPrecis(ieff,ife) = precis0(isol) 1321 end do 1322 end do 1323 end if 1324 1325 !Covariate 1326 if ( modele(ic,2) > 0 ) then 1327 ieff=ieff+1 1328 incsol%eqtl_print(ieff)=.false. 1329 incsol%groupeName(ieff) = 'Covariates' 1330 incsol%nbParameterGroup(ieff)=modele(ic,2) 1331 1332 do ief=1,modele(ic,2) 1333 isol=isol+1 1334 if ( vecsol0(isol) ) then 1335 indest=indest+1 1336 incsol%paramaterValue(ieff,ief) = par0_t(indest) 1337 else 1338 incsol%paramaterValue(ieff,ief) = 0.d0 1339 end if 1340 incsol%parameterName(ieff,ief) = trim(namecov(modele(ic,3+modele(ic,1)+ief))) 1341 incsol%parameterVecsol(ieff,ief) = vecsol0(isol) 1342 incsol%parameterPrecis(ieff,ief) = precis0(isol) 1343 end do 1344 end if 1345 1346 end subroutine set_solution_hypothesis0
set_solution_hypothesis1
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
set_solution_hypothesis1
DESCRIPTION
NOTES
SOURCE
1358 subroutine set_solution_hypothesis1(ic,incsol) 1359 integer ,intent(in) :: ic 1360 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1361 1362 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 1363 integer :: ntlev,nbtp,jef,lp,indest,km 1364 integer :: isol,ipar,nlevel,ief,i,ife 1365 1366 real(kind=dp) ,dimension(size(par1)) :: par1_t 1367 1368 ! remise a l echelle des variables 1369 do ipar=1,size(par1) 1370 par1_t(ipar)=par1(ipar)*sigt(ic) 1371 end do 1372 1373 allocate (incsol%sig(1,np)) 1374 1375 incsol%hypothesis=1 1376 ! Mean, Polygenic family, QTL effect 1377 nteff = 3 1378 if ( count(estime(ic,:)) > 0 ) nteff = 5 1379 1380 !Fixed effect and covariate 1381 if ( modele(ic,1) > 0 ) nteff = nteff+1 1382 if ( modele(ic,2) > 0 ) nteff = nteff+1 1383 1384 maxNbPar = max(np,count(estime(ic,:))) 1385 !max numbre de covariable ? 1386 maxNbPar = max(maxNbPar,modele(ic,2)) 1387 1388 nlevel=0 1389 !max nombre de niveau pour un effet fixe ? 1390 do i=1,modele(ic,1) 1391 nlevel=nlev(modele(ic,3+i))+nlevel 1392 end do 1393 maxNbPar = max(maxNbPar,nlevel) 1394 1395 ntlev=1 1396 nbtp = 3 + modele(ic,1)+modele(ic,2)+modele(ic,3) 1397 do jef=1,modele(ic,3) 1398 ntlev=ntlev*nlev(modele(ic,nbtp+jef)) 1399 end do 1400 1401 !max nombre de niveau pour un effet fixe en interaction avec le qtl ? 1402 maxNbPar = max(maxNbPar,ntlev*np) 1403 maxNbPar = max(maxNbPar,ntlev*count(estime(ic,:))) 1404 1405 allocate (incsol%groupeName(nteff)) 1406 allocate (incsol%nbParameterGroup(nteff)) 1407 allocate (incsol%eqtl_print(nteff)) 1408 allocate (incsol%parameterName(nteff,maxNbPar)) 1409 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1410 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1411 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 1412 allocate (incsol%qtl_groupeName(1,1)) 1413 incsol%parameterPrecis=0.d0 1414 incsol%parameterVecsol=.true. 1415 incsol%eqtl_print=.true. 1416 1417 1418 do ip=1,np 1419 incsol%sig(1,ip) = par1_t(ip) 1420 end do 1421 1422 ieff=1 1423 incsol%groupeName(ieff) = 'General Mean' 1424 incsol%nbParameterGroup(ieff)=1 1425 incsol%parameterName(ieff,1) ='General Mean' 1426 incsol%paramaterValue(ieff,1) = par1_t(np+1)+xmut(ic) 1427 incsol%parameterVecsol(ieff,1) = vecsol1(1) 1428 incsol%parameterPrecis(ieff,1) = precis1(1) 1429 1430 ieff=ieff+1 1431 incsol%qtl_groupeName(1,1)=ieff 1432 incsol%groupeName(ieff) = 'Sire QTL effects' 1433 incsol%nbParameterGroup(ieff)=np*ntlevp 1434 1435 indest = np+1 1436 isol=1 1437 ife=0 1438 do ip=1,np 1439 do lp=1,ntlev 1440 ife=ife+1 1441 isol=isol+1 1442 if ( vecsol1(isol) ) then 1443 indest=indest+1 1444 incsol%paramaterValue(ieff,ife) = par1_t(indest) 1445 else 1446 incsol%paramaterValue(ieff,ife) = 0.d0 1447 end if 1448 1449 incsol%parameterName(ieff,ife) ='Sire '//trim(pere(ip))//" "//trim(str(lp)) 1450 incsol%parameterVecsol(ieff,ife) = vecsol1(isol) 1451 incsol%parameterPrecis(ieff,ife) = precis1(isol) 1452 end do 1453 end do 1454 1455 if ( count(estime(ic,:)) > 0 ) then 1456 ieff = ieff +1 1457 incsol%groupeName(ieff) = 'Dam QTL effects' 1458 incsol%nbParameterGroup(ieff)=namest(ic)*ntlevp 1459 ifem=0 1460 do jm=1,nm 1461 if ( estime(ic,jm)) then 1462 do lp=1,ntlev 1463 isol=isol+1 1464 ifem=ifem+1 1465 if ( vecsol1(isol) ) then 1466 indest=indest+1 1467 incsol%paramaterValue(ieff,ifem) = par1_t(indest) 1468 else 1469 incsol%paramaterValue(ieff,ifem) = 0.d0 1470 end if 1471 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" "//trim(str(lp)) 1472 incsol%parameterVecsol(ieff,ifem) = vecsol1(isol) 1473 incsol%parameterPrecis(ieff,ifem) = precis1(isol) 1474 end do 1475 end if 1476 end do 1477 end if 1478 1479 1480 ieff=ieff+1 1481 incsol%groupeName(ieff) = 'Sire polygenic effects' 1482 incsol%nbParameterGroup(ieff)=np 1483 1484 do ip=1,np 1485 isol=isol+1 1486 if ( vecsol1(isol) ) then 1487 indest=indest+1 1488 incsol%paramaterValue(ieff,ip) = par1_t(indest) 1489 else 1490 incsol%paramaterValue(ieff,ip) = 0.d0 1491 end if 1492 1493 incsol%parameterName(ieff,ip) ='Sire '//trim(pere(ip)) 1494 incsol%parameterVecsol(ieff,ip) = vecsol1(isol) 1495 incsol%parameterPrecis(ieff,ip) = precis1(isol) 1496 end do 1497 1498 if ( count(estime(ic,:)) > 0 ) then 1499 ieff = ieff +1 1500 incsol%groupeName(ieff) = 'Dam polygenic effects' 1501 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 1502 ifem=0 1503 do ip=1,np 1504 do jm=nmp(ip)+1,nmp(ip+1) 1505 if ( estime(ic,jm)) then 1506 isol=isol+1 1507 ifem=ifem+1 1508 if ( vecsol1(isol) ) then 1509 indest=indest+1 1510 incsol%paramaterValue(ieff,ifem) = par1_t(indest) 1511 else 1512 incsol%paramaterValue(ieff,ifem) = 0.d0 1513 end if 1514 1515 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1516 incsol%parameterVecsol(ieff,ifem) = vecsol1(isol) 1517 incsol%parameterPrecis(ieff,ifem) = precis1(isol) 1518 end if 1519 end do 1520 end do 1521 end if 1522 1523 !Fixed effect 1524 if ( modele(ic,1) > 0 ) then 1525 ieff=ieff+1 1526 incsol%eqtl_print(ieff)=.false. 1527 incsol%groupeName(ieff) = 'Fixed effects' 1528 incsol%nbParameterGroup(ieff)=nlevel 1529 ife=0 1530 do ief=1,nbef 1531 do i=1,nlev(modele(ic,3+ief)) 1532 isol=isol+1 1533 ife=ife+1 1534 if ( vecsol1(isol) ) then 1535 indest=indest+1 1536 incsol%paramaterValue(ieff,ife) = par1_t(indest) 1537 else 1538 incsol%paramaterValue(ieff,ife) = 0.d0 1539 end if 1540 incsol%parameterName(ieff,ife) = trim(namefix(modele(ic,3+ief)))//' Level '//trim(str(i)) 1541 incsol%parameterVecsol(ieff,ife) = vecsol1(isol) 1542 incsol%parameterPrecis(ieff,ife) = precis1(isol) 1543 end do 1544 end do 1545 end if 1546 1547 !Covariate 1548 if ( modele(ic,2) > 0 ) then 1549 ieff=ieff+1 1550 incsol%eqtl_print(ieff)=.false. 1551 incsol%groupeName(ieff) = 'Covariates' 1552 incsol%nbParameterGroup(ieff)=modele(ic,2) 1553 1554 do ief=1,modele(ic,2) 1555 isol=isol+1 1556 if ( vecsol1(isol) ) then 1557 indest=indest+1 1558 incsol%paramaterValue(ieff,ief) = par1_t(indest) 1559 else 1560 incsol%paramaterValue(ieff,ief) = 0.d0 1561 end if 1562 incsol%parameterName(ieff,ief) = trim(namecov(modele(ic,3+modele(ic,1)+ief))) 1563 incsol%parameterVecsol(ieff,ief) = vecsol1(isol) 1564 incsol%parameterPrecis(ieff,ief) = precis1(isol) 1565 end do 1566 end if 1567 1568 end subroutine set_solution_hypothesis1
test_lin
[ Top ] [ m_qtlmap_analyse_modlin ] [ Subroutines ]
NAME
test_lin
DESCRIPTION
Test des differents effets de nuisance du modele par une LRT compare a une chi2
NOTES
SOURCE
1023 subroutine test_lin(chr,ic,est_moy,supnbnivest,fmax,nposx) 1024 use m_qtlmap_analyse_lin_gen, only : prepinc,contingence,nbef,nbco,meff,mcov,mint,nbnivest 1025 1026 logical , intent(in) :: est_moy 1027 integer , intent(in) :: chr,ic 1028 integer , intent(in) :: supnbnivest 1029 real (kind=dp) , intent(in) :: fmax 1030 integer , intent(in) :: nposx 1031 1032 ! 1033 ! Divers 1034 logical itest 1035 integer iuser(1) 1036 real (kind=dp) ,dimension(:) , allocatable :: par,borni,borns 1037 real (kind=dp) :: user(1),prob,xlrt_t,f1 1038 integer :: nbint,iecd,iecq,ief,ibound,npar,ip,i,ifail 1039 integer :: nbreduit 1040 logical , dimension(:,:,:),pointer :: filter_inc 1041 1042 write(nficout,600) 1043 600 format(//,80('*')/'test of the effets of the model',//) 1044 1045 ! 1046 ! preparation de la matrice d'incidence 1047 ! 1048 call prepinc(chr,nposx,ic,"LA ") 1049 ! 1050 ! on recalcule la vraisemblance en retirant les effets du mod�les un � un 1051 ! 1052 ! 1053 ! on commence par les effets hierarchis�s dans les effets qtl 1054 ! 1055 current_chr = chr 1056 nbef=modele(ic,1) 1057 nbco=modele(ic,2) 1058 nbint=modele(ic,3) 1059 1060 iecd=0 1061 iecq=0 1062 1063 allocate ( par( np + ntnivmax ) ) 1064 allocate ( borni( np + ntnivmax ) ) 1065 allocate ( borns( np + ntnivmax ) ) 1066 !MODIF - OPMIZATION 1067 allocate (filter_inc(np,nm, np + ntnivmax )) 1068 !FIN MODIF - OPMIZATION 1069 1070 1071 do ief=1,nbint+nbef+nbco 1072 meff=0 1073 mcov=0 1074 mint=0 1075 if(ief.le.nbef)meff=ief 1076 if(ief.gt.nbef.and.ief.le.(nbef+nbco))mcov=ief-nbef 1077 if(ief.gt.(nbef+nbco))mint=ief-(nbef+nbco) 1078 ! 1079 ! recherche des elements estimables � partir de la d�composition de choleski 1080 ! comptage et recodification de ces elements 1081 ! 1082 1083 itest=.true. 1084 call contingence(ic,1,itest,est_moy) 1085 call set_filter_optim(ic,.true.,.false.,ntnivmax,ntniv,vecsol,xinc,filter_inc) 1086 1087 ! 1088 ! la r�dustion du nombre d'effet estim�e est calcul�e 1089 ! 1090 nbreduit=supnbnivest-nbnivest 1091 ! 1092 ! Parametres de maximisation 1093 ibound=0 1094 npar=np+nbnivest 1095 1096 do ip=1,np 1097 borni(ip)=1.do-6 1098 borns(ip)=1.d6 1099 end do 1100 do i=np+1,npar 1101 borni(i)=-1d6 1102 borns(i)=1.d6 1103 par(i)=0.d0 1104 end do 1105 ! 1106 ! Point de depart (on reprend les point d'arriv�e pr�c�dents 1107 ! 1108 do i=1,np 1109 par(i)=par1(i) 1110 end do 1111 1112 do i=np+1,npar 1113 par(i)=0.d0 1114 end do 1115 ! 1116 ! Optimisation de la vraisemblance a la position dx 1117 ifail=1 1118 1119 ! call minimizing_funct(npar,ibound,funct_1qtl_modlin,borni,borns,par,f1,iuser,user,ifail) 1120 call minimizing_funct_family(npar,ibound,funct_1qtl_modlin_family,filter_inc,fm1,fp1,borni,borns,par,f1,iuser,user,ifail) 1121 ! if (ifail.ne.0) print *,'Code retour e04jyf H1 : ',ifail 1122 ! 1123 xlrt_t=-2.d0*(fmax-f1) 1124 if(xlrt_t.le.1.do-8)xlrt_t=0.d0 1125 ! 1126 ! impression des tests des effets du mod�le 1127 ! 1128 1129 nbef=modele(ic,1) 1130 nbco=modele(ic,2) 1131 nbint=modele(ic,3) 1132 1133 if(ief.le.(nbef+nbco).and.iecd.eq.0) then 1134 iecd=1 1135 write(nficout,610) 1136 610 format(' Direct effects'/) 1137 1138 write(nficout,601) 1139 601 format('Tested effect df. Likelihood p-value'/ & 1140 ' ratio '/) 1141 end if 1142 1143 if(ief.le.nbef) then 1144 ifail=0 1145 ! prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nlev(modele(ic,3+ief))),ifail) 1146 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nbreduit),ifail) 1147 1148 write(nficout,602) trim(namefix(modele(ic,3+ief))),nbreduit,xlrt_t,prob 1149 1150 end if 1151 602 format(a15,3x,i2,6x,f8.2,7x,f5.3) 1152 1153 if(ief.gt.nbef.and.ief.le.(nbef+nbco)) then 1154 ifail=0 1155 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(1.0),ifail) 1156 1157 write(nficout,603) trim(namecov(modele(ic,3+ief))),xlrt_t,prob 1158 1159 end if 1160 603 format(a15,4x,'1',6x,f8.2,7x,f5.3) 1161 1162 if(ief.gt.(nbef+nbco).and.iecq.eq.0)then 1163 iecq=1 1164 write(nficout,611) 1165 611 format(/' Intra qtl effects'/) 1166 write(nficout,601) 1167 end if 1168 1169 if(ief.gt.(nbef+nbco)) then 1170 ifail=0 1171 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nbreduit),ifail) 1172 write(nficout,602) trim(namefix(modele(ic,3+ief))),nbreduit,xlrt_t,prob 1173 ! write(*,602) trim(namefix(modele(ic,3+ief))),nbreduit,xlrt_t,prob 1174 end if 1175 1176 end do 1177 1178 deallocate ( par ) 1179 deallocate ( borni ) 1180 deallocate ( borns ) 1181 deallocate (filter_inc) 1182 1183 return 1184 end subroutine test_lin