Revision 858421cbbaf89fbe446d88489ffe8534a90cb0f4 authored by Venkatraman E. Seshan on 20 November 2009, 00:00:00 UTC, committed by Gabor Csardi on 20 November 2009, 00:00:00 UTC
1 parent af5dcd9
calogrank.f
c This subroutine is adapted from the code for the JRSS-B paper
c it does the conditional logrank in a 3 group setting
c need to change it to a general k-group setting
c create k indicator variables and the k vector of U-statistics
c
c n = number of observations, ng = number of groups
c p = number of covariates (maximum of 3) -- moved to R wrapper?
c osts, ogrp, ocov are the status, group and covariate variables
c a0 and a1 (per grp) vectors correspond to the same in the paper
c vectors Vii, Vij, Vji, Vjj are the group specific U-statistics
c Vidot, Vdotj, Vijij are the row, column and product sums
c igrp, jgrp are the group indicator for observation i & j
c lrmn, lrvar are the test statistic vector and variance
c bb is the bandwidth for the kernel
subroutine uclrst(n, ng, p, osts, ogrp, ocov, a0, a1, xi, xj,
1 Vii, Vij, Vji, Vjj, Vidot, Vdotj, Vijij, igrp, jgrp, lrmn,
2 lrvar, bb)
integer n, ng, p, ogrp(n)
double precision osts(n), ocov(n,p), a0(n), a1(ng,n), xi(p),
1 xj(p), Vii(ng), Vij(ng), Vji(ng), Vjj(ng), Vidot(ng,n),
2 Vdotj(ng,n), Vijij(ng,ng), igrp(ng), jgrp(ng), lrmn(ng),
3 lrvar(ng,ng), bb(p)
integer i,j,ic,g,gr,gc
double precision rn,kernel,kwtij
external kernel
rn = dfloat(n)
c calculate the a0 and a1 vector for all time points
c the vij calculations depend only on these
do 25 i = 1,n
a0(i) = 0.0
do 21 g = 1,ng
a1(g,i) = 0.0
21 continue
do 23 j = i,n
do 22 ic = 1,p
xi(ic) = ocov(i,ic)
xj(ic) = ocov(j,ic)
22 continue
kwtij = kernel(p,xj,xi,bb)
c a1 changes only for the group for jth observation
a0(i) = a0(i) + kwtij
a1(ogrp(j), i) = a1(ogrp(j), i) + kwtij
23 continue
a0(i) = a0(i)/rn
do 24 g = 1,ng
a1(g,i) = a1(g,i)/rn
24 continue
25 continue
c now calculate Vij to compute the statistic and its variance
do 40 i = 1,n
do 31 ic = 1,p
xi(ic) = ocov(i,ic)
31 continue
c set the indicator of the ith observation
igrp(ogrp(i)) = 1.0
do 39 j = i,n
do 32 ic = 1,p
xj(ic) = ocov(j,ic)
32 continue
kwtii = kernel(p,xi,xi,bb)
kwtij = kernel(p,xj,xi,bb)
c set the indicator of the jth observation
jgrp(ogrp(j)) = 1.0
do 33 g = 1,ng
c compute the Vij's
Vii(g) = osts(i)*(igrp(g) - a1(g,i)/a0(i))*
1 (1.0-kwtii/a0(i))/rn
if (j .gt. i) then
Vij(g) = osts(i)*(igrp(g) - a1(g,i)/a0(i) -
1 kwtij*(jgrp(g) - a1(g,i)/a0(i))/a0(i))/rn
Vji(g) = osts(j)*(jgrp(g) - a1(g,j)/a0(j))/rn
Vjj(g) = osts(j)*(jgrp(g) - a1(g,j)/a0(j))*
1 (1.0-kwtii/a0(j))/rn
else
Vij(g) = Vii(g)
Vji(g) = Vii(g)
Vjj(g) = Vii(g)
endif
33 continue
c unset the indicator of the jth observation
jgrp(ogrp(j)) = 0.0
c compute the row and column sums of Vij's
do 34 g = 1,ng
lrmn(g) = lrmn(g) + Vij(g)
Vidot(g,i) = Vidot(g,i) + Vij(g)
Vdotj(g,j) = Vdotj(g,j) + Vij(g)
if(j .gt. i) then
lrmn(g) = lrmn(g) + Vji(g)
Vidot(g,j) = Vidot(g,j) + Vji(g)
Vdotj(g,i) = Vdotj(g,i) + Vji(g)
endif
34 continue
c compute the Vij product contributions to the variance
if (i.lt.j) then
c these are the i neq j terms
do 36 gr = 1,ng
do 35 gc = gr,ng
if (gr .eq. gc) then
Vijij(gr,gr) = Vijij(gr,gr) +
1 (Vij(gr) + Vji(gr))*(Vij(gr) +
2 Vji(gr) + 2*Vii(gr) + 2*Vjj(gr))
else
Vijij(gr,gc) = Vijij(gr,gc) +
1 (Vij(gr) + Vji(gr))*(Vij(gc) + Vji(gc)) +
2 (Vij(gr) + Vji(gr))*(Vii(gc) + Vjj(gc)) +
3 (Vij(gc) + Vji(gc))*(Vii(gr) + Vjj(gr))
4
endif
35 continue
36 continue
else
c these are the ii terms
do 38 gr = 1,ng
do 37 gc = gr,ng
if (gr .eq. gc) then
Vijij(gr,gr) = Vijij(gr,gr) + 3*Vii(gr)**2
else
Vijij(gr,gc) = Vijij(gr,gc) +
1 3*Vii(gr)*Vii(gc)
endif
37 continue
38 continue
endif
39 continue
c unset the indicator of the ith observation
igrp(ogrp(i)) = 0.0
40 continue
c finish the variance computation
do 50 gr = 1,ng
do 49 gc = gr,ng
do 48 i = 1,n
lrvar(gr,gc) = lrvar(gr,gc) + (Vidot(gr,i)+Vdotj(gr,i))*
1 (Vidot(gc,i)+Vdotj(gc,i))
48 continue
lrvar(gr,gc) = lrvar(gr,gc) - Vijij(gr,gc)
if (gc .gt. gr) lrvar(gc,gr) = lrvar(gr,gc)
49 continue
50 continue
do 60 i = 1,n
if (osts(i) .eq. 1) then
igrp(ogrp(i)) = igrp(ogrp(i)) + 1.0
do 55 g = 1,ng
jgrp(g) = jgrp(g) + a1(g,i)/a0(i)
55 continue
endif
60 continue
return
end
double precision function kernel(p,xx,xi,bb)
integer i,p
double precision xx(p),xi(p),bb(p),uu
kernel = 1.0
do 10 i = 1,p
uu = ((xx(i)-xi(i))/bb(i))**2
kernel=kernel*(1.0/bb(i))*exp(-uu/2)
10 continue
return
end
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...