spatial AFT model issues

12 views
Skip to first unread message

Rachel Carroll

unread,
Jun 14, 2021, 7:35:06 PMJun 14
to nimble-users
Hello,

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.

Rachel

#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


Chris Paciorek

unread,
Jun 15, 2021, 8:52:46 PMJun 15
to Rachel Carroll, nimble-users
hi Rachel,

It looks like you are specifying constants correctly when you create the model but you then specify the constants again and don't provide 'spt' when you call nimbleMCMC with only the code and not the model. So the model is being recreated within nimbleMCMC using only "I" as constants. Please try either specifying the full set of constants when passing in the code to nimbleModel. Alternatively pass the model into nimbleMCMC and not the code, and don't specify constants in the call to nimbleMCMC, only in the call to nimbleModel.

Side note: In nimble you can use "dnorm(some_mean, sd = some_sd)" and not have to create the tauu intermediate.

Side note #2: We don't allow you to specify 'zeroes' on the left hand side twice (once to set as zero and once with dpois), hence you need to set the value of 'zeroes' as part of 'data' as you have done.

-chris

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/10d37ef7-29db-4dd5-876e-3e945bf89f88n%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages