← Back to branch summary

~ubuntu-branches/ubuntu/lucid/python-scipy/lucid

« back to all changes in this revision

Viewing changes to Lib/sandbox/models/cox.py

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/models/cox.py
 
1
import shutil
 
2
import tempfile
 
3
import numpy as N
 
4
from scipy.sandbox.models import survival, model
 
5
 
 
6
class discrete:
 
7
 
 
8
    """
 
9
    A simple little class for working with discrete random vectors.
 
10
    """
 
11
 
 
12
    def __init__(self, x, w=None):
 
13
        self.x = N.squeeze(x)
 
14
        if self.x.shape == ():
 
15
            self.x = N.array([self.x])
 
16
        self.n = self.x.shape[0]
 
17
        if w is None:
 
18
            w = N.ones(self.n, N.float64) 
 
19
        else:
 
20
            if w.shape[0] != self.n:
 
21
                raise ValueError, 'incompatible shape for weights w'
 
22
            if N.any(N.less(w, 0)):
 
23
                raise ValueError, 'weights should be non-negative'
 
24
        self.w = w / w.sum()
 
25
 
 
26
    def mean(self, f=None):
 
27
        if f is None:
 
28
            fx = self.x
 
29
        else:
 
30
            fx = f(self.x)
 
31
        return (fx * self.w).sum()
 
32
 
 
33
    def cov(self):
 
34
        mu = self.mean()
 
35
        dx = self.x - N.multiply.outer(mu, self.x.shape[1])
 
36
        return N.dot(dx, N.transpose(dx))
 
37
 
 
38
class observation(survival.right_censored):
 
39
 
 
40
    def __getitem__(self, item):
 
41
        if self.namespace is not None:
 
42
            return self.namespace[item]
 
43
        else:
 
44
            return getattr(self, item)
 
45
 
 
46
    def __init__(self, time, delta, namespace=None):
 
47
        self.namespace = namespace
 
48
        survival.right_censored.__init__(self, time, delta)
 
49
 
 
50
    def __call__(self, formula, time=None, **extra):
 
51
        return formula(namespace=self, time=time, **extra)
 
52
 
 
53
class coxph(model.likelihood_model):
 
54
 
 
55
    def __init__(self, subjects, formula, time_dependent=False):
 
56
        self.subjects, self.formula = subjects, formula
 
57
        self.time_dependent = time_dependent
 
58
        self.initialize(self.subjects)
 
59
 
 
60
    def initialize(self, subjects):
 
61
 
 
62
        self.failures = {}
 
63
        for i in range(len(subjects)):
 
64
            s = subjects[i]
 
65
            if s.delta:
 
66
                if not self.failures.has_key(s.time):
 
67
                    self.failures[s.time] = [i]
 
68
                else:
 
69
                    self.failures[s.time].append(i)
 
70
        
 
71
        self.failure_times = self.failures.keys()
 
72
        self.failure_times.sort()
 
73
 
 
74
    def cache(self):
 
75
        if self.time_dependent:
 
76
            self.cachedir = tempfile.mkdtemp()
 
77
 
 
78
        self.design = {}
 
79
        self.risk = {}
 
80
        first = True
 
81
        
 
82
        for t in self.failures.keys():
 
83
            if self.time_dependent:
 
84
                d = N.array([s(self.formula, time=t)
 
85
                             for s in self.subjects]).astype('<f8')
 
86
                dshape = d.shape
 
87
                dfile = file(tempfile.mkstemp(dir=self.cachedir)[1], 'w')
 
88
                d.tofile(dfile)
 
89
                dfile.close()
 
90
                del(d)
 
91
                self.design[t] = N.memmap(dfile.name,
 
92
                                          dtype=N.dtype('<f8'),
 
93
                                          shape=dshape)
 
94
            elif first:
 
95
                d = N.array([s(self.formula, time=t)
 
96
                             for s in self.subjects]).astype(N.float64)
 
97
                self.design[t] = d
 
98
            else:
 
99
                self.design[t] = d
 
100
            self.risk[t] = N.compress([s.atrisk(t) for s in self.subjects],
 
101
                                      N.arange(self.design[t].shape[0]),axis=-1)
 
102
    def __del__(self):
 
103
        shutil.rmtree(self.cachedir, ignore_errors=True)
 
104
 
 
105
    def logL(self, b, ties='breslow'):
 
106
 
 
107
        logL = 0
 
