I am running into issues with the following model. The first was to do with the zeros trick that works in BUGS but not nimble. I *think* I get around that by reading in a vector of zeros per another thread response. The second issue I run in to is with the indexing of my spatial random effect. It's a simple uncorrelated random effect in this model but I get the errors below. Please advise.
#AFT attempt nimble
library(nimble)
ModelCode <- nimbleCode({
for (i in 1:I){ #individuals
temp[i]<-(log(timeF[i]+0.01)-log(lam[i]))/sigma
log(lam[i])<-beta[1]+beta[2]*x1[i]+beta[3]*x2[i]+beta[4]*x3[i]+beta[5]*x4[i]+
beta[6]*x5[i]+beta[7]*x6[i]+beta[8]*x7[i]+beta[9]*x8[i]+
beta[10]*x9[i]+beta[11]*x10[i]+beta[12]*x11[i]+beta[13]*x12[i]+
beta[14]*x13[i]+beta[15]*x14[i]+beta[16]*x15[i]+
u[spt[i]]
#logistic distn surv and density
s[i]<-1/(1+exp(temp[i]))
f[i]<-exp(temp[i])*pow(s[i],2)
#log like
L[i]<-indic[i]*log(f[i]/(sigma*(timeF[i]+0.01)))+(1-indic[i])*log(s[i])
#Poisson zeros trick
# zeros[i]<-0
new[i]<- -L[i]
zeros[i]~dpois(new[i])
}#close i loop
for (j in 1:16){beta[j]~dnorm(0,0.001)}
tauu<-pow(sdu,-2)
sdu~dunif(0,10)
sigma~dunif(0.01,10)
for (j in 1:64){u[j]~dnorm(0,tauu)}#close j loop
})
#data
databrca=read.csv("C:\\Users\\carrollr\\OneDrive - UNC-Wilmington\\Student Projects\\Spring 19\\Data\\BrCadat.csv")
AFTdata=list(timeF=databrca$timeF,indic=databrca$indic,
x1=databrca$x1,x2=databrca$x2,x3=databrca$x3,
x4=databrca$x4,x5=databrca$x5,x6=databrca$x6,
x7=databrca$x7,x8=databrca$x8,x9=databrca$x9,
x10=databrca$x10,x11=databrca$x11,x12=databrca$x12,
x13=databrca$x13,x14=databrca$x14,x15=databrca$x15,#)
zeros=rep(0,length(databrca$timeF)))
#predictor order: x1=agedx,x2=race (AA), x3=grade high,x4=igrade unkn
#x5=surg yes,x6=surg unkn, x7=married,x8=prev married,
#x9=marital stat unkn,x10=eppr -+,x11=erpr +-,x12=erpr ++,
#x13=erpr unkn,x14=radition yes, x15=radiation unkn
init=list(beta=rep(0,16),sdu=5,u=rep(0,64),sigma=1)
## build the model
rModel <- nimbleModel(ModelCode, constants = list(I = length(databrca$timeF),spt=databrca$spt),
data <- AFTdata, check = FALSE )
#run model
mcmc.out <- nimbleMCMC(code=ModelCode, constants = list(I = I),
data <- AFTdata,
inits = init,
nchains = 2, niter = 5000, summary = T, WAIC = T)
names(mcmc.out)
mcmc.out$summary$all.chains#parameter estimates
mcmc.out$WAIC
...
Warning: dynamic index out of
bounds: exp(beta[1] + beta[2] * x1[i] + beta[3] * x2[i] + beta[4] * x3[i] +
beta[5] * x4[i] + beta[6] * x5[i] + beta[7] * x6[i] + beta[8] * x7[i] + beta[9]
* x8[i] + beta[10] * x9[i] + beta[11] * x10[i] + beta[12] * x11[i] + beta[13] *
x12[i] + beta[14] * x13[i] + beta[15] * x14[i] + beta[16] * x15[i] + u[spt[i]])
NAs were detected in model
variable: spt.
NaNs were detected in model
variables: lam, temp, s, f, L, new, logProb_zeros.
model building finished.
Warning message:
In warnRHSonlyDynIdx() :
Detected use of non-constant indexes: spt[1],
spt[2], spt[3], spt[4], spt[5], spt[6], spt[7], spt[8], spt[9], spt[10],
spt[11], spt[12], spt[13], spt[14], spt[15], spt[16], spt[17], spt[18],
spt[19], spt[20], spt[21], spt[22], spt[23], spt[24], spt[25], spt[26],
spt[27], spt[28], spt[29], spt[30], spt[31], spt[32], spt[33], spt[34],
spt[35], spt[36], spt[37], spt[38], spt[39], spt[40], spt[41], spt[42],
spt[43], spt[44], spt[45], spt[46], spt[47], spt[48], spt[49], spt[50], ... For
computational efficiency we recommend specifying these in 'constants'.
> spt=databrca$spt
> spt[1]
[1] 60
> summary(spt)
Min. 1st Qu.
Median Mean 3rd Qu. Max.
1.00
25.00 38.00 35.96
48.00 64.00
> is.integer(spt)
[1] TRUE
