0

我无法使用 Rjags 和 runjags复制此纸模型http://dx.doi.org/10.1080/00031305.2017.1328374 。这是我正在使用的模型

library(rjags)
library(runjags)
modelo= "
  model{
for(i in 1:n){
#
# Likelihood part
#
censored[i] ~ dinterval(survival[i], censoring_limits[i])
# To deal with censoring:
# censored equals 1 for Benedict XVI as of 25/12/2016,
# 0 for the other popes.
# and the survival times of the other popes.
# censoring_limits equals 11.7 for Benedict XVI,
# and values (e.g. 32) greater than or equal to
# the survival times of the other popes.
survival[i] ~ dweib(r, mu[i]) # Basic Weibull assumption
mu[i] <- exp(beta[i]) # Defining beta as log(mu)
beta[i] <- beta_0 + beta_1*x_1[i] + beta_2*x_2[i]
# beta = log(mu) is a linear function of the covariates
}
#
####################################################
#
# Priors
#
beta_0 ~ dnorm(0.0, 1.0e-4) # Prior on beta_0 is normal with low precision
beta_1 ~ dnorm(0.0, 1.0e-4) # Prior on beta_1 is normal with low precision
beta_2 ~ dnorm(0.0, 1.0e-4) # Prior on beta_2 is normal with low precision
#
r ~ dexp(0.001) # Prior on r
#
####################################################
#
# Define the alphas
#
alpha_0 <- - beta_0 / r
alpha_1 <- - beta_1 / r
alpha_2 <- - beta_2 / r
#
# Percentage increases
#
percentage_increase_age <- 100*(exp(alpha_1) - 1)
percentage_increase_year <- 100*(exp(alpha_2) - 1)
#
####################################################
#
# Posterior median at various covariate values
#
beta_med <- beta_0 + beta_1*x_1_new + beta_2*x_2_new
# New values need to be supplied
t_median <- pow(log(2) * exp(-beta_med), 1 / r)
#
####################################################
#
# Predictive distribution of age at the new values
#
beta_Francis <- beta_0 + beta_1*age_Francis + beta_2*year_Francis
# Values of age_Francis and year_Francis need to be provided
mu_Francis <- exp(beta_Francis)
survival_Francis~ dweib(r, mu_Francis)T(present_length, upper_length)
# Take into account the current pontificate length
# Also specify a sensible upper bound
age_Francis_predictive <- age_at_election + survival_Francis
# Work also with age
}"
writeLines(modelo, con = "TEMPmodel.txt")

我用来提供模型的数据库位于文章链接的补充部分（在顶部）。这是我将数据提供给模型的方式

writeLines(modelo, con = "TEMPmodel.txt")
base1<-read.csv("C:/Users/jgush/Downloads/popes_25_December_2016.csv")
base1<-base1[-c(1),]
#base1$Survival[1]<-NA
base1$Censoring_limits<-ifelse(base1$Censored==1,11.7,32)
datalist<-list("survival"=base1$Survival,"censoring_limits"=base1$Censoring_limits,"n"=length(base1$Year.Elected),
               "censored"=base1$Censored,"age_at_election"=76.29,"age_Francis"=80.1,
               "year_Francis"=2013,"x_1_new"=75,"x_2_new"=1666,"present_length"=3.79,"upper_length"=17,
               "x_1"=base1$Age.Election,"x_2"=base1$Year.Elected)
modelo.con.jags <- jags.model(file = "TEMPmodel.txt", 
                              data = datalist,
                              n.adapt = 500) # Burning-period

而且我总是收到下一个错误“update.jags（object，n.iter，...）中的错误：节点生存错误[1]无法计算日志密度”我迷路了，我不知道为什么发生错误任何帮助将不胜感激。

4

0 回答 0