108
        for t in self.failures.keys():
 
109
            fail = self.failures[t]
 
110
            d = len(fail)
 
111
            risk = self.risk[t]
 
112
            Zb = N.dot(self.design[t], b)
 
113
 
 
114
            logL += Zb[fail].sum()
 
115
 
 
116
            if ties == 'breslow':
 
117
                s = N.exp(Zb[risk]).sum()
 
118
                logL -= N.log(N.exp(Zb[risk]).sum()) * d
 
119
            elif ties == 'efron':
 
120
                s = N.exp(Zb[risk]).sum()
 
121
                r = N.exp(Zb[fail]).sum()
 
122
                for j in range(d):
 
123
                    logL -= N.log(s - j * r / d)
 
124
            elif ties == 'cox':
 
125
                raise NotImplementedError, 'Cox tie breaking method not implemented'
 
126
            else:
 
127
                raise NotImplementedError, 'tie breaking method not recognized'
 
128
        return logL
 
129
 
 
130
    def score(self, b, ties='breslow'):
 
131
 
 
132
        score = 0
 
133
        for t in self.failures.keys():
 
134
            fail = self.failures[t]
 
135
            d = len(fail)
 
136
            risk = self.risk[t]
 
137
            Z = self.design[t]
 
138
 
 
139
            score += Z[fail].sum()
 
140
 
 
141
            if ties == 'breslow':
 
142
                w = N.exp(N.dot(Z, b))
 
143
                rv = discrete(Z[risk], w=w[risk])
 
144
                score -= rv.mean() * d
 
145
            elif ties == 'efron':
 
146
                w = N.exp(N.dot(Z, b))
 
147
                score += Z[fail].sum()
 
148
                for j in range(d):
 
149
                    efron_w = w
 
150
                    efron_w[fail] -= i * w[fail] / d
 
151
                    rv = discrete(Z[risk], w=efron_w[risk])
 
152
                    score -= rv.mean()
 
153
            elif ties == 'cox':
 
154
                raise NotImplementedError, 'Cox tie breaking method not implemented'
 
155
            else:
 
156
                raise NotImplementedError, 'tie breaking method not recognized'
 
157
        # FIXME: score is an int. it has no shape
 
158
        # is it that we shouldn't be using an int above
 
159
        # or that we shouldn't be looking at shape here
 
160
        if score.shape == ():
 
161
            score = N.array([score])
 
162
        return score
 
163
 
 
164
    def information(self, b, ties='breslow'):
 
165
 
 
166
        info = 0
 
167
        score = 0
 
168
        for t in self.failures.keys():
 
169
            fail = self.failures[t]
 
170
            d = len(fail)
 
171
            risk = self.risk[t]
 
172
            Z = self.design[t]
 
173
 
 
174
            if ties == 'breslow':
 
175
                w = N.exp(N.dot(Z, b))
 
176
                rv = discrete(Z[risk], w=w[risk])
 
177
                info += rv.cov()
 
178
            elif ties == 'efron':
 
179
                w = N.exp(N.dot(Z, b))
 
180
                score += Z[fail].sum()
 
181
                for j in range(d):
 
182
                    efron_w = w
 
183
                    efron_w[fail] -= i * w[fail] / d
 
184
                    rv = discrete(Z[risk], w=efron_w[risk])
 
185
                    info += rv.cov()
 
186
            elif ties == 'cox':
 
187
                raise NotImplementedError, 'Cox tie breaking method not implemented'
 
188
            else:
 
189
                raise NotImplementedError, 'tie breaking method not recognized'
 
190
        return score
 
191
 
 
192
if __name__ == '__main__':
 
193
    import numpy.random as R
 
194
    n = 100
 
195
    X = N.array([0]*n + [1]*n)
 
196
    b = 0.4
 
197
    lin = 1 + b*X
 
198
    Y = R.standard_exponential((2*n,)) / lin
 
199
    delta = R.binomial(1, 0.9, size=(2*n,))
 
200
 
 
201
    subjects = [observation(Y[i], delta[i]) for i in range(2*n)]
 
202
    for i in range(2*n):
 
203
        subjects[i].X = X[i]
 
204
 
 
205
    import formula as F
 
206
    x = F.quantitative('X')
 
207
    f = F.formula(x)
 
208
 
 
209
    c = coxph(subjects, f)
 
210
 
 
211
    c.cache()
 
212
#    c.newton([0.4])