./blupsurv/src/bivrec.f

1. fortran/fmkmrs

fortran/fmkmrs [ Methods ]

NAME

   fmkmrs --- terms for profile likelihood

FUNCTION

Compute terms needed in the profile likelihood and gradient. See also the R function makemrs.

SYNOPSIS

51       subroutine fmkmrs(betahat,index,times,Z,as,Uijmat,d,ncovs,nr,ns,
52      $  m,maxj,mrs,mrsgr,mkgr)

INPUTS

    betahat    regression parameters (length p)
    index      d x 6 matrix of indices i,j,k,r,smin,smax 
    delta      length d vector of event indicators (double)
    times      length d vector of actual time span for each row (double)
    Z          d x ncovs matrix of covariates (as a double row vector)
    as         nr x ns matrix of discretization breakpoints (double) 
    Uijmat     m x maxj matrix of frailty estimates (double row vector)
    d          number of rows of Z
    ncovs      number of covariates
    nr         number of strata
    ns         number of breakpoints
    m          number of clusters
    maxj       largest cluster size
    mrs        nr x ns matrix for output
    mkgr       integer indicator for whether to make gradient

OUTPUTS

    mrsgr      nr x ns x ncovs gradient for output

SOURCE

54       integer d,ncovs,nr,ns,nk,m,maxj,mkgr
55       integer index(d,6)
56       double precision betahat(ncovs)
57       double precision times(d)
58       double precision Z(d,ncovs)
59       double precision Uijmat(m,maxj)
60       double precision lp,time,pred1,pred2,thisA
61       double precision mrs(nr,ns)
62       double precision mrsgr(nr,ns,ncovs)
63       double precision as(nr,ns)
64       
65       integer ind,i,j,k,r,s,smax,smin,maxs    
66       parameter(icol=1, jcol=2, kcol=3, ircol=4, ismincol=5, ismaxcol=6)
67       
68 !     clear output arrays
69       call dzero(mrs,nr*ns)
70       call dzero(mrsgr,nr*ns*ncovs)
71 
72       maxs=(ns-1)*95/100      
73       do 100 ind=1,d
74         i=index(ind,icol)
75         j=index(ind,jcol)
76         k=index(ind,kcol)
77         r=index(ind,ircol)
78         smin=index(ind,ismincol)
79         smax=index(ind,ismaxcol)
80         time=times(ind)
81         lp=ddot(ncovs,Z(ind,1),d,betahat,1)
82 !       ignore the last 5% of intervals
83 !        if(smax.GT.maxs) then
84 !            smax=maxs
85 !        endif
86         pred1=Uijmat(i,j)*exp(lp)  
87         do 110 s=smin,smax
88             thisA=0.d0
89             call A(thisA,time,as,r,s,nr,ns)
90             pred2=pred1*thisA
91             mrs(r,s)=mrs(r,s)+pred2
92             if(mkgr.EQ.1) then
93                 call daxpy(ncovs,pred2,Z(ind,1),d,mrsgr(r,s,1),nr*ns)
94             endif
95  110    continue
96                 
97  100  continue
98       end

TABLE OF CONTENTS

fortran/fmkmrs [ Methods ]