m_qtlmap_analyse_modlin_ldla
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_analyse_modlin_ldla
DESCRIPTION
NOTES
SEE ALSO
f0
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
f0
DESCRIPTION
value of the likelihood under H0
fm0
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
fm0
DESCRIPTION
value of the likelihood by full-sib family under H0
DIMENSIONS
nm
fm1
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
fm1
DESCRIPTION
value of the likelihood by full-sib family under H1
DIMENSIONS
nm
fp0
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
fp0
DESCRIPTION
value of the likelihood by sire family under H0
DIMENSIONS
np
fp1
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
fp1
DESCRIPTION
value of the likelihood by sire family under H1
DIMENSIONS
np
sig1
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
sig1
DESCRIPTION
The standart deviation under H0
DIMENSIONS
np
xmu1g
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
xmu1g
DESCRIPTION
The genral mean
xmu1m
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
xmu1m
DESCRIPTION
The polygenic mean for each full sib family
DIMENSIONS
nm
xmu1p
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Variables ]
NAME
xmu1p
DESCRIPTION
The polygenic mean for each sire family under H0
DIMENSIONS
np
end_analyse_modlin_ldla
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
end_analyse_modlin_ldla
DESCRIPTION
SOURCE
178 subroutine end_analyse_modlin_ldla 179 deallocate (sig1) 180 deallocate (xmu1p) 181 deallocate (xmu1m) 182 deallocate (fp0) 183 deallocate (fp1) 184 deallocate (fm0) 185 deallocate (fm1) 186 end subroutine end_analyse_modlin_ldla
funct_0qtl_modlin_ldla_family
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
funct_0qtl_modlin_ldla_family
DESCRIPTION
NOTES
Calcul de la vraisemblance et de ses derivees partielles sous H0
SOURCE
443 subroutine funct_0qtl_modlin_ldla_family(ip,jm,n,x,f,iuser,user) 444 use m_qtlmap_analyse_lin_gen, only : nbfem,nbef,ntnifix,nbco,covdir,nivdir,corniv,vecsol 445 446 implicit none 447 integer , intent(in) :: ip,jm,n 448 ! 449 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 450 ! 451 real (kind=dp) ,dimension(n), intent(in) :: x 452 real (kind=dp) , intent(inout) :: f 453 integer , dimension(2) :: iuser 454 real (kind=dp) ,dimension(1), intent(in) :: user 455 456 ! Tableaux dimensionnes selon nm, le nombre de meres 457 real (kind=dp) ,dimension(nm) :: effm 458 459 integer ::jnit,nm1,nm2,nd1,nd2,kkd,ilev,ief,ico,imoy 460 integer :: nb_var 461 real (kind=dp) :: sig,var,vmere,vpf,v 462 ! 463 464 !****************************************************************************** 465 466 !****************** 467 ! passage des arguments par iuser 468 ! iuser(1) = nb_var 469 ! iuser(2) = imoy 470 nb_var = iuser(1) 471 imoy = iuser(2) 472 473 jnit=3 474 if (nbfem.eq.0)jnit=2 475 if(iuser(1)== 0) jnit=jnit-1 476 477 478 f=0.d0 479 if(variance_homo) then 480 sig=x(1) 481 var=sig*sig 482 else 483 sig=x(ip) 484 var=sig*sig 485 end if 486 487 do kkd=ndm(jm)+1,ndm(jm+1) 488 ! 489 if(presentc(mod_ic,kkd) ) then 490 v=y(mod_ic,kkd) 491 ! 492 if(imoy ==1 .and. vecsol(nivdir(kkd,1))) v=v-x(nb_var+corniv(nivdir(kkd,1))) 493 ! 494 if(vecsol(nivdir(kkd,imoy+1))) v=v-x(nb_var+corniv(nivdir(kkd,imoy+1))) 495 ! 496 if(estime(mod_ic,jm)) then 497 if(vecsol(nivdir(kkd,imoy+2))) v=v-x(nb_var+corniv(nivdir(kkd,imoy+2))) 498 end if 499 500 do ief=jnit+1,jnit+nbef 501 if(vecsol(nivdir(kkd,ief))) v=v-x(nb_var+corniv(nivdir(kkd,ief))) 502 end do 503 504 do ico=1,nbco 505 if(vecsol(ntnifix+ico)) v=v-covdir(kkd,ico)*x(nb_var+corniv(ntnifix+ico)) 506 end do 507 508 f=f+v*v*cd(mod_ic,kkd)/var + dlog(var) 509 end if!presentc 510 end do!kkd 511 return 512 end subroutine funct_0qtl_modlin_ldla_family
funct_1qtl_modlin_ldla
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
funct_1qtl_modlin_ldla
DESCRIPTION
NOTES
Calcul de la vraisemblance et de ses derivees partielles sous Hypoth�se 1QTL 1carac
SOURCE
860 subroutine funct_1qtl_modlin_ldla(n,x,f,iuser,user) 861 implicit none 862 863 integer , intent(in) :: n 864 ! 865 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 866 ! 867 real (kind=dp) ,dimension(n), intent(in) :: x 868 real (kind=dp) , intent(inout) :: f 869 integer , dimension(8) :: iuser 870 real (kind=dp) ,dimension(1), intent(in) :: user 871 872 ! Tableaux dimensionnes selon nm, le nombre de meres 873 real (kind=dp) ,dimension(nm) :: effm 874 875 876 ! 877 ! Divers 878 real (kind=dp) :: sig,var,vmere,vpf,v 879 real (kind=dp) :: sigc,varc,vart,det 880 integer :: ip,nm1,nm2,jm,ngeno1,ngeno2,ig 881 integer :: nd1,nd2,kd,kkd,ilev,ief,ico 882 integer :: i_haplo 883 integer :: imp_funct,i 884 885 !****************************************************************************** 886 ! 887 !****************** 888 ! passage des arguments par iuser 889 ! iuser(1) = nb_var 890 ! iuser(2) = n 891 ! iuser(3) = imoy 892 ! iuser(4) = imoyld 893 ! iuser(5) = init 894 ! iuser(6) = jnit 895 896 897 f=0.d0 898 if(variance_homo) then 899 sig=x(1) 900 var=sig*sig 901 end if 902 903 do ip=1,np 904 if(.not. variance_homo) then 905 sig=x(ip) 906 var=sig*sig 907 end if 908 909 fp1(ip)=0.d0 910 do jm=nmp(ip)+1,nmp(ip+1) 911 vmere=0.d0 912 det=0.d0 913 do ig=ngenom(current_chr,jm)+1,ngenom(current_chr,jm+1) 914 vpf=0.d0 915 do kd=ngend(current_chr,ig)+1,ngend(current_chr,ig+1) 916 kkd=ndesc(current_chr,kd) 917 918 if(presentc(mod_ic,kkd)) then 919 920 v=y(mod_ic,kkd) 921 if(estime_moy) v=v-x(mod_nb_var+corniv(nivdir(kkd,1))) 922 ! 923 ! effets haplotypes 924 if(is_ld .or. is_ldla ) then 925 do i_haplo=1,nb_haplo_reduit 926 if(vecsol(nivdir(kkd,mod_imoy+i_haplo))) & 927 v=v-pb_haplo(kkd,i_haplo)*x(mod_nb_var+corniv(nivdir(kkd,mod_imoy+i_haplo))) 928 end do 929 end if ! option 930 931 if(.not. is_ld)then 932 ! effet qtl pere 933 if(vecsol(nivdir(kkd,mod_imoyld+1))) v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_imoyld+1)))*ppt(kd) 934 ! 935 ! effet qtl mere 936 if(.not. is_ldjh) then 937 if (estime(mod_ic,jm)) then 938 if (vecsol(nivdir(kkd,mod_imoyld+2))) v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_imoyld+2)))*pmt(kd) 939 940 end if 941 end if 942 ! effet qtl mere ou qtl pere transmis par la mere (modele JH) 943 944 if(is_ldjh .and. iuser(7) ==1) then 945 if(vecsol(nivdir(kkd,mod_imoy+2)).and. nivdir(kkd,mod_imoy+2) /= 0) & 946 v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_imoy+2)))*(1-2*modulo(shared_haplo(kkd),2)) 947 end if 948 end if ! option 949 ! 950 ! effet polygenique pere 951 if(vecsol(nivdir(kkd,mod_init))) v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_init))) 952 ! 953 ! effet polygenique mere 954 if(estime(mod_ic,jm) ) then 955 if ( vecsol(nivdir(kkd,mod_init+1))) v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_init+1))) 956 end if 957 ! 958 ! effets fixes 959 do ief=mod_jnit+1,mod_jnit+nbef 960 if(vecsol(nivdir(kkd,ief))) v=v-x(mod_nb_var+corniv(nivdir(kkd,ief))) 961 end do 962 ! 963 ! covariables 964 do ico=1,nbco 965 if(vecsol(ntnifix+ico)) v=v-covdir(kkd,ico)*x(mod_nb_var+corniv(ntnifix+ico)) 966 end do 967 968 vpf=vpf+v*v*cd(mod_ic,kkd)/var 969 det=det+dlog(var) 970 end if !presentc 971 end do ! kd 972 973 vmere=vmere+probg(current_chr,ig)*dexp(-0.5d0*vpf) 974 end do !ig 975 976 977 if (vmere == 0) then 978 fm1(jm)=INIFINY_REAL_VALUE 979 else 980 fm1(jm)=-2.d0*dlog(vmere)+det 981 end if 982 983 fp1(ip)=fp1(ip)+fm1(jm) 984 end do !jm 985 f=f+fp1(ip) 986 end do !ip 987 988 return 989 end subroutine funct_1qtl_modlin_ldla
funct_1qtl_modlin_ldla_family
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
funct_1qtl_modlin_ldla_family
DESCRIPTION
NOTES
SOURCE
1001 subroutine funct_1qtl_modlin_ldla_family(ip,jm,n,x,f,iuser,user) 1002 implicit none 1003 1004 integer , intent(in) :: n,ip,jm 1005 ! 1006 ! Tableaux dimensionnes selon n, le nombre de parametres a estimer 1007 ! 1008 real (kind=dp) ,dimension(n), intent(in) :: x 1009 real (kind=dp) , intent(inout) :: f 1010 integer , dimension(8) :: iuser 1011 real (kind=dp) ,dimension(1), intent(in) :: user 1012 1013 1014 ! 1015 ! Divers 1016 real (kind=dp) :: sig,var,vmere,vpf,v 1017 real (kind=dp) :: sigc,varc,vart,det 1018 integer :: nm1,nm2,ngeno1,ngeno2,ig 1019 integer :: nd1,nd2,kd,kkd,ilev,ief,ico 1020 integer :: i_haplo 1021 integer :: imp_funct,i 1022 !****************************************************************************** 1023 ! 1024 !****************** 1025 ! passage des arguments par iuser 1026 ! iuser(1) = nb_var 1027 ! iuser(2) = n 1028 ! iuser(3) = imoy 1029 ! iuser(4) = imoyld 1030 ! iuser(5) = init 1031 ! iuser(6) = jnit 1032 1033 f=0.d0 1034 if(variance_homo) then 1035 sig=x(1) 1036 var=sig*sig 1037 end if 1038 1039 1040 if(.not. variance_homo) then 1041 sig=x(ip) 1042 var=sig*sig 1043 end if 1044 1045 vmere=0.d0 1046 det=0.d0 1047 do ig=ngenom(current_chr,jm)+1,ngenom(current_chr,jm+1) 1048 vpf=0.d0 1049 do kd=ngend(current_chr,ig)+1,ngend(current_chr,ig+1) 1050 kkd=ndesc(current_chr,kd) 1051 1052 if(presentc(mod_ic,kkd)) then 1053 1054 v=y(mod_ic,kkd) 1055 if(estime_moy) v=v-x(mod_nb_var+corniv(nivdir(kkd,1))) 1056 ! 1057 ! effets haplotypes 1058 if(is_ld .or. is_ldla ) then 1059 do i_haplo=1,nb_haplo_reduit 1060 if(vecsol(nivdir(kkd,mod_imoy+i_haplo))) & 1061 v=v-pb_haplo(kkd,i_haplo)*x(mod_nb_var+corniv(nivdir(kkd,mod_imoy+i_haplo))) 1062 end do 1063 end if ! option 1064 1065 if(.not. is_ld)then 1066 ! effet qtl pere 1067 if(vecsol(nivdir(kkd,mod_imoyld+1))) v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_imoyld+1)))*ppt(kd) 1068 ! 1069 ! effet qtl mere 1070 if(.not. is_ldjh) then 1071 if (estime(mod_ic,jm)) then 1072 if (vecsol(nivdir(kkd,mod_imoyld+2))) v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_imoyld+2)))*pmt(kd) 1073 1074 end if 1075 end if 1076 ! effet qtl mere ou qtl pere transmis par la mere (modele JH) 1077 1078 if(is_ldjh .and. iuser(7) ==1) then 1079 if(vecsol(nivdir(kkd,mod_imoy+2)).and. nivdir(kkd,mod_imoy+2) /= 0) & 1080 v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_imoy+2)))*(1-2*modulo(shared_haplo(kkd),2)) 1081 end if 1082 end if ! option 1083 ! 1084 ! effet polygenique pere 1085 if(vecsol(nivdir(kkd,mod_init))) v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_init))) 1086 ! 1087 ! effet polygenique mere 1088 if(estime(mod_ic,jm) ) then 1089 if ( vecsol(nivdir(kkd,mod_init+1))) v=v-x(mod_nb_var+corniv(nivdir(kkd,mod_init+1))) 1090 end if 1091 ! 1092 ! effets fixes 1093 do ief=mod_jnit+1,mod_jnit+nbef 1094 if(vecsol(nivdir(kkd,ief))) v=v-x(mod_nb_var+corniv(nivdir(kkd,ief))) 1095 end do 1096 ! 1097 ! covariables 1098 do ico=1,nbco 1099 if(vecsol(ntnifix+ico)) v=v-covdir(kkd,ico)*x(mod_nb_var+corniv(ntnifix+ico)) 1100 end do 1101 1102 vpf=vpf+v*v*cd(mod_ic,kkd)/var 1103 det=det+dlog(var) 1104 end if !presentc 1105 end do ! kd 1106 1107 vmere=vmere+probg(current_chr,ig)*dexp(-0.5d0*vpf) 1108 end do !ig 1109 1110 1111 if (vmere == 0) then 1112 f=INIFINY_REAL_VALUE 1113 else 1114 f=-2.d0*dlog(vmere)+det 1115 end if 1116 1117 return 1118 end subroutine funct_1qtl_modlin_ldla_family
init_analyse_modlin_ldla
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
init_analyse_modlin_ldla
DESCRIPTION
SOURCE
151 subroutine init_analyse_modlin_ldla 152 integer :: stat 153 154 allocate (sig1(np),STAT=stat) 155 call check_allocate(stat,'sig1 [m_qtlmap_analyse_modlin_ldla]') 156 allocate (xmu1p(np),STAT=stat) 157 call check_allocate(stat,'xmu1p [m_qtlmap_analyse_modlin_ldla]') 158 allocate (xmu1m(nm),STAT=stat) 159 call check_allocate(stat,'xmu1m [m_qtlmap_analyse_modlin_ldla]') 160 allocate (fp0(np),STAT=stat) 161 call check_allocate(stat,'fp0 [m_qtlmap_analyse_modlin_ldla]') 162 allocate (fp1(np),STAT=stat) 163 call check_allocate(stat,'fp1 [m_qtlmap_analyse_modlin_ldla]') 164 allocate (fm0(nm),STAT=stat) 165 call check_allocate(stat,'fm0 [m_qtlmap_analyse_modlin_ldla]') 166 allocate (fm1(nm),STAT=stat) 167 call check_allocate(stat,'fm1 [m_qtlmap_analyse_modlin_ldla]') 168 169 end subroutine init_analyse_modlin_ldla
opti_0qtl_modlin_ldla
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
opti_0qtl_modlin_ldla
DESCRIPTION
NOTES
Calcul de la vraisemblance 0 QTL, 1 caractere , effets parasites inclus
SOURCE
197 subroutine opti_0qtl_modlin_ldla(ic,est_moy,option_var) 198 use m_qtlmap_analyse_lin_gen, only : contingence,precision,vecsol,precis, & 199 nbnivest,ntniv,par0,precis0,vecsol0,xx 200 201 implicit none 202 203 integer , intent(in) :: ic 204 logical , intent(in) :: est_moy 205 character(len=4) , intent(in) :: option_var 206 ! 207 ! Divers 208 209 integer :: iuser(2) 210 integer ,dimension(:),allocatable :: iw 211 real (kind=dp) ,dimension(:),allocatable :: par,borni,borns,w 212 real (kind=dp) :: user(1) 213 real (kind=dp) :: f,var_yQy 214 logical , dimension(:,:,:),allocatable :: filter_inc 215 216 integer :: npar,ibound,ip,ix,ifail,i,indest,ntot 217 integer :: km,jm,imoy 218 integer :: nb_var,nparx 219 logical itest,details 220 ! 221 !****************************************************************************** 222 ! 223 ! recherche des elements estimables � partir de la d�composition de choleski 224 ! comptage et recodification de ces elements (on commence par mettre � 0 les 225 ! compteurs destin�s aux tests des effes fix�s 226 ! 227 nb_var=1 228 variance_homo = (option_var == 'homo') 229 230 if(option_var == 'hete') nb_var=np 231 ! if(option_anal == 'LD ' .or. option_anal == 'LDLA')nb_var=2*nb_var 232 233 estime_moy = est_moy 234 235 itest=.false. 236 details=.false. 237 call contingence_ldla(ic,0,itest,1,0,var_yQy,details,est_moy,"LA ",.false.,.false.) 238 call precision(xx,precis) 239 240 !****************** 241 ! passage des arguments par iuser 242 ! iuser(1) = nb_var 243 ! iuser(2) = imoy 244 iuser(1) = nb_var 245 iuser(2) = 0 246 if(est_moy) iuser(2) = 1 247 imoy=iuser(2) 248 249 250 ibound=0 251 npar=nb_var+nbnivest 252 253 allocate (borni(npar)) 254 allocate (borns(npar)) 255 allocate (par(npar)) 256 257 258 !MODIF - OPMIZATION 259 allocate (filter_inc(np,nm, npar)) 260 call set_filter_optim(ic,(option_var == 'hete'),.false.,ntnivmax,ntniv,vecsol,xinc,filter_inc) 261 !FIN MODIF - OPMIZATION 262 263 do ip=1,nb_var 264 borni(ip)=SIG_MIN 265 borns(ip)=SIG_MAX 266 end do 267 do ix=nb_var+1,npar 268 borni(ix)=XMU_MIN 269 borns(ix)=XMU_MAX 270 end do 271 ! 272 ! Point de depart 273 par=0.d0 274 par(np+1)=0.d0 275 do ip=1,np 276 if (option_var == 'hete' ) then 277 par(ip)=sig0(ip) 278 else 279 par(1)=sig0(ip)+par(1) 280 end if 281 par(np+1)=par(np+1)+xmu0p(ip) 282 end do 283 par(np+1)=par(np+1)/dble(np) 284 if (option_var == 'homo' ) par(1) = par(1) / dble(np) 285 do ix=nb_var+2,npar 286 par(ix)=0.d0 287 end do 288 ! filter_inc=.true. 289 ! 290 ! Optimisation de la vraisemblance 291 ifail=1 292 mod_ic = ic ! pour le mode lineaire.... 293 ! call minimizing_funct(npar,ibound,funct_0qtl_modlin_ldla,borni,borns,par,f,iuser,user,ifail) 294 call minimizing_funct_family(npar,ibound,funct_0qtl_modlin_ldla_family,filter_inc,fm0,fp0,borni,borns,par,f,iuser,user,ifail) 295 ! if (ifail.ne.0) print *,'Code retour optimizing 0 QTL : ',ifail 296 297 f0=f 298 do i=1,npar 299 par0(i)=par(i) 300 end do 301 do i=1,ntniv 302 vecsol0(i)=vecsol(i) 303 precis0(i)=precis(i) 304 end do 305 306 if ( option_var == 'hete' ) then 307 do ip = 1,np 308 sig1(ip)=par(ip) 309 end do 310 else 311 sig1 =par(1) 312 end if 313 314 xmu1g=par(nb_var+1) 315 indest=1 316 do ip = 1,np 317 if (vecsol(1+ip)) then 318 indest=indest+1 319 xmu1p(ip)=par(nb_var+indest) 320 end if 321 end do 322 323 ntot=nb_var+1 324 km=0 325 do jm=1,nm 326 if (estime(ic,jm))then 327 km=km+1 328 if(vecsol(ntot+km)) then 329 indest=indest+1 330 xmu1m(jm)=par(nb_var+indest) 331 end if 332 end if 333 end do 334 335 deallocate (borni) 336 deallocate (borns) 337 deallocate (par) 338 deallocate (filter_inc) 339 340 return 341 end subroutine opti_0qtl_modlin_ldla
opti_1qtl_modlin_ldla
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
opti_1qtl_modlin_ldla
DESCRIPTION
Calcul de la statistique de test le long du chromosome
NOTES
SOURCE
524 subroutine opti_1qtl_modlin_ldla(ic,lrtsol,fmax,supnbnivest,est_moy,option_var,option_anal,hdam) 525 integer , intent(in) :: ic 526 type(TYPE_LRT_SOLUTION) , intent(out) :: lrtsol 527 real (kind=dp) , intent(out) :: fmax 528 integer , intent(out) :: supnbnivest 529 logical , intent(in) :: est_moy 530 character(len=4) , intent(in) :: option_var,option_anal 531 logical ,intent(in) :: hdam 532 ! 533 534 ! Divers 535 real (kind=dp) ,dimension(:),allocatable :: val,par,borni,borns,par2 536 real (kind=dp) :: user(1) 537 real (kind=dp) :: xlrt_t,dx,f1,var_yQy,weight 538 real (kind=dp) :: xlrt_t21,xlrt_t20,f2,f21 539 540 logical itest,stvecsol(size(vecsol1)) 541 integer :: ip,i,n,ilong,ibound,npar,nparx,ix,j,ifail,iuser(0),chr 542 integer :: ii,imoy,kd,nstart,nend 543 logical :: details,hsire 544 integer :: nb_var,jm 545 546 integer :: i_haplo,ig,kkd 547 548 logical , dimension(:,:,:),allocatable :: filter_inc 549 550 variance_homo = (option_var == 'homo') 551 is_la = (option_anal == 'LA ') 552 is_ld = (option_anal == 'LD ') 553 is_ldla = (option_anal == 'LDLA') 554 is_ldjh = (option_anal == 'LDJH') 555 hsire = ( is_ld .or. is_ldla .or. is_ldjh ) 556 557 ! open(1212,file='lrt_LDJH') 558 nb_var=1 559 if(option_var == 'hete') nb_var=np 560 ! if(option_anal == 'LD ' .or. option_anal == 'LDLA')nb_var=2*nb_var 561 562 nparx = 12*(nb_var+ntnivmax)+((nb_var+ntnivmax)*((nb_var+ntnivmax)-1)/2) 563 allocate ( val(nb_var+ntnivmax ) ) 564 allocate ( par(nb_var + ntnivmax) ) 565 allocate ( par2(nb_var + ntnivmax) ) 566 allocate ( borni(nb_var + ntnivmax) ) 567 allocate ( borns(nb_var + ntnivmax) ) 568 iuser=0 569 570 estime_moy = est_moy 571 572 open(66,file='compALvsJM') 573 !MODIF - OPMIZATION 574 allocate (filter_inc(np,nm, nb_var + ntnivmax)) 575 !FIN MODIF - OPMIZATION 576 577 578 ! 579 !****************************************************************************** 580 ! Calcul de la vraisemblance sous H1 581 582 ! 583 ! initialisation 584 ! on utilisera : 585 ! PAR (vecteur des param�tres � optimiser), 586 ! CORNIV (vecteur des positions, parmi NTNIV, des effets extimables), 587 ! SOLVEC (vecteur disant l'estimabilit� des effets) 588 ! VAL le vecteur complet (ntniv positions) des valeurs des niveaux des effets 589 ! STVAL,STVECSOL et STCORNIV les copies de VAL, VECSOL et CORNIV � l'it�ration 590 ! n-1 quand on analyse la position n 591 ! 592 ! Point de depart 593 ! 594 595 do ip=1,np 596 par(ip)=sig1(ip) 597 end do 598 599 par(np+1)=xmu1g 600 601 do i=np+2,size(par) 602 par(i)=0.d0 603 end do 604 605 do i=np+1,size(val) 606 val(i-np)=par(i) 607 end do 608 ! 609 do i=1,size(vecsol) 610 vecsol(i)=.true. 611 end do 612 613 call new(1,lrtsol) 614 lrtsol%lrtmax(0)=-1.d75 615 lrtsol%nxmax(0)=1 616 lrtsol%chrmax(0)=1 617 lrtsol%lrt1=0.d0 618 619 mod_ic = ic 620 do chr=1,nchr 621 current_chr = chr 622 ! call set_tab_ibs(chr) !peut etre a faire que pour LDJH 623 ! 624 ! Marche le long du chromosome 625 !n=0 626 627 nstart=1 628 do while (absi(chr,nstart)<=posi(chr,longhap/2)) 629 nstart=nstart+1 630 end do 631 632 633 634 nend=get_npo(chr) 635 do while (absi(chr,nend)>= posi(chr,nmk(chr)-longhap/2+1)) 636 nend = nend - 1 637 end do 638 639 ! nend=nstart 640 do n=nstart,nend 641 642 dx=absi(chr,n) 643 ! 644 ! recherche des haplotypes possibles 645 ! 646 if(option_anal == 'LD ' .or. option_anal == 'LDLA' )call set_haplo_for_ldla(chr,dx,n,hsire,hdam) 647 if(option_anal == 'LDJH' )then 648 call liste_shared_haplo(chr,dx,n) 649 do ip = 1,np 650 ! call bilan_shared_haplo(current_chr,dx,n,ip) 651 end do 652 end if 653 ! 654 ! on stocke les conditions en n 655 do i=1,ntniv 656 stvecsol(i)=vecsol(i) 657 end do 658 ! 659 ! preparation de la matrice d'incidence 660 ! 661 call prepinc(current_chr,n,ic,option_anal) 662 ! recherche des elements estimables � partir de la d�composition de choleski 663 ! comptage et recodification de ces elements (on commence par mettre � 0 les 664 ! compteurs destin�s aux tests des effes fix�s 665 ! 666 itest=.false. 667 details=.false. 668 call contingence_ldla(ic,1,itest,chr,n,var_yQy,details,est_moy,option_anal,hsire,hdam) 669 !MODIF - OPMIZATION 670 call set_filter_optim(ic,(option_var == 'hete'),(option_anal == 'LD ' .or. option_anal == 'LDLA'),& 671 ntnivmax,ntniv,vecsol,xinc,filter_inc) 672 !FIN MODIF - OPMIZATION 673 674 675 ! weight=dsqrt(dlog(var_yQy)) 676 ! print*,dx,weight 677 ! go to 111 678 !****************** 679 ! passage des arguments par iuser 680 ! iuser(1) = nb_var 681 ! iuser(2) = n 682 ! iuser(3) = imoy 683 ! iuser(4) = imoyld 684 ! iuser(5) = init 685 ! iuser(6) = jnit 686 mod_nb_var = nb_var 687 mod_pos = n 688 mod_imoy = 0 689 690 if(est_moy) mod_imoy = 1 691 692 mod_imoyld=mod_imoy 693 694 if(option_anal == 'LD ' .or. option_anal == 'LDLA')mod_imoyld=mod_imoy+nb_haplo_reduit 695 696 mod_init=1 ! +1 par rapport au dernier effet estime 697 698 if(mod_imoy==1)mod_init=mod_init+1 ! moyenne 699 700 if(option_anal /= 'LD ')mod_init=mod_init+1 ! qtl pere 701 702 if(option_anal == 'LDJH')mod_init=mod_init+1 ! qtl pere transmis par la mere 703 704 if(nbfem /= 0 .and. (option_anal == 'LA ' .or. option_anal == 'LDLA'))mod_init=mod_init+1 !qtl mere 705 ! if(option_anal == 'LD ' .or. option_anal == 'LDLA')iuser(5)=iuser(5)+nb_max_haplo-1 ! haplotype 706 if(option_anal == 'LD ' .or. option_anal == 'LDLA')mod_init=mod_init+nb_haplo_reduit ! haplotype 707 mod_jnit=mod_init 708 if(nbfem /= 0 )mod_jnit=mod_jnit+1 709 710 711 712 713 ! Parametres de maximisation 714 ibound=0 715 npar=nb_var+nbnivest 716 do ip=1,nb_var 717 borni(ip)=SIG_MIN 718 borns(ip)=SIG_MAX 719 end do 720 do i=nb_var+1,npar 721 borni(i)=XMU_MIN 722 borns(i)=XMU_MAX 723 end do 724 ! 725 ! Point de depart (on reprend les point d'arriv�e pr�c�dents 726 ! 727 do i=nb_var+1,npar 728 par(i)=0.d0 729 end do 730 731 ! j=nb_var 732 ! do i=1,ntniv 733 ! if(vecsol(i))then 734 ! j=j+1 735 ! par(j)=0.d0 736 ! if(stvecsol(i)) par(j)=val(i) 737 ! end if 738 ! end do 739 740 741 write(66,*)' a la position dx = ', dx 742 do ip=1,np 743 do jm=nmp(ip)+1,nmp(ip+1) 744 do ig=ngenom(current_chr,jm)+1,ngenom(current_chr,jm+1) 745 do kd=ngend(current_chr,ig)+1,ngend(current_chr,ig+1) 746 kkd=ndesc(current_chr,kd) 747 if(presentc(mod_ic,kkd)) then 748 write(66,*) n,ip,jm,ig,kd,animal(kkd),(trim(name_haplo_reduit(i_haplo)),& 749 pb_haplo(kkd,i_haplo),i_haplo=1,nb_haplo_reduit) 750 end if 751 end do !kd 752 end do !ig 753 end do !jm 754 end do !ip 755 756 ! Optimisation de la vraisemblance a la position dx pour le modele sans transmission par la mere 757 if(option_anal == 'LDJH') then 758 ifail=1 759 ! iuser(7)=0 760 ! iuser(8)=0 761 call minimizing_funct(npar,ibound,funct_1qtl_modlin_ldla,borni,borns,par,f2,iuser,user,ifail) 762 end if 763 ! Optimisation de la vraisemblance a la position dx 764 ifail=1 765 ! iuser(7)=1 766 ! call minimizing_funct(npar,ibound,funct_1qtl_modlin_ldla,borni,borns,par,f1,iuser,user,ifail) 767 call minimizing_funct_family(npar,ibound,funct_1qtl_modlin_ldla_family,filter_inc(:,:,:npar),& 768 fm1,fp1,borni,borns,par,f1,iuser,user,ifail) 769 770 xlrt_t=f0-f1 771 print*,n,dx,xlrt_t 772 773 ! if(xlrt_t < 0.d0) go to 111 774 775 if(xlrt_t < 0.d0) then 776 print*,'n, nb_haplo_reduit',n, nb_haplo_reduit 777 print *,'par',(par(i),i=nb_var+2,nb_var+1+nb_haplo_reduit) 778 end if 779 780 j=nb_var 781 do i=1,ntniv 782 if(vecsol(i)) then 783 j=j+1 784 val(i)=par(j) 785 else 786 val(i)=9999.d0 787 end if 788 end do 789 ! 790 791 ! 792 ! on met les profil / progeniteur dans les fichier ad hoc 793 ! 794 do ii=1,np 795 lrtsol%xlrp(chr,ii,n)=-1.d0*(fp1(ii)-fp0(ii)) 796 end do 797 do ii=1,nm 798 lrtsol%xlrm(chr,ii,n)=-1.d0*(fm1(ii)-fm0(ii)) 799 end do 800 801 ! on stocke la position si elle est meilleure que les pr�c�dents 802 ! 803 if(lrtsol%lrtmax(0) < xlrt_t) then 804 ! dxmax=dx 805 lrtsol%nxmax(0)=n 806 lrtsol%chrmax(0)=chr 807 lrtsol%lrtmax(0)=xlrt_t 808 fmax=f1 809 supnbnivest=nbnivest 810 do i=1,npar 811 par1(i)=par(i) 812 end do 813 814 end if 815 ! 816 ! on garde les valeurs du LRT pour pouvopir dessiner la courbe de vraisemblance 817 ! 818 lrtsol%lrt1(chr,n)=xlrt_t 819 111 continue 820 end do 821 822 end do 823 ! 824 ! on calcule la pr�cision des estimation au point correspondant 825 ! au LRT maximum 826 ! 827 ! hsire= true pour recuperer les haplotypes des peres 828 call set_haplo_for_ldla(lrtsol%chrmax(0),absi(lrtsol%chrmax(0),lrtsol%nxmax(0)),lrtsol%nxmax(0),.true.,hdam) 829 830 call prepinc(lrtsol%chrmax(0),lrtsol%nxmax(0),ic,option_anal) 831 details=.false. 832 call contingence_ldla(ic,1,itest,lrtsol%chrmax(0),lrtsol%nxmax(0),var_yQy,details,est_moy,option_anal,hsire,hdam) 833 call precision(xx,precis) 834 835 do i=1,ntniv 836 precis1(i)=precis(i) 837 vecsol1(i)=vecsol(i) 838 end do 839 840 deallocate ( val ) 841 deallocate ( par ) 842 deallocate ( borni ) 843 deallocate ( borns ) 844 deallocate (filter_inc) 845 ! close(1212) 846 847 return 848 end subroutine opti_1qtl_modlin_ldla
set_solution_hypothesis0
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
set_solution_hypothesis0
DESCRIPTION
NOTES
SOURCE
1299 subroutine set_solution_hypothesis0(ic,is_hetero,incsol) 1300 integer ,intent(in) :: ic 1301 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1302 logical ,intent(in) :: is_hetero 1303 1304 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 1305 integer :: indest,isol,ipar,nlevel,ief,i,ife,start 1306 1307 real(kind=dp) ,dimension(size(par0)) :: par0_t 1308 1309 ! remise a l echelle des variables 1310 do ipar=1,size(par0) 1311 par0_t(ipar)=par0(ipar)*sigt(ic) 1312 end do 1313 1314 allocate (incsol%sig(1,np)) 1315 1316 incsol%hypothesis=0 1317 ! Mean, Polygenic family 1318 nteff = 2 1319 if ( count(estime(ic,:)) > 0 ) nteff = 3 1320 1321 !Fixed effect and covariate 1322 if ( modele(ic,1) > 0 ) nteff = nteff+1 1323 if ( modele(ic,2) > 0 ) nteff = nteff+1 1324 1325 maxNbPar = max(np,count(estime(ic,:))) 1326 1327 !max numbre de covariable ? 1328 maxNbPar = max(maxNbPar,modele(ic,2)) 1329 1330 nlevel=0 1331 !max nombre de niveau pour un effet fixe ? 1332 do i=1,modele(ic,1) 1333 nlevel=nlev(modele(ic,3+i))+nlevel 1334 end do 1335 maxNbPar = max(maxNbPar,nlevel) 1336 1337 allocate (incsol%groupeName(nteff)) 1338 allocate (incsol%eqtl_print(nteff)) 1339 allocate (incsol%nbParameterGroup(nteff)) 1340 allocate (incsol%parameterName(nteff,maxNbPar)) 1341 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1342 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1343 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 1344 incsol%eqtl_print=.true. 1345 incsol%parameterPrecis=0.d0 1346 incsol%parameterVecsol=.true. 1347 1348 if ( is_hetero ) then 1349 do ip=1,np 1350 incsol%sig(1,ip) = par0_t(ip) 1351 end do 1352 start = np 1353 else 1354 incsol%sig(1,:np) = par0_t(1) 1355 start = 1 1356 end if 1357 1358 ieff=1 1359 incsol%groupeName(ieff) = 'General Mean' 1360 incsol%nbParameterGroup(ieff)=1 1361 incsol%parameterName(ieff,1) ='General Mean' 1362 incsol%paramaterValue(ieff,1) = par0_t(start+1)+xmut(ic) 1363 incsol%parameterVecsol(ieff,1) = vecsol0(1) 1364 incsol%parameterPrecis(ieff,1) = precis0(1) 1365 1366 1367 ieff=ieff+1 1368 incsol%groupeName(ieff) = 'Sire polygenic effects' 1369 incsol%nbParameterGroup(ieff)=np 1370 1371 indest = start+1 1372 isol=1 1373 do ip=1,np 1374 isol=isol+1 1375 if ( vecsol0(isol) ) then 1376 indest=indest+1 1377 incsol%paramaterValue(ieff,ip) = par0_t(indest) 1378 else 1379 incsol%paramaterValue(ieff,ip) = 0.d0 1380 end if 1381 1382 incsol%parameterName(ieff,ip) ='Sire '//trim(pere(ip)) 1383 incsol%parameterVecsol(ieff,ip) = vecsol0(isol) 1384 incsol%parameterPrecis(ieff,ip) = precis0(isol) 1385 end do 1386 1387 if ( count(estime(ic,:)) > 0 ) then 1388 ieff = ieff +1 1389 incsol%groupeName(ieff) = 'Dam polygenic effects' 1390 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 1391 ifem=0 1392 do ip=1,np 1393 do jm=nmp(ip)+1,nmp(ip+1) 1394 if ( estime(ic,jm)) then 1395 isol=isol+1 1396 ifem=ifem+1 1397 if ( vecsol0(isol) ) then 1398 indest=indest+1 1399 incsol%paramaterValue(ieff,ifem) = par0_t(indest) 1400 else 1401 incsol%paramaterValue(ieff,ifem) = 0.d0 1402 end if 1403 1404 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1405 incsol%parameterVecsol(ieff,ifem) = vecsol0(isol) 1406 incsol%parameterPrecis(ieff,ifem) = precis0(isol) 1407 end if 1408 end do 1409 end do 1410 end if 1411 1412 !Fixed effect 1413 if ( modele(ic,1) > 0 ) then 1414 ieff=ieff+1 1415 incsol%eqtl_print(ieff)=.false. 1416 incsol%groupeName(ieff) = 'Fixed effects' 1417 incsol%nbParameterGroup(ieff)=nlevel 1418 ife=0 1419 do ief=1,nbef 1420 do i=1,nlev(modele(ic,3+ief)) 1421 ife=ife+1 1422 isol=isol+1 1423 1424 if ( vecsol0(isol) ) then 1425 indest=indest+1 1426 incsol%paramaterValue(ieff,ife) = par0_t(indest) 1427 else 1428 incsol%paramaterValue(ieff,ife) = 0.d0 1429 end if 1430 incsol%parameterName(ieff,ife) = trim(namefix(modele(ic,3+ief)))//' Level '//trim(str(i)) 1431 incsol%parameterVecsol(ieff,ife) = vecsol0(isol) 1432 incsol%parameterPrecis(ieff,ife) = precis0(isol) 1433 end do 1434 end do 1435 end if 1436 1437 !Covariate 1438 if ( modele(ic,2) > 0 ) then 1439 ieff=ieff+1 1440 incsol%eqtl_print(ieff)=.false. 1441 incsol%groupeName(ieff) = 'Covariates' 1442 incsol%nbParameterGroup(ieff)=modele(ic,2) 1443 1444 do ief=1,modele(ic,2) 1445 isol=isol+1 1446 if ( vecsol0(isol) ) then 1447 indest=indest+1 1448 incsol%paramaterValue(ieff,ief) = par0_t(indest) 1449 else 1450 incsol%paramaterValue(ieff,ief) = 0.d0 1451 end if 1452 incsol%parameterName(ieff,ief) = trim(namecov(modele(ic,3+modele(ic,1)+ief))) 1453 incsol%parameterVecsol(ieff,ief) = vecsol0(isol) 1454 incsol%parameterPrecis(ieff,ief) = precis0(isol) 1455 end do 1456 end if 1457 1458 end subroutine set_solution_hypothesis0
set_solution_hypothesis1
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
set_solution_hypothesis1
DESCRIPTION
NOTES
SOURCE
1470 subroutine set_solution_hypothesis1(ic,is_hetero,qtl,hsire,hdam,incsol) 1471 integer ,intent(in) :: ic 1472 type(TYPE_INCIDENCE_SOLUTION) ,intent(inout) :: incsol 1473 logical ,intent(in) :: is_hetero,qtl,hsire,hdam 1474 1475 integer :: nteff,maxNbPar,jm,ifem,ip,ieff 1476 integer :: ntlev,nbtp,jef,lp,indest,km 1477 integer :: isol,ipar,nlevel,ief,i,ife,start,i_haplo 1478 integer :: jhr,j 1479 1480 real(kind=dp) ,dimension(size(par1)) :: par1_t 1481 character(LEN=400) nhr(2) 1482 1483 ! remise a l echelle des variables 1484 do ipar=1,size(par1) 1485 par1_t(ipar)=par1(ipar)*sigt(ic) 1486 end do 1487 1488 allocate (incsol%sig(1,np)) 1489 if(hsire) allocate (incsol%unknown_dam_sig(1,np)) 1490 1491 incsol%hypothesis=1 1492 ! Mean, Polygenic family 1493 nteff = 2 1494 if ( count(estime(ic,:)) > 0 ) nteff = nteff +1 1495 1496 ! Haplotype sire 1497 if ( hsire ) then 1498 nteff = nteff + 1 1499 end if 1500 1501 ! QTLEffect 1502 if ( qtl ) then 1503 nteff = nteff + 1 1504 if ( count(estime(ic,:)) > 0 ) nteff = nteff + 1 1505 end if 1506 1507 !Fixed effect and covariate 1508 if ( modele(ic,1) > 0 ) nteff = nteff+1 1509 if ( modele(ic,2) > 0 ) nteff = nteff+1 1510 1511 maxNbPar = max(np,count(estime(ic,:))) 1512 !max numbre de covariable ? 1513 maxNbPar = max(maxNbPar,modele(ic,2)) 1514 1515 nlevel=0 1516 !max nombre de niveau pour un effet fixe ? 1517 do i=1,modele(ic,1) 1518 nlevel=nlev(modele(ic,3+i))+nlevel 1519 end do 1520 maxNbPar = max(maxNbPar,nlevel) 1521 1522 ntlev=1 1523 nbtp = 3 + modele(ic,1)+modele(ic,2)+modele(ic,3) 1524 do jef=1,modele(ic,3) 1525 ntlev=ntlev*nlev(modele(ic,nbtp+jef)) 1526 end do 1527 1528 !max nombre de niveau pour un effet fixe en interaction avec le qtl ? 1529 maxNbPar = max(maxNbPar,ntlev*np) 1530 maxNbPar = max(maxNbPar,ntlev*count(estime(ic,:))) 1531 maxNbPar = max(maxNbPar,nb_haplo_reduit) 1532 1533 allocate (incsol%groupeName(nteff)) 1534 allocate (incsol%nbParameterGroup(nteff)) 1535 allocate (incsol%eqtl_print(nteff)) 1536 allocate (incsol%parameterName(nteff,maxNbPar)) 1537 allocate (incsol%paramaterValue(nteff,maxNbPar)) 1538 allocate (incsol%parameterVecsol(nteff,maxNbPar)) 1539 allocate (incsol%parameterPrecis(nteff,maxNbPar)) 1540 allocate (incsol%qtl_groupeName(1,1)) 1541 incsol%parameterPrecis=0.d0 1542 incsol%parameterVecsol=.true. 1543 incsol%eqtl_print=.true. 1544 1545 if ( is_hetero ) then 1546 do ip=1,np 1547 incsol%sig(1,ip) = par1_t(ip) 1548 ! if(hsire) incsol%unknown_dam_sig(1,ip) = par1_t(np+ip) 1549 end do 1550 ! start = 2*np 1551 start=np 1552 else 1553 incsol%sig(1,:np) = par1_t(1) 1554 ! if(hsire) incsol%unknown_dam_sig(1,:np) = par1_t(2) 1555 ! start = 2 1556 start=1 1557 end if 1558 1559 ieff=1 1560 incsol%groupeName(ieff) = 'General Mean' 1561 incsol%nbParameterGroup(ieff)=1 1562 incsol%parameterName(ieff,1) ='General Mean' 1563 incsol%paramaterValue(ieff,1) = par1_t(start+1)+xmut(ic) 1564 incsol%parameterVecsol(ieff,1) = vecsol1(1) 1565 incsol%parameterPrecis(ieff,1) = precis1(1) 1566 1567 indest = start+1 1568 isol=1 1569 1570 if ( hsire ) then 1571 ieff=ieff+1 1572 incsol%qtl_groupeName(1,1)=ieff ! a modifier.. 1573 incsol%groupeName(ieff) = 'Haplotypes effects' 1574 incsol%nbParameterGroup(ieff)=nb_haplo_reduit 1575 1576 do i_haplo=1,nb_haplo_reduit 1577 isol=isol+1 1578 if(vecsol1(isol)) then 1579 indest=indest+1 1580 incsol%paramaterValue(ieff,i_haplo) = par1_t(indest) 1581 else 1582 incsol%paramaterValue(ieff,i_haplo) = 0.d0 1583 end if 1584 incsol%parameterName(ieff,i_haplo)=trim(name_haplo_reduit(i_haplo))//"(race="//trim(race_haplo_reduit(i_haplo)) & 1585 //", freq="//trim(str(pb_haplo_reduit(i_haplo)))//")" 1586 incsol%parameterVecsol(ieff,i_haplo) = vecsol1(isol) 1587 incsol%parameterPrecis(ieff,i_haplo) = precis1(isol) 1588 end do 1589 end if 1590 1591 if ( qtl ) then 1592 ieff=ieff+1 1593 incsol%qtl_groupeName(1,1)=ieff 1594 incsol%groupeName(ieff) = 'Sire QTL effects' 1595 incsol%nbParameterGroup(ieff)=np*ntlevp 1596 1597 ife=0 1598 do ip=1,np 1599 do lp=1,ntlev 1600 ife=ife+1 1601 isol=isol+1 1602 if ( vecsol1(isol) ) then 1603 indest=indest+1 1604 incsol%paramaterValue(ieff,ife) = par1_t(indest) 1605 else 1606 incsol%paramaterValue(ieff,ife) = 0.d0 1607 end if 1608 1609 nhr= ' ' 1610 do jhr=1,2 1611 if(num_haplo_pere(ip,jhr,1) /= 0) nhr(jhr) =name_haplo_reduit(num_haplo_pere(ip,jhr,1)) 1612 end do 1613 1614 1615 ! incsol%parameterName(ieff,ife) ='Sire '//trim(pere(ip))//" "//trim(str(lp))//" - "& 1616 ! //'['//trim(name_haplo_reduit(num_haplo_pere(ip,1,1)))//'/'& 1617 ! //trim(name_haplo_reduit(num_haplo_pere(ip,2,1)))//']' 1618 1619 1620 incsol%parameterName(ieff,ife) ='Sire '//trim(pere(ip))//" "//trim(str(lp))//" - "& 1621 //'['//trim(nhr(1))//'/'//trim(nhr(2))//']' 1622 1623 1624 ! print*,'ip,num_haplo_pere(ip,2,1),name_haplo_reduit(num_haplo_pere(ip,2,1))',& 1625 ! ip,num_haplo_pere(ip,2,1),name_haplo_reduit(num_haplo_pere(ip,2,1)) 1626 1627 ! 1628 ! a voir avec ofi en fait il y a maintenant > 1 phases possibles 1629 ! 1630 incsol%parameterVecsol(ieff,ife) = vecsol1(isol) 1631 incsol%parameterPrecis(ieff,ife) = precis1(isol) 1632 end do 1633 end do 1634 1635 if ( count(estime(ic,:)) > 0 ) then 1636 ieff = ieff +1 1637 incsol%groupeName(ieff) = 'Dam QTL effects' 1638 incsol%nbParameterGroup(ieff)=namest(ic)*ntlevp 1639 ifem=0 1640 do jm=1,nm 1641 if ( estime(ic,jm)) then 1642 do lp=1,ntlev 1643 isol=isol+1 1644 ifem=ifem+1 1645 if ( vecsol1(isol) ) then 1646 indest=indest+1 1647 incsol%paramaterValue(ieff,ifem) = par1_t(indest) 1648 else 1649 incsol%paramaterValue(ieff,ifem) = 0.d0 1650 end if 1651 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" "//trim(str(lp)) 1652 incsol%parameterVecsol(ieff,ifem) = vecsol1(isol) 1653 incsol%parameterPrecis(ieff,ifem) = precis1(isol) 1654 end do 1655 end if 1656 end do 1657 end if 1658 end if 1659 1660 ieff=ieff+1 1661 incsol%groupeName(ieff) = 'Sire polygenic effects' 1662 incsol%nbParameterGroup(ieff)=np 1663 1664 do ip=1,np 1665 isol=isol+1 1666 if ( vecsol1(isol) ) then 1667 indest=indest+1 1668 incsol%paramaterValue(ieff,ip) = par1_t(indest) 1669 else 1670 incsol%paramaterValue(ieff,ip) = 0.d0 1671 end if 1672 1673 incsol%parameterName(ieff,ip) ='Sire '//trim(pere(ip)) 1674 incsol%parameterVecsol(ieff,ip) = vecsol1(isol) 1675 incsol%parameterPrecis(ieff,ip) = precis1(isol) 1676 end do 1677 1678 if ( count(estime(ic,:)) > 0 ) then 1679 ieff = ieff +1 1680 incsol%groupeName(ieff) = 'Dam polygenic effects' 1681 incsol%nbParameterGroup(ieff)=count(estime(ic,:)) 1682 ifem=0 1683 do ip=1,np 1684 do jm=nmp(ip)+1,nmp(ip+1) 1685 if ( estime(ic,jm)) then 1686 isol=isol+1 1687 ifem=ifem+1 1688 if ( vecsol1(isol) ) then 1689 indest=indest+1 1690 incsol%paramaterValue(ieff,ifem) = par1_t(indest) 1691 else 1692 incsol%paramaterValue(ieff,ifem) = 0.d0 1693 end if 1694 1695 incsol%parameterName(ieff,ifem)='Dam '//trim(mere(jm))//" [Sire "//trim(pere(ip))//"]" 1696 incsol%parameterVecsol(ieff,ifem) = vecsol1(isol) 1697 incsol%parameterPrecis(ieff,ifem) = precis1(isol) 1698 end if 1699 end do 1700 end do 1701 end if 1702 1703 !Fixed effect 1704 if ( modele(ic,1) > 0 ) then 1705 ieff=ieff+1 1706 incsol%eqtl_print(ieff)=.false. 1707 incsol%groupeName(ieff) = 'Fixed effects' 1708 incsol%nbParameterGroup(ieff)=nlevel 1709 ife=0 1710 do ief=1,nbef 1711 do i=1,nlev(modele(ic,3+ief)) 1712 isol=isol+1 1713 ife=ife+1 1714 if ( vecsol1(isol) ) then 1715 indest=indest+1 1716 incsol%paramaterValue(ieff,ife) = par1_t(indest) 1717 else 1718 incsol%paramaterValue(ieff,ife) = 0.d0 1719 end if 1720 incsol%parameterName(ieff,ife) = trim(namefix(modele(ic,3+ief)))//' Level '//trim(str(i)) 1721 incsol%parameterVecsol(ieff,ife) = vecsol1(isol) 1722 incsol%parameterPrecis(ieff,ife) = precis1(isol) 1723 end do 1724 end do 1725 end if 1726 1727 !Covariate 1728 if ( modele(ic,2) > 0 ) then 1729 ieff=ieff+1 1730 incsol%eqtl_print(ieff)=.false. 1731 incsol%groupeName(ieff) = 'Covariates' 1732 incsol%nbParameterGroup(ieff)=modele(ic,2) 1733 1734 do ief=1,modele(ic,2) 1735 isol=isol+1 1736 if ( vecsol1(isol) ) then 1737 indest=indest+1 1738 incsol%paramaterValue(ieff,ief) = par1_t(indest) 1739 else 1740 incsol%paramaterValue(ieff,ief) = 0.d0 1741 end if 1742 incsol%parameterName(ieff,ief) = trim(namecov(modele(ic,3+modele(ic,1)+ief))) 1743 incsol%parameterVecsol(ieff,ief) = vecsol1(isol) 1744 incsol%parameterPrecis(ieff,ief) = precis1(isol) 1745 end do 1746 end if 1747 1748 end subroutine set_solution_hypothesis1
test_lin_ldla
[ Top ] [ m_qtlmap_analyse_modlin_ldla ] [ Subroutines ]
NAME
test_lin_ldla
DESCRIPTION
Test des differents effets de nuisance du modele par une LRT compere a une chi2
NOTES
SOURCE
1130 subroutine test_lin_ldla (chr,ic,est_moy,supnbnivest,fmax,nposx,par1,option_var,option_anal) 1131 use m_qtlmap_analyse_lin_gen, only : prepinc,contingence,nbef,nbco,meff,mcov,mint,nbnivest 1132 1133 logical , intent(in) :: est_moy 1134 integer , intent(in) :: chr,ic 1135 integer , intent(in) :: supnbnivest 1136 real (kind=dp) , intent(in) :: fmax 1137 integer , intent(in) :: nposx 1138 real (kind=dp),dimension(3*np+2*nm), intent(in) :: par1 1139 character(len=4) , intent(in) :: option_var,option_anal 1140 1141 ! 1142 ! Divers 1143 logical itest 1144 integer iw((3*np+2*nm) +2),iuser(1) 1145 real (kind=dp) :: par(3*np+2*nm),borni(3*np+2*nm),borns(3*np+2*nm) 1146 real (kind=dp) ::w(12*(3*np+2*nm)+((3*np+2*nm)*((3*np+2*nm)-1)/2)) 1147 real (kind=dp) :: user(1),prob,xlrt_t,f1 1148 integer :: nbint,iecd,iecq,ief,ibound,npar,ip,i,ifail,liw,liwx 1149 integer :: lw,lwx,nbreduit 1150 1151 lwx = (12*(3*np+2*nm))+((3*np+2*nm)*((3*np+2*nm)-1)/2) 1152 liwx = (3*np+2*nm) +2 1153 write(nficout,600) 1154 600 format(//,80('*')/'test of the effets of the model',//) 1155 1156 variance_homo = (option_var == 'homo') 1157 is_la = (option_anal == 'LA ') 1158 is_ld = (option_anal == 'LD ') 1159 is_ldla = (option_anal == 'LDLA') 1160 is_ldjh = (option_anal == 'LDJH') 1161 1162 ! 1163 ! preparation de la matrice d'incidence 1164 ! 1165 call prepinc(chr,nposx,ic,option_anal) 1166 ! 1167 ! on recalcule la vraisemblance en retirant les effets du mod�les un � un 1168 ! 1169 ! 1170 ! on commence par les effets hierarchis�s dans les effets qtl 1171 ! 1172 current_chr = chr 1173 nbef=modele(ic,1) 1174 nbco=modele(ic,2) 1175 nbint=modele(ic,3) 1176 1177 iecd=0 1178 iecq=0 1179 1180 1181 do ief=1,nbint+nbef+nbco 1182 meff=0 1183 mcov=0 1184 mint=0 1185 if(ief.le.nbef)meff=ief 1186 if(ief.gt.nbef.and.ief.le.(nbef+nbco))mcov=ief-nbef 1187 if(ief.gt.(nbef+nbco))mint=ief-(nbef+nbco) 1188 ! 1189 ! recherche des elements estimables � partir de la d�composition de choleski 1190 ! comptage et recodification de ces elements 1191 ! 1192 1193 itest=.true. 1194 call contingence(ic,1,itest,est_moy) 1195 1196 ! 1197 ! la r�dustion du nombre d'effet estim�e est calcul�e 1198 ! 1199 nbreduit=supnbnivest-nbnivest 1200 ! 1201 ! Parametres de maximisation 1202 ibound=0 1203 npar=np+nbnivest 1204 do ip=1,np 1205 borni(ip)=1.do-6 1206 borns(ip)=1.d6 1207 end do 1208 do i=np+1,npar 1209 borni(i)=-1d6 1210 borns(i)=1.d6 1211 par(i)=0.d0 1212 end do 1213 ! 1214 ! Point de depart (on reprend les point d'arriv�e pr�c�dents 1215 ! 1216 do i=1,np 1217 par(i)=par1(i) 1218 end do 1219 1220 do i=np+1,npar 1221 par(i)=0.d0 1222 end do 1223 ! 1224 ! Optimisation de la vraisemblance a la position dx 1225 ifail=1 1226 liw=liwx 1227 lw=lwx 1228 1229 call minimizing_funct(npar,ibound,funct_1qtl_modlin_ldla,borni,borns,par,f1,iuser,user,ifail) 1230 ! 1231 xlrt_t=-2.d0*(fmax-f1) 1232 if(xlrt_t.le.1.do-8)xlrt_t=0.d0 1233 ! 1234 ! impression des tests des effets du mod�le 1235 ! 1236 1237 nbef=modele(ic,1) 1238 nbco=modele(ic,2) 1239 nbint=modele(ic,3) 1240 1241 if(ief.le.(nbef+nbco).and.iecd.eq.0) then 1242 iecd=1 1243 write(nficout,610) 1244 610 format(' Direct effects'/) 1245 1246 write(nficout,601) 1247 601 format('Tested effect df. Likelihood p-value'/ & 1248 ' ratio '/) 1249 end if 1250 1251 if(ief.le.nbef) then 1252 ifail=0 1253 ! prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nlev(modele(ic,3+ief))),ifail) 1254 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nbreduit),ifail) 1255 1256 write(nficout,602) trim(namefix(modele(ic,3+ief))),nbreduit,xlrt_t,prob 1257 ! & nlev(modele(ic,3+ief)),xlrt,prob 1258 end if 1259 602 format(a15,3x,i2,6x,f8.2,7x,f5.3) 1260 1261 if(ief.gt.nbef.and.ief.le.(nbef+nbco)) then 1262 ifail=0 1263 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(1.0),ifail) 1264 1265 write(nficout,603) trim(namecov(modele(ic,3+ief))),xlrt_t,prob 1266 end if 1267 603 format(a15,4x,'1',6x,f8.2,7x,f5.3) 1268 1269 if(ief.gt.(nbef+nbco).and.iecq.eq.0)then 1270 iecq=1 1271 write(nficout,611) 1272 611 format(/' Intra qtl effects'/) 1273 write(nficout,601) 1274 end if 1275 1276 if(ief.gt.(nbef+nbco)) then 1277 ifail=0 1278 prob=MATH_QTLMAP_G01ECF('U',xlrt_t,dble(nbreduit),ifail) 1279 write(nficout,602) trim(namefix(modele(ic,3+ief))),nbreduit,xlrt_t,prob 1280 1281 end if 1282 1283 end do 1284 1285 return 1286 end subroutine test_lin_ldla