TABLE OF CONTENTS


fortran/fmksens2 [ Methods ]

NAME

  fmksens2 --- compute sensitivity matrix

FUNCTION

Construct the sensitivity matrix at the end of the procedure. This is only for the regression coefficients. The full sensitivity is computed by fmksens2full See also the R function makeSens2.

SYNOPSIS

615       subroutine fmksens2(Smat,index1,index2,Z,Zd,alphars,alpharsd,
616      $ as,asd,betahat,betadhat,times1,times2,pi,qi,ri,si,wi,wij,zij,
617      $ sig2,sig2d,nu2,nu2d,theta,ncovs1,ncovs2,nr,ns,nsd,
618      $ d1,d2,m,Ji,maxj)

INPUTS

   Too many to list, but analogous to those in other routines in
   this file. Inputs appended with d refer to process 2.

OUTPUTS

   Smat    matrix of size ncovs1+ncovs2 x ncovs2+ncovs2
           sensitivity matrix

SOURCE

621 ! output of frailty estimation
622       double precision pi(m),qi(m),ri(m),si(m),wi(m),
623      $  wij(m,maxj),zij(m,maxj)
624 ! output of dispersion estimates
625       double precision sig2,sig2d,nu2,nu2d,theta
626 ! d1,d2: length of data matrix
627 ! m, Ji, maxj: # of clusters, individuals, max individuals
628       integer d1,d2,m,Ji(m),maxj,ncovs1,ncovs2,nr,ns,nsd
629 ! Sensitivity matrix (for output)
630       double precision Smat(ncovs1+ncovs2,ncovs1+ncovs2)
631 ! matrix of alpharh values
632       double precision as(nr,ns), asd(nr,nsd)
633       double precision alphars(nr,ns),alpharsd(nr,nsd)
634 ! matrix of coefficient estimates
635       double precision betahat(ncovs1),betadhat(ncovs2)
636 ! matrix of indices
637       integer index1(d1,6),index2(d2,6)
638 ! covariates: matrix of covariates for each index entry
639       double precision Z(d1,ncovs1)
640       double precision Zd(d2,ncovs2)
641       double precision times1(d1),times2(d2)
642 ! variables to hold the current values
643       double precision bZd,bZ,mu,eta
644       double precision w,wj,time,thisA
645       integer i,ind,r,iind0,iind1,iind0d,iind1d,s
646 ! Constants for column indices
647       parameter(icol=1, jcol=2, kcol=3, ircol=4, ismincol=5, ismaxcol=6)
648       
649       double precision S11(ncovs1,ncovs1)
650       double precision S12(ncovs1,ncovs2)
651       double precision S13(ncovs2,ncovs2)
652       double precision S21(ncovs1),S22(ncovs2)
653       double precision S31(ncovs1),S32(ncovs2)      
654       double precision S2(ncovs1+ncovs2)
655       double precision S3(ncovs1+ncovs2)
656       double precision Sii(ncovs1+ncovs2,ncovs1+ncovs2)
657 !      double precision Sinv(g+gd,g+gd)
658       
659       double precision Sm(maxj)
660       double precision Se(maxj)
661       double precision Smx(maxj,ncovs1)
662       double precision Sex(maxj,ncovs2)
663       
664       double precision phi(maxj),Z2(ncovs1),Z2d(ncovs2)
665       
666 ! indices for the start and end of the i portion of the data matrix      
667       iind0=0
668       iind1=0
669       iind0d=0
670       iind1d=0
671       
672       call dzero(Smat,(ncovs1+ncovs2)*(ncovs1+ncovs2))
673       
674 ! Loop over the values of i
675       do 200 i=1,m
676       
677         call dzero(S11,ncovs1*ncovs1)
678         call dzero(S12,ncovs1*ncovs2)
679         call dzero(S13,ncovs2*ncovs2)
680         call dzero(S21,ncovs1)
681         call dzero(S31,ncovs1)
682         call dzero(S22,ncovs2)
683         call dzero(S32,ncovs2)       
684         call dzero(S2,ncovs1+ncovs2)
685         call dzero(S3,ncovs1+ncovs2)
686         call dzero(Sm,maxj)
687         call dzero(Se,maxj)
688         call dzero(Smx,maxj*ncovs1)
689         call dzero(Sex,maxj*ncovs2)
690         call dzero(Z2,ncovs1)
691         call dzero(Z2d,ncovs2)
692 
693         iind0=iind1+1
694         iind1=iind1+1
695         iind0d=iind1d+1
696         iind1d=iind1d+1
697 
698       
699 !     while loop to find start and end indices in ij for cluster i
700  210    if((index1(iind1,icol).eq.i).and.(iind1.le.d1)) then
701             iind1=iind1+1
702             goto 210
703         endif      
704  211    if((index2(iind1d,icol).eq.i).and.(iind1d.le.d2)) then
705             iind1d=iind1d+1
706             goto 211
707         endif     
708 
709 !       Compute the Sm, Se vectors and Smx,Sex matrix
710 !       Sm=Sum(mu_ij*) for j=1..Ji
711 !       Smx=Sum(mu_ijkh*x_ijkh) for j=1..Ji
712         iind1=iind1-1                
713         do 220 ind=iind0,iind1
714             j=index1(ind,jcol)
715             smin=index1(ind,ismincol)
716             smax=index1(ind,ismaxcol)
717             r=index1(ind,ircol)
718             time=times1(ind)
719             bZ=ddot(ncovs1,Z(ind,1),d1,betahat,1)
720             if(smax .GT. 0) then
721                 do 230 s=smin,smax
722                     thisA=0.d0
723                     call A(thisA,time,as,r,s,nr,ns)
724                     mu=exp(bZ)*alphars(r,s)*thisA
725                     Sm(j)=Sm(j)+mu                        
726                     call daxpy(ncovs1,mu,Z(ind,1),d1,Smx(j,1),maxj)
727  230            continue
728             endif  
729  220    continue 
730         iind1d=iind1d-1                
731         do 221 ind=iind0d,iind1d              
732             j=index2(ind,jcol)
733             smin=index2(ind,ismincol)
734             smax=index2(ind,ismaxcol)
735             r=index2(ind,ircol)
736             time=times2(ind)
737             bZ=ddot(ncovs2,Zd(ind,1),d2,betadhat,1)
738             if((smax .GT. 0)) then
739                 do 240 s=smin,smax
740                     thisA=0.d0
741                     call A(thisA,time,asd,r,s,nr,nsd)
742                     eta=exp(bZ)*alpharsd(r,s)*thisA
743                     Se(j)=Se(j)+eta                        
744                     call daxpy(ncovs2,eta,Zd(ind,1),d2,Sex(j,1),maxj)
745  240            continue
746             endif   
747  221    continue 
748         
749  
750 ! Compute phi
751         do 300 j=1,Ji(i)
752             phi(j)=(1.0d0+nu2*Sm(j))*(1.0d0+nu2d*Se(j))
753      $      -theta*theta*Sm(j)*Se(j)   
754  300    continue    
755             
756       
757 ! Compute the various components of the matrix (pass 2)
758         do 400 ind=iind0,iind1
759             j=index1(ind,jcol)
760             smin=index1(ind,ismincol)
761             smax=index1(ind,ismaxcol)
762             r=index1(ind,ircol)
763             time=times1(ind)
764             bZ=ddot(ncovs1,Z(ind,1),d1,betahat,1)
765             call dcopy(ncovs1,Z(ind,1),d1,Z2(1),1)
766             
767             if(smax .GT. 0) then
768                 do 410 s=smin,smax
769                     thisA=0.d0
770                     call A(thisA,time,as,r,s,nr,ns)
771                     mu=exp(bZ)*alphars(r,s)*thisA
772       
773                     call dger(ncovs1,ncovs1,mu,Z2,1,Z2,1,S11,ncovs1)
774                     
775                     w=1.0d0/(1.0d0+wij(i,j)*Sm(j))
776                     call daxpy(ncovs1,w*mu,Z(ind,1),d1,S21(1),1)
777                     
778                     w= -theta*Se(j)/phi(j)
779                     call daxpy(ncovs1,w*mu,Z(ind,1),d1,S31(1),1)
780  410            continue
781             endif
782  400    continue
783         do 401 ind=iind0d,iind1d
784             j=index2(ind,jcol)
785             smin=index2(ind,ismincol)
786             smax=index2(ind,ismaxcol)
787             r=index2(ind,ircol)
788             time=times2(ind)
789             bZ=ddot(ncovs2,Zd(ind,1),d2,betadhat,1)
790             call dcopy(ncovs,Zd(ind,1),d2,Z2(1),1)
791             if((smax .GT. 0)) then
792                 do 420 s=smin,smax
793                     thisA=0.d0
794                     call A(thisA,time,asd,r,s,nr,nsd)
795                     eta=exp(bZ)*alpharsd(r,s)*thisA
796 
797                     call dger(ncovs2,ncovs2,eta,Z2,1,Z2,1,S13,ncovs2)
798                     
799                     w= -theta*Sm(j)/phi(j)
800                     call daxpy(ncovs2,w*eta,Zd(ind,1),d2,S22(1),1)
801                     
802                     w=1.0d0/(1.0d0+zij(i,j)*Se(j))
803                     call daxpy(ncovs2,w*eta,Zd(ind,1),d2,S32(1),1)
804  420            continue
805             endif
806  401    continue
807   
808           
809 
810 ! Finish computation of the components via appropriate rank 1 transforms
811         do 500 j=1,Ji(i)
812          w=-wij(i,j)/(1.0d0+wij(i,j)*Sm(j))
813        call dger(ncovs1,ncovs1,w,Smx(j,1),maxj,Smx(j,1),maxj,S11,ncovs1)
814          w=-theta/phi(j)
815        call dger(ncovs1,ncovs2,w,Smx(j,1),maxj,Sex(j,1),maxj,S12,ncovs1)
816          w=-zij(i,j)/(1.0d0+zij(i,j)*Se(j))
817        call dger(ncovs2,ncovs2,w,Sex(j,1),maxj,Sex(j,1),maxj,S13,ncovs2)
818  500    continue 
819         call dcopy(ncovs1,S21,1,S2,1)
820         call dcopy(ncovs2,S22,1,S2(ncovs1+1),1)
821         call dcopy(ncovs1,S31,1,S3,1)
822         call dcopy(ncovs2,S32,1,S3(ncovs1+1),1)
823  
824  
825       
826 ! Concatenate the matrices S11,S12,S13 correctly into Si
827 ! At the same time, apply the minus sign
828         do 610 jj=1,ncovs1
829             do 620 ii=1,ncovs1
830                 Sii(ii,jj)=-S11(ii,jj)
831  620        continue
832  610    continue        
833         do 630 jj=1,ncovs2
834             do 640 ii=1,ncovs1
835                 Sii(ii,jj+ncovs1)=-S12(ii,jj)
836                 Sii(jj+ncovs1,ii)=-S12(ii,jj)
837  640        continue
838  630    continue
839         do 650 jj=1,ncovs2
840             do 660 ii=1,ncovs2
841                 Sii(ii+ncovs1,jj+ncovs1)=-S13(ii,jj)
842  660        continue
843  650    continue
844  
845 
846 ! Compute the ith sensitivity matrix with several rank 1 operations
847         wj=wi(i)*sig2*(1.0d0+sig2d*si(i))
848         call dsyr('u',ncovs1+ncovs2,wj,S2,1,Sii,ncovs1+ncovs2)
849         wj=-wi(i)*sig2*sig2d*qi(i)
850         call dsyr2('u',ncovs1+ncovs2,wj,S2,1,S3,1,Sii,ncovs1+ncovs2)
851         wj=wi(i)*sig2d*(1.0d0+sig2*pi(i))
852         call dsyr('u',ncovs1+ncovs2,wj,S3,1,Sii,ncovs1+ncovs2)
853         
854 
855 ! Add this matrix to the total
856         call daxpy((ncovs1+ncovs2)*(ncovs1+ncovs2),1.d0,Sii,1,Smat,1)
857  
858  200  continue    
859 
860 ! Make the matrix symmetrical again (copy U to L)
861       do 700 j=1,ncovs1+ncovs2-1
862          do 710 i=j+1,ncovs1+ncovs2
863             Smat(i,j)=Smat(j,i)
864  710     continue
865  700  continue      
866       
867       end