TABLE OF CONTENTS
fortran/fmkfrail [ Methods ]
NAME
fmkfrail --- compute BLUP frailty estimates
FUNCTION
Compute cluster- and subject-level frailty BLUPs under the auxiliary Poisson model. See also the R function fupdatefrailties3.
SYNOPSIS
490 subroutine fmkfrail(index,indexd,delta,deltad,Z,Zd,alphars, 491 $ alpharsd,as,asd,betahat,betadhat,m,Ji,jimax,Jicum,jisum,d1, 492 $ d2,ncovs1,ncovs2,nr,ns,nsd,sigma2,sigma2d,nu2,nu2d,theta,Uihat, 493 $ Vihat,Uijhat,Vijhat,pi,qi,ri,si,piprime,qiprime,riprime,siprime, 494 $ pij,qij,rij,sij,pijprime,qijprime,rijprime,sijprime,wi,wij,zij)
INPUTS
Too many to list, but analogous to those in other routines in this file. Inputs appended with d refer to process 2.
OUTPUTS
Uihat vector of cluster-level frailties for process 1
Vihat vector of cluster-level frailties for process 2
Uijhat vector of subject-level frailties for process 1
Vijhat vector of subject-level frailties for process 2
p,q,r,s, etc intermediate values
SOURCE
496 integer d1,d2,ncovs1,ncovs2,nr,ns,nsd,m,jimax 497 integer Ji(m),Jicum(m) 498 integer index(d1,6), indexd(d2,6) 499 double precision delta(d1), deltad(d2) 500 double precision Z(d1,ncovs1), Zd(d2,ncovs2) 501 double precision alphars(nr,ns),alpharsd(nr,nsd) 502 double precision as(nr,ns), asd(nr,nsd) 503 double precision betahat(ncovs1), betadhat(ncovs2) 504 double precision sigma2,sigma2d,nu2,nu2d,theta 505 double precision Uihat(m),Vihat(m),Uijhat(jisum),Vijhat(jisum) 506 double precision pi(m),qi(m),ri(m),si(m),piprime(m),qiprime(m), 507 $ riprime(m),siprime(m),pij(m,jimax),qij(m,jimax),rij(m,jimax), 508 $ sij(m,jimax),pijprime(m,jimax),qijprime(m,jimax), 509 $ rijprime(m,jimax),sijprime(m,jimax) 510 double precision wi(m),wij(m,jimax),zij(m,jimax) 511 512 integer i,j 513 double precision Smu(m,jimax),Seta(m,jimax),Sdelta(m,jimax), 514 $ SDeltad(m,jimax) 515 double precision Smuij,Setaij,Sdeltaij,Sdeltaijd 516 517 call dzero(Smu,m*jimax) 518 call dzero(Seta,m*jimax) 519 call dzero(Sdelta,m*jimax) 520 call dzero(SDeltad,m*jimax) 521 522 ! call dblepr("Smu",-1,Smu,m*jimax) 523 524 call fsmuij(betahat,index,delta,Z,alphars,as,d1,ncovs1,nr,ns, 525 $ m,jimax,Smu,Sdelta) 526 call fsmuij(betadhat,indexd,deltad,Zd,alpharsd,asd,d2,ncovs2,nr, 527 $ nsd,m,jimax,Seta,SDeltad) 528 529 do 100 i=1,m 530 ! compute intermediate values 531 do 110 j=1,Ji(i) 532 Smuij=Smu(i,j) 533 Setaij=Seta(i,j) 534 Sdeltaij=Sdelta(i,j) 535 SDeltaijd=SDeltad(i,j) 536 537 wij(i,j)=nu2-theta*theta*Setaij/(1.0d0+nu2d*Setaij) 538 zij(i,j)=nu2d-theta*theta*Smuij/(1.0d0+nu2*Smuij) 539 pij(i,j)=Smuij/(1.0d0+wij(i,j)*Smuij) 540 qij(i,j)=-theta*Smuij*Setaij/((1.0d0+nu2*Smuij)* 541 $ (1.0d0+nu2d*Setaij)-theta*theta*Smuij*Setaij) 542 rij(i,j)=qij(i,j) 543 sij(i,j)=Setaij/(1.0d0+zij(i,j)*Setaij) 544 pijprime(i,j)=Sdeltaij/(1.0d0+wij(i,j)*Smuij) 545 qijprime(i,j)=-theta*Smuij*SDeltaijd/((1.0d0+nu2*Smuij)* 546 $ (1.0d0+nu2d*Setaij)-theta*theta*Smuij*Setaij) 547 rijprime(i,j)=-theta*Sdeltaij*Setaij/((1.0d0+nu2*Smuij)* 548 $ (1.0d0+nu2d*Setaij)-theta*theta*Smuij*Setaij) 549 sijprime(i,j)=SDeltaijd/(1.0d0+zij(i,j)*Setaij) 550 551 piprime(i)=piprime(i)+pijprime(i,j) 552 qiprime(i)=qiprime(i)+qijprime(i,j) 553 riprime(i)=riprime(i)+rijprime(i,j) 554 siprime(i)=siprime(i)+sijprime(i,j) 555 pi(i)=pi(i)+pij(i,j) 556 qi(i)=qi(i)+qij(i,j) 557 ri(i)=ri(i)+rij(i,j) 558 si(i)=si(i)+sij(i,j) 559 110 continue 560 ! compute cluster-level frailties 561 wi(i)=1.0d0/((1.0d0+sigma2*pi(i))*(1.0d0+sigma2d*si(i)) 562 $ -sigma2*sigma2d*ri(i)*ri(i)) 563 Uihat(i)=1.0d0+sigma2*wi(i)*((1.0d0+sigma2d*si(i))* 564 $ (piprime(i)-pi(i)+qiprime(i)-qi(i))-sigma2d*ri(i)* 565 $ (riprime(i)-ri(i)+siprime(i)-si(i))) 566 Vihat(i)=1.0d0+sigma2d*wi(i)*((1.0d0+sigma2*pi(i))* 567 $ (riprime(i)-ri(i)+siprime(i)-si(i))-sigma2*ri(i)* 568 $ (piprime(i)-pi(i)+qiprime(i)-qi(i))) 569 570 if(Uihat(i).LE.0.01) then 571 Uihat(i)=0.01 572 endif 573 if(Vihat(i).LE.0.01) then 574 Vihat(i)=0.01 575 endif 576 ! compute subject-level frailties 577 do 120 j=1,Ji(i) 578 Uijhat(Jicum(i)+j)=Uihat(i)-(nu2*pij(i,j)+theta*rij(i,j)) 579 $ *Uihat(i)-(nu2*qij(i,j)+theta*sij(i,j))*Vihat(i) 580 $ +nu2*(pijprime(i,j)+qijprime(i,j)) 581 $ +theta*(rijprime(i,j)+sijprime(i,j)) 582 Vijhat(Jicum(i)+j)=Vihat(i)-(theta*qij(i,j)+nu2d*sij(i,j)) 583 $ *Vihat(i)-(theta*pij(i,j)+nu2d*rij(i,j))*Uihat(i) 584 $ +theta*(pijprime(i,j)+qijprime(i,j)) 585 $ +nu2d*(rijprime(i,j)+sijprime(i,j)) 586 if(Uijhat(Jicum(i)+j).LE.0.01) then 587 Uijhat(Jicum(i)+j)=0.01 588 endif 589 if(Vijhat(Jicum(i)+j).LE.0.01) then 590 Vijhat(Jicum(i)+j)=0.01 591 endif 592 120 continue 593 100 continue 594 595 end