m_qtlmap_incidence_linear
[ Top ] [ ANALYSE ] [ Modules ]
NAME
m_qtlmap_incidence_linear
SYNOPSIS
Model Homoscedastic and Heteroscedastic Linear analysis
Implementation of a likelihood ratio test linearised for Homoscedastic and Heteroscedastic Linear
DESCRIPTION
i : sire index
j : Dam index
k : Progeny index
N i : Number of progeny of Sire i
~
SIG i : Sqrt Variance of Sire family i under H0
^
SIG i : Sqrt Variance of Sire family i under H1 at a position dx
General Case (Heteroscedastic model) :
-----------------------------------
| B = (X'.V-1.X)-1 . X' . V-1.Y
System | ^
| SIG² i = SUM [ Y ijk - X ijk . B]² CD² ijk / Ni
becomes in homoscedastic model :
------------------------------
| B = (X'.X)-1 . X'.Y
System | ^
| SIG² = SUM [ Y ijk - X ijk . B]² CD² ijk / Ni
~ ^ ~ ^
LRT = -2*Log(VL) + 2*Log(VL) = SUM i[ Ni * (Log(SIG²) - Log(SIG²)) ]
NOTES
SEE ALSO
model_lin_h0
[ Top ] [ m_qtlmap_incidence_linear ] [ Subroutines ]
NAME
model_lin_h0
DESCRIPTION
compute under H0 the LRT and the solution
INPUTS
xinc : incidence matrix incidenceDesc : description of incidence matrix (see INCIDENCE_TYPE) workstruct : configuration and information about analysis (see INCIDENCE_GEN_STRUCT) performPrecision : boolean to compute for each parameter the precision (fill the variable embeded INCIDENCE_TYPE%precis)
OUTPUTS
sigsquareEstime : the variance of each sire family
Bestim : the solution of the estimation
NOTES
SOURCE
257 subroutine model_lin_h0(xinc,incidenceDesc,workstruct,sigsquareEstime,Bestim,performPrecision) 258 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 259 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xinc 260 type(INCIDENCE_GEN_STRUCT) ,intent(in) :: workstruct 261 real (kind=dp) , dimension(np) ,intent(out) :: sigsquareEstime 262 real (kind=dp) , dimension(incidenceDesc%ntnivmax) ,intent(out) :: Bestim 263 logical ,intent(in) :: performPrecision!,tConf 264 ! real (kind=dp),dimension(ntnivmax,ntnivmax) ,intent(out) :: tempForConfusion 265 266 real (kind=dp),dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) :: temp 267 real (kind=dp) :: startsig(np) 268 integer :: hypothesis,ip 269 270 !print *,ic 271 !call debug_write_incidence(xinc,incidenceDesc) 272 select case (workstruct%type_model) 273 case (LINEAR_HOMOSCEDASTIC) 274 call modele_homoscedastic(xinc,incidenceDesc,sigsquareEstime,Bestim,performPrecision,temp) 275 ! print *,"HOMO:",osigsquare 276 case (LINEAR_HETEROSCEDASTIC) 277 do ip=1,np 278 startsig(ip) = incidenceDesc%dataset%lSires(ip)%sig0(1) 279 end do 280 call modele_heteroscedastic(xinc,incidenceDesc,startsig,sigsquareEstime,Bestim,performPrecision,temp) 281 ! print *,"HETERO:",osigsquare 282 case default 283 call stop_application("Devel Error : Unknown linear type : "//trim(str(workstruct%type_model))) 284 end select 285 ! print *,"H0:",sigsquareEstime 286 287 end subroutine model_lin_h0
model_lin_hn
[ Top ] [ m_qtlmap_incidence_linear ] [ Subroutines ]
NAME
model_lin_hn
DESCRIPTION
compute under Hn the LRT and the solution
INPUTS
xinc : incidence matrix incidenceDesc : description of incidence matrix (see INCIDENCE_TYPE) curPos : current point to test (see POSITION_LRT_INCIDENCE) workstruct : configuration and information about analysis (see INCIDENCE_GEN_STRUCT) performPrecision : boolean to compute for each parameter the precision (fill the variable embeded INCIDENCE_TYPE%precis) tConf : boolean to get an array used by the confusion method invlrt : boolean to compute the LRT
OUTPUTS
sigsquareEstime : the variance of each sire family
Bestim : the solution of the estimation
tempForConfusion : the buffer array for the confusion function
NOTES
SOURCE
310 subroutine model_lin_hn(xinc,incidenceDesc,curPos,workstruct,sigsquareEstime,Bestim,& 311 performPrecision,tConf,tempForConfusion,invlrt) 312 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 313 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xinc 314 type(POSITION_LRT_INCIDENCE) ,intent(inout) :: curPos 315 type(INCIDENCE_GEN_STRUCT) ,intent(in) :: workstruct 316 real (kind=dp) , dimension(incidenceDesc%ntnivmax) ,intent(out) :: Bestim 317 real (kind=dp) , dimension(np) ,intent(out) :: sigsquareEstime 318 logical ,intent(in) :: performPrecision,tConf 319 real (kind=dp),dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(out) :: tempForConfusion 320 logical ,intent(in) :: invlrt 321 322 real (kind=dp),dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) :: temp 323 real (kind=dp) :: startsig(np) 324 integer :: hypothesis,ip,kd1,kd2 325 326 !print *,ic 327 !call debug_write_incidence(xinc,incidenceDesc) 328 329 select case (workstruct%type_model) 330 case (LINEAR_HOMOSCEDASTIC) 331 call modele_homoscedastic(xinc,incidenceDesc,sigsquareEstime,Bestim,performPrecision,temp) 332 case (LINEAR_HETEROSCEDASTIC) 333 startsig=sqrt(workstruct%sigsquare(workstruct%nqtl,:,1)) 334 call modele_heteroscedastic(xinc,incidenceDesc,startsig,sigsquareEstime,Bestim,performPrecision,temp) 335 case default 336 call stop_application("Devel Error : Unknown linear type : "//trim(str(workstruct%type_model))) 337 end select 338 339 curPos%lrt=0.d0 340 341 !compute LRT 342 343 if (invLrt) then 344 do ip=1,np 345 kd1=incidenceDesc%dataset%lSires(ip)%half_sib%firstKd 346 kd2=incidenceDesc%dataset%lSires(ip)%half_sib%lastkd 347 if (kd2>kd1) then 348 do hypothesis=1,workstruct%nqtl 349 curPos%lrt(hypothesis)=curPos%lrt(hypothesis) + & 350 (kd2-kd1+1)*(log(sigsquareEstime(ip))-log(workstruct%sigsquare(hypothesis,ip,1))) ! (IQ-1) / NQTL QTL 351 end do 352 end if 353 end do 354 else 355 do ip=1,np 356 kd1=incidenceDesc%dataset%lSires(ip)%half_sib%firstKd 357 kd2=incidenceDesc%dataset%lSires(ip)%half_sib%lastkd 358 if (kd2>kd1) then 359 do hypothesis=1,workstruct%nqtl 360 361 curPos%lrt(hypothesis)=curPos%lrt(hypothesis) + & 362 (kd2-kd1+1)*(log(workstruct%sigsquare(hypothesis,ip,1)) - log(sigsquareEstime(ip))) ! (IQ-1) / NQTL QTL 363 end do 364 end if 365 end do 366 end if 367 368 if (tConf) then 369 tempForConfusion = temp 370 end if 371 end subroutine model_lin_hn
modele_heteroscedastic
[ Top ] [ m_qtlmap_incidence_linear ] [ Subroutines ]
NAME
modele_heteroscedastic
DESCRIPTION
From a incidence matrix, get the B estimation and the variance of each sires in a heteroscedastic model
Bs = (X'.V-1.X)-1 . X'.Vs -1 .Y
SIGs = SUM CD * [ Y - X . Bs]² / N
i jk ijk ijk ijk i
NOTES
SOURCE
156 subroutine modele_heteroscedastic(xinc,incidenceDesc,startsig,osigsquare,Bestim,performPrecision,tempForConfusion) 157 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 158 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xinc 159 real (kind=dp) , dimension(np) ,intent(in) :: startsig 160 real (kind=dp) , dimension(np) ,intent(out) :: osigsquare 161 real (kind=dp) , dimension(incidenceDesc%ntnivmax) ,intent(out) :: Bestim 162 logical ,intent(in) :: performPrecision 163 real (kind=dp),dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(out) :: tempForConfusion 164 165 !local 166 real (kind=dp) , dimension(incidenceDesc%ntnivmax) :: RHS,YTemp 167 real (kind=dp) ,dimension(incidenceDesc%dataset%nkd,incidenceDesc%dataset%nkd) :: Vinv 168 real (kind=dp) ,dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) :: XVX 169 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) :: xincreduit 170 real (kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) :: XX,Mr 171 real (kind=dp) , dimension(incidenceDesc%ntniv,incidenceDesc%ntniv) :: triang 172 integer :: ip,ERROR,it,i,j,kd1,kd2 173 real (kind=dp) :: f,XB(incidenceDesc%dataset%nkd),lrtv,diff,lrt,oldlrt,sig(np) 174 175 RHS=0.d0 176 sig=startsig 177 diff = 2*EPS_LINEAR_HETEROSCEDASTIC 178 lrt=INIFINY_REAL_VALUE 179 it=0 180 181 do while ( diff > EPS_LINEAR_HETEROSCEDASTIC ) 182 if ( it < MAX_LINEAR_ITERATION) then 183 it=it+1 184 call set_VInv( incidenceDesc, sig,Vinv ) 185 call model_XT_V_X(xinc,incidenceDesc,Vinv,XVX) 186 187 call estim_cholesky(XVX,incidenceDesc,incidenceDesc%ntniv,triang) 188 call set_corrxinc(xinc,incidenceDesc,xincreduit) 189 190 call set_RHS_V(xincreduit,Vinv,incidenceDesc,RHS) 191 192 do i=1,incidenceDesc%nbnivest 193 YTemp(i) = RHS(i) 194 do j=i-1,1,-1 195 YTemp(i) = YTemp(i) - YTemp(j)*triang(i,j) 196 end do 197 YTemp(i) = YTemp(i) / triang(i,i) 198 end do 199 200 do i=incidenceDesc%nbnivest,1,-1 201 Bestim(i) = YTemp(i) 202 do j=i+1,incidenceDesc%nbnivest 203 Bestim(i) = Bestim(i) - Bestim(j)*triang(j,i) 204 end do 205 Bestim(i) = Bestim(i) / triang(i,i) 206 end do 207 208 XB = matmul(xincreduit(:incidenceDesc%dataset%nkd,:incidenceDesc%nbnivest),Bestim(:incidenceDesc%nbnivest)) 209 XB = incidenceDesc%dataset%Y(1,:)-XB 210 XB = XB*XB*incidenceDesc%dataset%CD(1,:)*incidenceDesc%dataset%CD(1,:) 211 212 do ip=1,np 213 kd1=incidenceDesc%dataset%lSires(ip)%half_sib%firstKd 214 kd2=incidenceDesc%dataset%lSires(ip)%half_sib%lastkd 215 osigsquare(ip)=sum(XB(kd1:kd2)) / (kd2-kd1+1) 216 end do 217 218 oldlrt=lrt 219 lrt = 0 220 do ip=1,np 221 kd1=incidenceDesc%dataset%lSires(ip)%half_sib%firstKd 222 kd2=incidenceDesc%dataset%lSires(ip)%half_sib%lastkd 223 lrt = lrt +(kd2-kd1+1)*(log(osigsquare(ip)) - log(sig(ip)*sig(ip))) 224 sig(ip)=sqrt(osigsquare(ip)) 225 end do 226 diff = abs(lrt-oldlrt) 227 else 228 diff = EPS_LINEAR_HETEROSCEDASTIC 229 end if 230 end do 231 232 233 if (performPrecision) then 234 call get_precision(XVX,tempForConfusion,incidenceDesc) 235 end if 236 237 end subroutine modele_heteroscedastic
modele_homoscedastic
[ Top ] [ m_qtlmap_incidence_linear ] [ Subroutines ]
NAME
modele_homoscedastic
DESCRIPTION
From a incidence matrix, get the B estimation and the variance of each sires in a homoscedastic model
Bs = (X'.X)-1 . X'.Y
SIGs = SUM [ Y - X . Bs]² / N
i jk ijk ijk i
NOTES
SOURCE
83 subroutine modele_homoscedastic(xinc,incidenceDesc,osigsquare,Bestim,performPrecision,tempForConfusion) 84 type(INCIDENCE_TYPE) , intent(inout) :: incidenceDesc 85 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax), intent(in) :: xinc 86 real (kind=dp) , dimension(np) ,intent(out) :: osigsquare 87 real (kind=dp) , dimension(incidenceDesc%ntnivmax) ,intent(out) :: Bestim 88 logical ,intent(in) :: performPrecision 89 real (kind=dp),dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) ,intent(out) :: tempForConfusion 90 91 !local 92 real (kind=dp) , dimension(incidenceDesc%ntnivmax) :: RHS,YTemp 93 real (kind=dp) , dimension(nd,incidenceDesc%ntnivmax) :: xincreduit 94 real (kind=dp) , dimension(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax) :: XX!,Mr 95 real (kind=dp) , dimension(incidenceDesc%ntniv,incidenceDesc%ntniv) :: triang 96 integer :: ip,kd1,kd2,jm,i,ERROR,jjm,j 97 real (kind=dp) :: f,XB(incidenceDesc%dataset%nkd),sigsq 98 99 ! create X'.X matrix from incidence matrix 100 call model_XT_X(xinc,incidenceDesc,XX) 101 102 ! Check all parameters to remove from the estimation 103 call estim_cholesky(XX,incidenceDesc,incidenceDesc%ntniv,triang) 104 ! compute the precision of each parameter 105 if (performPrecision) call get_precision(XX,tempForConfusion,incidenceDesc) 106 ! call debug_write_incidence(xinc,incidenceDesc) 107 108 call set_corrxinc(xinc,incidenceDesc,xincreduit) 109 call set_RHS(xincreduit,incidenceDesc,RHS) 110 111 ! print *,"RHS****" 112 ! print *,RHS(:incidenceDesc%nbnivest) 113 114 do i=1,incidenceDesc%nbnivest 115 YTemp(i) = RHS(i) 116 do j=i-1,1,-1 117 YTemp(i) = YTemp(i) - YTemp(j)*triang(i,j) 118 end do 119 YTemp(i) = YTemp(i) / triang(i,i) 120 ! print *,'Temp(',i,')=',YTemp(i) 121 end do 122 123 do i=incidenceDesc%nbnivest,1,-1 124 Bestim(i) = YTemp(i) 125 do j=i+1,incidenceDesc%nbnivest 126 Bestim(i) = Bestim(i) - Bestim(j)*triang(j,i) 127 end do 128 Bestim(i) = Bestim(i) / triang(i,i) 129 end do 130 131 XB = matmul(xincreduit(:incidenceDesc%dataset%nkd,:incidenceDesc%nbnivest),Bestim(:incidenceDesc%nbnivest)) 132 XB = incidenceDesc%dataset%Y(1,:)-XB 133 XB = XB*XB 134 sigsq=sum(XB) 135 136 ! do i=1,10 137 ! print *,XB(i),incidenceDesc%dataset%Y(1,i) 138 ! end do 139 sigsq = sigsq / incidenceDesc%dataset%nkd 140 osigsquare = sigsq 141 142 end subroutine modele_homoscedastic
opti_0qtl_linear
[ Top ] [ m_qtlmap_incidence_linear ] [ Subroutines ]
NAME
opti_0qtl_linear
DESCRIPTION
<DEV>pour l interaction...voir opti_0qtl
NOTES
SOURCE
382 subroutine opti_0qtl_linear( ic, & ! the trait 383 incsol , & ! 384 sigsquareEstime, & 385 typeModel, & ! type of model homoscedastic or heteroscedastic 386 maxqtl) ! maximum qtl to find 387 388 integer ,intent(in) :: ic,typeModel,maxqtl 389 type(TYPE_INCIDENCE_SOLUTION) ,intent(out) :: incsol 390 real (kind=dp) ,dimension(np),intent(out) :: sigsquareEstime 391 392 integer :: chr,ip,i,hypothesis 393 394 real (kind=dp) ,dimension(:,:),pointer :: xinc 395 real (kind=dp) , dimension(:) ,pointer :: Bestim 396 type(INCIDENCE_TYPE) :: incidenceDesc 397 character(len=500) :: additional_info 398 real (kind=dp) ,dimension(0:maxqtl-1) :: listDx 399 integer ,dimension(0:maxqtl-1) :: listN 400 real (kind=dp) ,dimension(0:maxqtl-1) :: dxmax 401 integer ,dimension(1) :: listnteff 402 type(INCIDENCE_GEN_STRUCT) :: workstruct 403 404 workstruct%type_incidence=INCIDENCE_TYPE_CLASSIC 405 406 call init_workstruct(workstruct,INCIDENCE_TYPE_CLASSIC,LINEAR_HOMOSCEDASTIC,0,.true.,.false.,.false.) 407 !initialisation of incidence matrix 408 call init_incidence(ic,0,incidenceDesc,workstruct%type_incidence) 409 allocate (xinc(nd,incidenceDesc%ntnivmax)) 410 xinc=0.d0 411 allocate (Bestim(incidenceDesc%ntnivmax)) 412 Bestim=0.d0 413 !add general mean to estim 414 call add_general_mean(xinc,incidenceDesc) 415 !add polygenic effect 416 call add_polygenic_effect(xinc,incidenceDesc) 417 !add fixed effect and covariate 418 call add_effcov(xinc,incidenceDesc,0,0) 419 !call debug_write_incidence(xinc,incidenceDesc) 420 !call model analysis 421 sigsquareEstime=0.d0 422 call model_lin_h0(xinc,incidenceDesc,workstruct,sigsquareEstime,Bestim,.true.) 423 424 call set_solution(0,workstruct,sigsquareEstime,Bestim,incidenceDesc,incsol,1,listnteff) 425 !PRINT *,sigsquareEstime 426 427 call end_incidence(incidenceDesc) 428 deallocate (xinc) 429 deallocate (Bestim) 430 call release_ws(workstruct) 431 432 end subroutine opti_0qtl_linear
opti_2qtl_linear_interaction
[ Top ] [ m_qtlmap_incidence_linear ] [ Subroutines ]
NAME
opti_2qtl_linear_interaction
DESCRIPTION
<DEV>
NOTES
SOURCE
444 subroutine opti_2qtl_linear_interaction(ic,workstruct,incsol,lrtsol,type) 445 446 integer , intent(in) :: ic,type 447 type(INCIDENCE_GEN_STRUCT) , intent(in) :: workstruct 448 type(TYPE_INCIDENCE_SOLUTION) ,intent(out) :: incsol 449 type(TYPE_LRT_SOLUTION) ,intent(out) :: lrtsol 450 451 integer :: n,init_n1 452 integer :: n1,ii,iq,npar2,ip,chr,chr2,nteff1,nteff2,nteffinter,listnteff(2) 453 real (kind=dp) :: dx,dx1,f,lrt,sigsquare_t(np) 454 455 real (kind=dp) ,dimension(:,:),pointer :: xinc,tempConfusion 456 real (kind=dp) , dimension(:) ,pointer :: Bestim 457 real (kind=dp) , dimension(np) :: sigsquareEstime 458 type(INCIDENCE_TYPE) :: incidenceDesc 459 460 type(POSITION_LRT_INCIDENCE) :: curPos 461 462 call new(2,lrtsol) 463 464 !initialisation of incidence matrix 465 call init_incidence(ic,2,incidenceDesc,workstruct%type_incidence) 466 allocate (xinc(nd,incidenceDesc%ntnivmax)) 467 xinc=0.d0 468 !Position Allocation 469 allocate (curPos%listN(2)) 470 allocate (curPos%listChr(2)) 471 allocate (curPos%lrt(2)) 472 473 lrtsol%lrtmax=-INIFINY_REAL_VALUE 474 475 allocate (Bestim(incidenceDesc%ntnivmax)) 476 allocate (tempConfusion(incidenceDesc%ntnivmax,incidenceDesc%ntnivmax)) 477 Bestim=0.d0 478 chr=1;n=1;n1=n+1 479 !add general mean to estim 480 call add_general_mean(xinc,incidenceDesc) 481 nteff1=2 482 !add qtl effect/interaction at position n to estim 483 call add_qtleffect(nteff1,n,chr,n,xinc,incidenceDesc,0,.true.) 484 nteff2=incidenceDesc%nteff+1 485 !add qtl effect/interaction at position n+1 to estim 486 call add_qtleffect(nteff2,n1,chr,n+1,xinc,incidenceDesc,0,.true.) 487 488 !add polygenic effect 489 call add_polygenic_effect(xinc,incidenceDesc) 490 !add fixed effect and covariate 491 call add_effcov(xinc,incidenceDesc,0,0) 492 493 !add interaction 494 nteffinter=incidenceDesc%nteff+1 495 call add_qtlinteraction(nteffinter,chr,chr,n,n1,xinc,incidenceDesc) 496 !call debug_write_incidence(xinc,incidenceDesc) 497 498 !Initalisation of the maximum finded 499 500 do chr=1,nchr 501 curPos%listChr(1)=chr 502 do n=1,get_npo(chr) 503 curPos%listN(1)=n 504 !change qtl effect/interaction at position on the 2nd effect (the third is qtl dam) to estim 505 call change_qtleffect(nteff1,1,chr,n,xinc,incidenceDesc,0,.true.) 506 do chr2=chr,nchr 507 curPos%listChr(2)=chr2 508 init_n1=n+1 509 if ( chr2 /= chr) then 510 init_n1=1 511 end if 512 do n1=init_n1,get_npo(chr2) 513 curPos%listN(2)=n1 514 !change qtl effect/interaction at the nteff2 th position to estim 515 call change_qtleffect(nteff2,2,chr2,n1,xinc,incidenceDesc,0,.true.) 516 call change_interaction_effect(nteffinter,chr,chr2,n,n1,xinc,incidenceDesc) 517 518 ! call debug_write_incidence(xinc,incidenceDesc) 519 ! print *,"pos:",n 520 !call model 521 call model_lin_hn(xinc,incidenceDesc,curPos,workstruct,sigsquareEstime,Bestim,.false.,.false.,tempConfusion,.false.) 522 ! !compute LRT 523 ! lrt = 0.d0 524 ! !print *,sigsquare_t 525 ! do ip=1,np 526 ! lrt=lrt + effp(ip)*(log(sigsquare(ip)) - log(sigsquare_t(ip))) 527 ! end do 528 ! print *,n,n1,"LRT :",curPos%lrt 529 ! stop 530 531 lrtsol%lrt1_2(chr,chr2,n,n1)=lrt 532 lrtsol%lrt0_2(chr,chr2,n,n1)=lrt 533 534 ! print *,absi(chr,n),absi(chr2,n1),lrt 535 !print *,incidenceDesc%vecsol 536 537 !print *,incidenceDesc%vecsol(:incidenceDesc%ntniv) 538 ! stop 539 540 541 if ( lrtsol%lrtmax(0) < lrt ) then 542 lrtsol%nxmax(0)=n 543 lrtsol%nxmax(1)=n1 544 lrtsol%chrmax(0)=chr 545 lrtsol%chrmax(1)=chr2 546 lrtsol%lrtmax=lrt 547 end if 548 end do 549 end do 550 end do 551 end do 552 553 call change_qtleffect(nteff1,1,lrtsol%chrmax(0),lrtsol%nxmax(0),xinc,incidenceDesc,0,.true.) 554 call change_qtleffect(nteff2,2,lrtsol%chrmax(1),lrtsol%nxmax(1),xinc,incidenceDesc,0,.true.) 555 call change_interaction_effect(nteffinter,lrtsol%chrmax(0),lrtsol%chrmax(1),lrtsol%nxmax(0),lrtsol%nxmax(1),& 556 xinc,incidenceDesc) 557 558 call model_lin_hn(xinc,incidenceDesc,curPos,workstruct,sigsquareEstime,Bestim,.false.,.false.,tempConfusion,.false.) 559 ! call debug_write_incidence(xinc,incidenceDesc) 560 561 ! print *,"Bestime:" 562 563 564 listnteff(1)=nteff1 565 listnteff(2)=nteff2 566 call set_solution(2,workstruct,sigsquareEstime,Bestim,incidenceDesc,incsol,2,listnteff) 567 568 call end_incidence(incidenceDesc) 569 deallocate (xinc) 570 deallocate(Bestim) 571 deallocate (tempConfusion) 572 deallocate (curPos%listN) 573 deallocate (curPos%listChr) 574 deallocate (curPos%lrt) 575 576 end subroutine opti_2qtl_linear_interaction