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