完整的笔记本在这里。问题出在最后的 Cox 模型中。其余的同意这篇论文。
背景。W 是一个共同的弱点。我有 48 个州的 430 个区。我希望每个状态的值都相同，所以我只需将 48 个值分割成W.value430 个（共享）弱点。
谁能发现我的错误？如果我省略了该W.value[stfips]部分，那么跟踪会更新。例如，
fitted = np.exp(np.dot(X, b1))# + W.value[stfips])
我正在打印beta.value，它肯定会发生变化。我W.value[stfips]在 Weibull 模型中使用具有共同弱点的符号，它工作正常。Cox 模型是正确的 AFAIK，并且与这个问题相同，除非我遗漏了一些错误。
我正在使用 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)
运行这个我得到类似下面打印出来的 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
但对于
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]])