Full notebook is here. The problem is in the last Cox model at the end. The rest agree with the paper.
Background. W is a shared frailty. I have 430 districts that are in 48 states. I want the value for each state to be the same, so I just slice into the 48 values in W.value to produce the 430 (shared) frailties.
Can anyone spot my error? If I omit the W.value[stfips] part then the trace updates. E.g.,
fitted = np.exp(np.dot(X, b1))# + W.value[stfips])
I'm printing beta.value and it's certainly changing. I use the W.value[stfips] notation in the Weibull model with shared frailties, and it works fine. The Cox model is correct AFAIK and the same as in this question unless there is some error I'm missing.
I'm using PyMC 2.3.2.
def cox_spatialshared_model():
(obs_t, t, pscenter, hhcenter, ncomact, rleader,
dleader, inter1, inter2, fail, adj) = load_data_cox()
T = len(t) - 1
nsubj = len(obs_t)
stfips = (np.array(
[41, 41, 41, 41, 41, 41, 41, 45, 45, 45, 45, 35, 35, 35, 35,
35, 35, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
24, 24, 24, 24, 24, 24, 24, 24, 24, 30, 30, 30, 30, 30, 30,
18, 18, 18, 18, 18, 18, 27, 47, 47, 47, 47, 47, 47, 47, 47,
47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47,
43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 13, 13, 13, 13,
13, 8, 8, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26,
26, 26, 26, 26, 26, 26, 26, 21, 21, 21, 21, 21, 21, 21, 21,
21, 21, 32, 32, 32, 32, 31, 31, 31, 31, 31, 31, 46, 46, 46,
46, 46, 46, 46, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 29,
29, 29, 29, 29, 29, 29, 29, 3, 3 , 48, 48, 48, 48, 48, 48,
48, 48, 48, 48, 48, 48, 48, 48, 48, 10, 10, 10, 10, 10, 10,
10, 10, 34, 34, 34, 34, 34, 34, 34, 34, 34, 42, 42, 42, 42,
42, 2, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 4, 15,
15, 15, 12, 12, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
20, 20, 40, 40, 40, 22, 22, 16, 16, 16, 16, 16, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
16, 16, 16, 16, 16, 16, 16, 16, 25, 25, 25, 25, 25, 25, 25,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 36, 36, 36,
36, 36, 36, 11, 11, 11, 11, 11, 17, 17, 17, 17, 17, 17, 17,
17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 19,
19, 44, 44, 44, 44, 44, 44, 5, 38, 38, 38, 38, 38, 38, 38,
38, 38, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39,
39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39,
39, 39, 23, 23, 23, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33,
33, 9, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 7, 7,
28, 28, 28, 6]) - 1).tolist() # zero-based indexing
# neighbor weights for each state
ww = np.zeros((48, 48))
for i, row in enumerate(adj):
ww[i, row] = 1.
# make row-stochastic / row sums to 1
ww /= ww.sum(1)[:,None]
# risk set equals one if obs_t >= t
Y = np.array([[int(obs >= time) for time in t] for obs in obs_t])
# counting process. jump = 1 if obs_t \in [t[j], t[j+1])
dN = np.array([[Y[i, j]*(t[j + 1] >= obs_t[i]) * fail[i]
for i in range(nsubj)] for j in range(T)])
c = Gamma('c', .0001, .00001, value=.01)
r = Gamma('r', .001, .0001, value=.01)
dL0_star = r*np.array([t[j+1] - t[j] for j in range(len(t) - 1)])
mu = dL0_star * c # prior mean hazard
dL0 = Gamma('dL0', mu, c, value=np.ones(len(t)-1))
beta = Normal('beta', np.zeros(7), np.ones(7)*.00001,
value=np.array([-.36, -.26, -.29, -.22, -.61, -9.73, -.23]))
tau = Gamma('tau', alpha=.01, beta=.01, value=.001)
sigma = 1./(tau**.5)
theta = 1./tau
@stochastic
def car(t=tau, value=np.zeros(48)):
# the conditional mean is the average of the neighbors random effects
mu_ = np.dot(ww, value)
# scale precision for the number of neighbors
taux = t*(ww > 0).sum(1)
return normal_like(value, mu_, taux)
# zero-centered W
W = Lambda('W', lambda mu=car : mu - np.mean(mu))
X = np.column_stack((pscenter, hhcenter, ncomact, rleader, dleader, inter1, inter2))
@deterministic
def idt(b1=beta, dl0=dL0):
print beta.value
fitted = np.exp(np.dot(X, b1) + W.value[stfips])
yy = (Y[:,:T] * np.outer(fitted, dl0))
return np.transpose(yy)
dn_like = Poisson('dn_like', idt, value=dN, observed=True)
return locals()
variables = cox_spatialshared_model()
m3 = MCMC(variables, db='pickle', dbname='junk.pkl')
m3.use_step_method(AdaptiveMetropolis, m3.beta)
m3.sample(10)
Running this I get something like the following printed out for beta.value
[-0.36 -0.26 -0.29 -0.22 -0.61 -9.73 -0.23]
[-0.4373514 -0.29374146 -0.22849882 -0.26092336 -0.51110185 -9.71171661
-0.26297276]
[-0.46275125 -0.1619391 -0.25272671 -0.22495446 -0.51698795 -9.93903304
-0.18825614]
[-0.37534779 -0.23732345 -0.28472978 -0.18660755 -0.62196907 -9.59429074
-0.13662035]
[-0.43710784 -0.28943668 -0.4024984 -0.21735101 -0.46731924 -9.47100265
-0.18356244]
[ -0.36404428 -0.33324418 -0.33465088 -0.30341687 -0.64940603
-10.11364253 -0.23847174]
[-0.43687098 -0.23435138 -0.25747498 -0.29602854 -0.68939179 -9.91790107
-0.21975309]
[ -0.37618078 -0.25929464 -0.24237221 -0.23415328 -0.63732179
-10.51364586 -0.26740793]
[-0.36737311 -0.34901568 -0.25329154 -0.17279257 -0.56187575 -9.45392942
-0.18361278]
[-0.43960484 -0.28093066 -0.22398154 -0.18963571 -0.6667351 -9.45784103
-0.23457689]
[-0.32375244 -0.32216695 -0.40709281 -0.26659056 -0.58418152 -9.93182569
-0.19932304]
[-----------------100%-----------------] 10 of 10 complete in 0.1 sec
But for
m3.beta.trace()
array([[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23],
[-0.36, -0.26, -0.29, -0.22, -0.61, -9.73, -0.23]])
For posterity. Thanks to Chris Fonnesbeck for pointing out that the problem was that I did not give W as an argument to idt. This function should be
@deterministic
def idt(b1=beta, dl0=dL0, W=W):
print beta.value
fitted = np.exp(np.dot(X, b1) + W[stfips])
yy = (Y[:,:T] * np.outer(fitted, dl0))
return np.transpose(yy)
And everything should work as expected.