1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
/* $Id: coxscore.c 11357 2009-09-04 15:22:46Z therneau $
**
** Compute the score residuals for a Cox model
**
** Input
** nx number of subjects
** nvarx number of variables in the covariance matrix
** y matrix of time and status values
** strata =1 for the last obs of each strata
** covar2 the matrix of covariates, rows=variables, columns=subjects
** (the S executive stores matrices in the Fortran ordering)
** score the vector of subject scores, i.e., exp(beta*z)
** weights case weight
** method ==1 for efron method
**
** Output
** resid a matrix of the same shape as x
**
** Scratch
** scratch, from which a and a2 are carved
**
** Data must be sorted by strata, ascending time within strata, death before
** censor within time.
*/
#include <stdio.h>
#include "survS.h"
#include "survproto.h"
void coxscore(Sint *nx, Sint *nvarx, double *y,
double *covar2, Sint *strata, double *score,
double *weights, Sint *method, double *resid2,
double *scratch)
{
int i,j, k;
double temp;
int n, nvar;
double deaths;
int dd;
double *time, *status;
double *a, *a2;
double denom=0, e_denom;
double risk;
double **covar;
double **resid;
double hazard, meanwt;
double downwt, temp2;
double mean;
n = *nx;
nvar = *nvarx;
time = y;
status = y+n;
a = scratch;
a2 = a+nvar;
/*
** Set up the ragged array
*/
covar= dmatrix(covar2, n, nvar);
resid= dmatrix(resid2, n, nvar);
e_denom=0;
deaths=0;
meanwt=0;
for (i=0; i<nvar; i++) a2[i] =0;
strata[n-1] =1; /*failsafe */
for (i=n-1; i >=0; i--) {
if (strata[i]==1) {
denom =0;
for (j=0; j<nvar; j++) a[j] =0;
}
risk = score[i] * weights[i];
denom += risk;
if (status[i]==1) {
deaths++;
e_denom += risk;
meanwt += weights[i];
for (j=0; j<nvar; j++) a2[j] += risk*covar[j][i];
}
for (j=0; j<nvar; j++) {
a[j] += risk * covar[j][i];
resid[j][i] =0;
}
if (deaths>0 && (i==0 || strata[i-1]==1 || time[i]!=time[i-1])){
/* last obs of a set of tied death times */
if (deaths <2 || *method==0) {
hazard = meanwt/denom;
for (j=0; j<nvar; j++) {
temp = (a[j]/denom); /* xbar */
for (k=i; k<n; k++) {
temp2 = covar[j][k] - temp;
if (time[k]==time[i] && status[k]==1)
resid[j][k] += temp2;
resid[j][k] -= temp2* score[k] * hazard;
if (strata[k]==1) break;
}
}
}
else { /* the harder case */
meanwt /= deaths;
for (dd=0; dd<deaths; dd++) {
downwt = dd/deaths;
temp = denom - downwt* e_denom;
hazard = meanwt/temp;
for (j=0; j<nvar; j++) {
mean = (a[j] - downwt*a2[j])/ temp;
for (k=i; k<n; k++) {
temp2 = covar[j][k] - mean;
if (time[k]==time[i] && status[k]==1) {
resid[j][k] += temp2/deaths;
resid[j][k] -= temp2 * score[k] * hazard *
(1 - downwt);
}
else resid[j][k]-= temp2*score[k] * hazard;
if (strata[k]==1) break;
}
}
}
}
e_denom =0;
deaths =0;
meanwt =0;
for (j=0; j<nvar; j++) a2[j] =0;
}
}
}
|
