Я не могу воспроизвести результаты предоставленного кода BUGS с помощью PyMC. Модель BUGS представляет собой модель мультипликативной интенсивности Андерсена-Гилла Кокса-PH.
model
{
# Set up data
for(i in 1:Nsubj) {
for(j in 1:T) {
# risk set = 1 if obs.t >= t
Y[i,j] <- step(obs.t[i] - t[j] + eps)
# counting process jump = 1 if obs.t in [ t[j], t[j+1] )
# i.e. if t[j] <= obs.t < t[j+1]
dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * FAIL[i]
}
Useless[i] <- pscenter[i] + hhcenter[i] + ncomact[i]
+ rleader[i] + dleader[i] + inter1[i] + inter2[i]
}
# Model
for(j in 1:T) {
for(i in 1:Nsubj) {
dN[i, j] ~ dpois(Idt[i, j]) # Likelihood
Idt[i, j] <- Y[i, j] * exp(beta[1]*pscenter[i] + beta[2]*
hhcenter[i] + beta[3]*ncomact[i] + beta[4]*rleader[i] + beta[5]*dleader[i] + beta[6]*inter1[i] + beta[7]*inter2[i]) * dL0[j] # Intensity
}
dL0[j] ~ dgamma(mu[j], c)
mu[j] <- dL0.star[j] * c # prior mean hazard
}
c ~ dgamma(0.0001, 0.00001)
r ~ dgamma(0.001, 0.0001)
for (j in 1 : T) { dL0.star[j] <- r * (t[j + 1] - t[j]) }
# next line indicates number of covariates and is for the corresponding betas
for(i in 1:7) {beta[i] ~ dnorm(0.0,0.00001)}
}
Я использую следующие начальные значения
list(beta=c(-.36,-.26,-.29,-.22,-.61,-9.73,-.23), c=0.01, r=0.01, dL0=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
Я использую одну цепочку (пока) и 5000 итераций для прожига. Я запускаю оценку для 10 000 дополнительных итераций и получаю те же точечные оценки, что и в статье. Они также близки к более ранним частотным оценкам.
OpenBUGS> samplesStats('beta')
mean sd MC_error val2.5pc median val97.5pcstart sample
beta[1] 3.466 0.8906 0.03592 1.696 3.48 5.175 501 9500
beta[2] -0.04155 0.06253 0.002487 -0.1609 -0.04355 0.08464 501 9500
beta[3] -0.009709 0.07353 0.002008 -0.1544 -0.01052 0.1365 501 9500
beta[4] 0.3535 0.1788 0.004184 -0.01523 0.3636 0.6724 501 9500
beta[5] 0.08454 0.1652 0.004261 -0.2464 0.08795 0.3964 501 9500
beta[6] -4.109 1.325 0.05224 -6.617 -4.132 -1.479 501 9500
beta[7] 0.1413 0.08594 0.003381 -0.03404 0.1423 0.3031 501 9500
OpenBUGS> samplesStats('c')
mean sd MC_error val2.5pc median val97.5pcstart sample
c 4.053 1.08 0.02896 2.202 3.974 6.373 1001 10000
OpenBUGS> samplesStats('r')
mean sd MC_error val2.5pc median val97.5pcstart sample
r 0.01162 0.002929 7.846E-5 0.007387 0.01119 0.01848 1001 10000
Я попытался воспроизвести это в PyMC 2.3.2 со следующим кодом. Полный код репликации доступен здесь
def cox_model(dta):
(t, obs_t, pscenter, hhcenter, ncomact, rleader,
dleader, inter1, inter2, fail) = load_data_cox()
T = len(t) - 1
nsubj = len(obs_t)
# 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]*int(t[j+1] >= obs_t[i])*fail[i] for i in range(nsubj)] for j in range(T)])
c = Gamma('c', .0001, .00001, value=.1)
r = Gamma('r', .001, .0001, value=.1)
dL0_star = r*np.array([t[j+1] - t[j] for j in range(T)])
mu = dL0_star * c # prior mean hazard
dL0 = Gamma('dL0', mu, c, value=np.ones(T))
beta = Normal('beta', np.zeros(7), np.ones(7)*.00001,
value=np.array([-.36, -.26, -.29, -.22, -.61, -9.73, -.23]))
@deterministic
def idt(b1=beta, dl0=dL0):
mu_ = [[Y[i,j] * np.exp(b1[0]*pscenter[i] + b1[1]*hhcenter[i] +
b1[2]*ncomact[i] + b1[3]*rleader[i] +
b1[4]*dleader[i] + b1[5]*inter1[i] +
b1[6]*inter2[i])*dl0[j] for i in range(nsubj)]
for j in range(T)] # intensity
return mu_
dn_like = Poisson('dn_like', idt, value=dN, observed=True)
return locals()
m = MCMC(cox_model())
m.sample(15000)
Впрочем, я и близко не подхожу к одинаковым точечным оценкам. я получаю что-то вроде
beta:
Mean SD MC Error 95% HPD interval
------------------------------------------------------------------
-0.537 1.094 0.099 [-2.549 1.492]
0.276 0.048 0.004 [ 0.184 0.36 ]
-1.092 0.385 0.038 [-1.559 -0.371]
-1.461 0.746 0.073 [-2.986 -0.496]
-1.865 0.382 0.038 [-2.471 -1.329]
3.778 1.539 0.133 [ 1.088 6.623]
-0.449 0.109 0.01 [-0.661 -0.26 ]
Posterior quantiles:
2.5 25 50 75 97.5
|---------------|===============|===============|---------------|
-2.892 -1.274 -0.385 0.268 1.253
0.191 0.244 0.278 0.305 0.374
-1.553 -1.434 -1.179 -0.793 -0.258
-3.132 -1.856 -1.196 -0.904 -0.526
-2.471 -2.199 -1.864 -1.632 -1.201
1.287 2.685 3.601 4.72 7.262
-0.714 -0.519 -0.445 -0.368 -0.273
Самое тревожное, что знаки разные. Я подумал, может быть, это просто проблема сходимости, поэтому я запустил его за ночь с 50 000 итераций без особых изменений. Может быть, в моей модели PyMC есть какая-то ошибка или разница, особенно со спецификацией dL0?
Я пробовал с разными начальными значениями. Я пытался запустить модель для многих итераций. Я сосредоточил априорные оценки на точечных оценках из BUGS.
