TABLE OF CONTENTS
fortran/fmksens2full [ Methods ]
NAME
fmksens2full --- compute sensitivity matrix
FUNCTION
Construct the full sensitivity matrix at the end of the procedure. This computes sensitivity to both regression coefficient and baseline hazard parameters See also the R function makeSens2.
SYNOPSIS
888 subroutine fmksens2full(Smat,index1,index2,Z,Zd,alphars,alpharsd, 889 $ as,asd,betahat,betadhat,times1,times2,pi,qi,ri,si,wi,wij,zij, 890 $ sig2,sig2d,nu2,nu2d,theta,ncovs1,ncovs2,nr,ns,nsd,np,npd, 891 $ d1,d2,m,Ji,maxj,Kcum,Kdcum)
INPUTS
Too many to list, but analogous to those in other routines in this file. Inputs appended with d refer to process 2. np, npd number of parrameters for each of the processes
OUTPUTS
Smat matrix of size ncovs1+ncovs2 x ncovs2+ncovs2
sensitivity matrix
SOURCE
894 ! output of frailty estimation 895 double precision pi(m),qi(m),ri(m),si(m),wi(m), 896 $ wij(m,maxj),zij(m,maxj) 897 ! output of dispersion estimates 898 double precision sig2,sig2d,nu2,nu2d,theta 899 ! d1,d2: length of data matrix 900 ! m, Ji, maxj: # of clusters, individuals, max individuals 901 integer d1,d2,m,Ji(m),maxj,ncovs1,ncovs2,nr,ns,nsd,np,npd 902 ! Sensitivity matrix (for output) 903 double precision as(nr,ns), asd(nr,nsd) 904 double precision Smat(np+npd,np+npd) 905 ! matrix of alpharh values 906 double precision alphars(nr,ns),alpharsd(nr,nsd) 907 ! matrix of coefficient estimates 908 double precision betahat(ncovs1),betadhat(ncovs2) 909 ! matrix of indices 910 integer index1(d1,6),index2(d2,6) 911 ! covariates: matrix of covariates for each index entry 912 double precision Z(d1,ncovs1) 913 double precision Zd(d2,ncovs2) 914 double precision times1(d1),times2(d2) 915 integer Kcum(nr+1), Kdcum(nr+1) 916 ! variables to hold the current values 917 double precision bZd,bZ,mu,eta,w,wj,time,thisA 918 integer i,ind,r,iind0,iind1,iind0d,iind1d,s,bst,bdst,sind 919 ! real*4 timex(2),ExecTime(3),LastTime(3),tottime,tottimex 920 921 ! Constants for column indices 922 parameter(icol=1, jcol=2, kcol=3, ircol=4, ismincol=5, ismaxcol=6) 923 924 double precision S11(np,np) 925 double precision S12(np,npd) 926 double precision S13(npd,npd) 927 double precision S21(np),S22(npd) 928 double precision S31(np),S32(npd) 929 double precision S2(np+npd) 930 double precision S3(np+npd) 931 double precision Sii(np+npd,np+npd) 932 ! double precision Sinv(g+gd,g+gd) 933 934 double precision Sm(maxj) 935 double precision Se(maxj) 936 double precision Smx(maxj,np) 937 double precision Sex(maxj,npd) 938 939 double precision phi(maxj),Z2(np),Z2d(npd) 940 941 ! indices for the start and end of the i portion of the data matrix 942 iind0=0 943 iind1=0 944 iind0d=0 945 iind1d=0 946 947 bst=Kcum(nr+1) 948 bdst=Kdcum(nr+1) 949 950 ! call etime(timex,tottimex) 951 ! call etime(timex,tottime) 952 ! LastTime(1) = tottimex 953 ! LastTime(2) = timex(1) 954 ! LastTime(3) = timex(2) 955 956 call dzero(Smat,(np+npd)*(np+npd)) 957 958 ! Loop over the values of i 959 do 200 i=1,m 960 961 ! if(i.eq.1) then 962 ! call etime(timex,tottimex) 963 ! ExecTime(1) = tottimex-LastTime(1) 964 ! ExecTime(2) = timex(1)-LastTime(2) 965 ! ExecTime(3) = timex(2)-LastTime(3) 966 ! LastTime(1)=tottimex 967 ! LastTime(2)=timex(1) 968 ! LastTime(3)=timex(2) 969 ! call realpr('Startup',-1,ExecTime,3) 970 ! endif 971 972 call dzero(S11,np*np) 973 call dzero(S12,np*npd) 974 call dzero(S13,npd*npd) 975 call dzero(S21,np) 976 call dzero(S31,np) 977 call dzero(S22,npd) 978 call dzero(S32,npd) 979 call dzero(S2,np+npd) 980 call dzero(S3,np+npd) 981 call dzero(Sm,maxj) 982 call dzero(Se,maxj) 983 call dzero(Smx,maxj*np) 984 call dzero(Sex,maxj*npd) 985 call dzero(Z2,np) 986 call dzero(Z2d,npd) 987 988 iind0=iind1+1 989 iind1=iind1+1 990 iind0d=iind1d+1 991 iind1d=iind1d+1 992 993 ! while loop to find start and end indices in ij for cluster i 994 210 if((index1(iind1,icol).eq.i).and.(iind1.le.d1)) then 995 iind1=iind1+1 996 goto 210 997 endif 998 211 if((index2(iind1d,icol).eq.i).and.(iind1d.le.d2)) then 999 iind1d=iind1d+1 1000 goto 211 1001 endif 1002 1003 ! Compute the Sm, Se vectors and Smx,Sex matrix 1004 ! Sm=Sum(mu_ij*) for j=1..Ji 1005 ! Smx=Sum(mu_ijkh*x_ijkh) for j=1..Ji 1006 iind1=iind1-1 1007 do 220 ind=iind0,iind1 1008 j=index1(ind,jcol) 1009 smin=index1(ind,ismincol) 1010 smax=index1(ind,ismaxcol) 1011 r=index1(ind,ircol) 1012 time=times1(ind) 1013 bZ=ddot(ncovs1,Z(ind,1),d1,betahat,1) 1014 if(smax .GT. 0) then 1015 do 230 s=smin,smax 1016 thisA=0.d0 1017 call A(thisA,time,as,r,s,nr,ns) 1018 mu=exp(bZ)*alphars(r,s)*thisA 1019 Sm(j)=Sm(j)+mu 1020 Smx(j,Kcum(r)+s-1)=Smx(j,Kcum(r)+s-1)+mu 1021 call daxpy(ncovs1,mu,Z(ind,1),d1,Smx(j,bst),maxj) 1022 230 continue 1023 endif 1024 220 continue 1025 iind1d=iind1d-1 1026 do 221 ind=iind0d,iind1d 1027 j=index2(ind,jcol) 1028 smin=index2(ind,ismincol) 1029 smax=index2(ind,ismaxcol) 1030 r=index2(ind,ircol) 1031 time=times2(ind) 1032 bZ=ddot(ncovs2,Zd(ind,1),d2,betadhat,1) 1033 if((smax .GT. 0)) then 1034 do 240 s=smin,smax 1035 thisA=0.d0 1036 call A(thisA,time,asd,r,s,nr,nsd) 1037 eta=exp(bZ)*alpharsd(r,s)*thisA 1038 Se(j)=Se(j)+eta 1039 Sex(j,Kdcum(r)+s-1)=Sex(j,Kdcum(r)+s-1)+eta 1040 call daxpy(ncovs2,eta,Zd(ind,1),d2,Sex(j,bdst),maxj) 1041 240 continue 1042 endif 1043 221 continue 1044 1045 1046 ! Compute phi 1047 do 300 j=1,Ji(i) 1048 phi(j)=(1.0d0+nu2*Sm(j))*(1.0d0+nu2d*Se(j)) 1049 $ -theta*theta*Sm(j)*Se(j) 1050 300 continue 1051 1052 ! Compute the various components of the matrix (pass 2) 1053 do 400 ind=iind0,iind1 1054 j=index1(ind,jcol) 1055 smin=index1(ind,ismincol) 1056 smax=index1(ind,ismaxcol) 1057 r=index1(ind,ircol) 1058 time=times1(ind) 1059 bZ=ddot(ncovs1,Z(ind,1),d1,betahat,1) 1060 ! call dcopy(ncovs,Z(ind,1),d1,Z2(bst),1) 1061 1062 if(smax .GT. 0) then 1063 do 410 s=smin,smax 1064 thisA=0.d0 1065 call A(thisA,time,as,r,s,nr,ns) 1066 mu=exp(bZ)*alphars(r,s)*thisA 1067 sind=Kcum(r)+s-1 1068 1069 S11(sind,sind)=S11(sind,sind)+mu 1070 call daxpy(ncovs1,mu,Z(ind,1),d1,S11(sind,bst),np) 1071 ! call daxpy(ncovs,mu,Z(ind,1),d1,S11(bst,sind),1) 1072 ! call dger(np,np,mu,Z2,1,Z2,1,S11,np) 1073 call dsyr('U',ncovs1,mu,Z(ind,1),d1,S11(bst,bst),np) 1074 1075 w=1.0d0/(1.0d0+wij(i,j)*Sm(j)) 1076 S21(sind)=S21(sind)+w*mu 1077 call daxpy(ncovs1,w*mu,Z(ind,1),d1,S21(bst),1) 1078 1079 w= -theta*Se(j)/phi(j) 1080 S31(sind)=S31(sind)+w*mu 1081 call daxpy(ncovs1,w*mu,Z(ind,1),d1,S31(bst),1) 1082 410 continue 1083 endif 1084 400 continue 1085 do 401 ind=iind0d,iind1d 1086 j=index2(ind,jcol) 1087 smin=index2(ind,ismincol) 1088 smax=index2(ind,ismaxcol) 1089 r=index2(ind,ircol) 1090 time=times2(ind) 1091 bZ=ddot(ncovs2,Zd(ind,1),d2,betadhat,1) 1092 ! call dcopy(ncovs,Zd(ind,1),d2,Z2d(bdst),1) 1093 if((smax .GT. 0)) then 1094 do 420 s=smin,smax 1095 thisA=0.d0 1096 call A(thisA,time,asd,r,s,nr,nsd) 1097 eta=exp(bZ)*alpharsd(r,s)*thisA 1098 sind=Kdcum(r)+s-1 1099 1100 S13(sind,sind)=S13(sind,sind)+eta 1101 call daxpy(ncovs2,eta,Zd(ind,1),d2,S13(sind,bdst),npd) 1102 ! call daxpy(ncovs,eta,Zd(ind,1),d2,S13(bdst,sind),1) 1103 call dger(npd,npd,eta,Z2d,1,Z2d,1,S13,npd) 1104 call dsyr('U',ncovs2,eta,Zd(ind,1),d2,S13(bdst,bdst),npd) 1105 1106 1107 w= -theta*Sm(j)/phi(j) 1108 S22(sind)=S22(sind)+w*eta 1109 call daxpy(ncovs2,w*eta,Zd(ind,1),d2,S22(bdst),1) 1110 1111 w=1.0d0/(1.0d0+zij(i,j)*Se(j)) 1112 S32(sind)=S32(sind)+w*eta 1113 call daxpy(ncovs2,w*eta,Zd(ind,1),d2,S32(bdst),1) 1114 420 continue 1115 endif 1116 401 continue 1117 1118 1119 ! Finish computation of the components via appropriate rank 1 transforms 1120 do 500 j=1,Ji(i) 1121 w=-wij(i,j)/(1.0d0+wij(i,j)*Sm(j)) 1122 ! call dger(np,np,w,Smx(j,1),maxj,Smx(j,1),maxj,S11,np) 1123 call dsyr('U',np,w,Smx(j,1),maxj,S11,np) 1124 w=-theta/phi(j) 1125 call dger(np,npd,w,Smx(j,1),maxj,Sex(j,1),maxj,S12,np) 1126 w=-zij(i,j)/(1.0d0+zij(i,j)*Se(j)) 1127 ! call dger(npd,npd,w,Sex(j,1),maxj,Sex(j,1),maxj,S13,npd) 1128 call dsyr('U',npd,w,Sex(j,1),maxj,S13,npd) 1129 1130 500 continue 1131 call dcopy(np,S21,1,S2,1) 1132 call dcopy(npd,S22,1,S2(np+1),1) 1133 call dcopy(np,S31,1,S3,1) 1134 call dcopy(npd,S32,1,S3(np+1),1) 1135 1136 ! Concatenate the matrices S11,S12,S13 correctly into Si 1137 ! At the same time, apply the minus sign 1138 do 610 jj=1,np 1139 do 620 ii=1,jj 1140 Sii(ii,jj)=-S11(ii,jj) 1141 620 continue 1142 610 continue 1143 do 630 jj=1,npd 1144 do 640 ii=1,np 1145 Sii(ii,jj+np)=-S12(ii,jj) 1146 ! Sii(jj+np,ii)=-S12(ii,jj) 1147 640 continue 1148 630 continue 1149 do 650 jj=1,npd 1150 do 660 ii=1,jj 1151 Sii(ii+np,jj+np)=-S13(ii,jj) 1152 660 continue 1153 650 continue 1154 1155 1156 ! Compute the ith sensitivity matrix with several rank 1 operations 1157 wj=wi(i)*sig2*(1.0d0+sig2d*si(i)) 1158 call dsyr('u',np+npd,wj,S2,1,Sii,np+npd) 1159 wj=-wi(i)*sig2*sig2d*qi(i) 1160 call dsyr2('u',np+npd,wj,S2,1,S3,1,Sii,np+npd) 1161 wj=wi(i)*sig2d*(1.0d0+sig2*pi(i)) 1162 call dsyr('u',np+npd,wj,S3,1,Sii,np+npd) 1163 1164 ! Add this matrix to the total 1165 call daxpy((np+npd)*(np+npd),1.d0,Sii,1,Smat,1) 1166 1167 200 continue 1168 1169 ! Make the matrix symmetrical again (copy U to L) 1170 do 700 j=1,np+npd-1 1171 do 710 i=j+1,np+npd 1172 Smat(i,j)=Smat(j,i) 1173 710 continue 1174 700 continue 1175 1176 end