[本週專題]WinBUGS: Undefined Real Result?

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Dear lister
I have run a zero-inflated Poisson model in WinBUGS without problems, and now I am trying to run its equivalent negative binomial. However, I get an "undefined real result" trap message over and over. Here is the data set http://www.megafileupload.com/en/file/474809/data1.html. The trap arises from the Pr(Y > 0) part of the model since the Pr(Y = 0) runs without problems. I've tried different inits and priors to no avail.
Thanks for your help!

1. #Data handling
   
2. library(R2WinBUGS)
   
3. setwd('') # working directory
   
4. dget('data1')
   
5. data1 <- data1[order(data1$C), ]
   
6.     n = nrow(data1)
   
7.     n0 = length(which(data1$C == 0))
   
8. 
   
9. sink("ZINB.txt")
   
10. cat("
    
11. model {
    
12. 
    
13. # Pr(Y=0|x,z) # probability of zeros
    
14. K <- 10000
    
15. for(i in 1:n0){
    
16. zeros0 <- 0
    
17. zeros0 ~ dpois(mu0)
    
18. mu0 <- -log(p00 + (1-p00)*pow(overd, r)) + K
    
19. 
    
20. overd <- r/(r + lambda0)
    
21. p00 <- min( 0.999999, max(.000001,p0) ) # avoid stackoverflow
    
22. logit(p0) <- a0 + b0.A*A + b0.M*M + b0.cat*cat + b0.d*d
    
23. log(lambda0)<- a + logh + b.A*A + b.M*M + b.cebo*cebo + b.descart*descart + b.d*d
    
24. }
    
25. 
    
26. # Pr(Y>0|x,z) # probability of positive counts
    
27. for(j in (n0+1):n){
    
28. zeros[j] <- 0
    
29. zeros[j] ~ dpois(mu[j])
    
30. mu[j] <- -log((1-p1[j])*fnb[j]) + K
    
31. 
    
32. lfnb[j] <-  r*log(overd[j]) + C[j]*log(1 - overd[j]) + loggam(C[j] + r) - loggam(r) - loggam(C[j] + 1) # - log(1 - pow(overd[j], r))
    
33. fnb[j] <- exp( lfnb[j] )
    
34. overd[j] <- r/(r + exp(lambda[j]))
    
35. p1[j] <- min( 0.999999, max(.000001,p[j]) )
    
36. logit(p[j]) <- a0 + b0.A*A[j] + b0.M*M[j] + b0.cat*cat[j] + b0.d*d[j] + b0.descart*descart[j]
    
37. log(lambda[j])<- a + logh[j] + b.A*A[j] + b.M*M[j] + b.cebo*cebo[j] + b.descart*descart[j] + b.d*d[j] + b.dcos*dcos[j]
    
38. pred[j] <- (1 - p1[j])*lambda[j] # predicted values (see Zuur et al 2009)
    
39. var[j] <- (1 - p1[j])*(lambda[j] + lambda[j]*lambda[j]/r) + lambda[j]*lambda[j]*(p1[j]*p1[j] + p1[j]) # variance (see Zuur et al 2009)
    
40. res[j] <- (C[j] - lambda[j]*(1 - p[j]))/sqrt( var[j] ) # Pearson residuals
    
41. disp[j] <- pred[j]/var[j] # dispersion
    
42. }
    
43. 
    
44. # Vague priors for model coefficients
    
45. r ~ dunif(0,10)
    
46. a0 ~ dnorm(0,0.001)
    
47. a ~ dnorm(0,0.001)
    
48. b.A ~ dnorm(0,0.001)
    
49. b0.A ~ dnorm(0,0.001)
    
50. b.M ~ dnorm(0,0.001)
    
51. b0.M ~ dnorm(0,0.001)
    
52. b.d ~ dnorm(0,0.001)
    
53. b0.d ~ dnorm(0,0.001)
    
54. b.descart ~ dnorm(.0,0.001)
    
55. b0.descart ~ dnorm(0,0.001)
    
56. b0.cat ~ dnorm(0,0.001)
    
57. b.dcos ~ dnorm(0,0.001)   
    
58. b.cebo ~ dnorm(0,0.001)
    
59. }
    
60. ",fill=TRUE)
    
61. sink()
    
62. 
    
63. # Data passed to WinBUGS
    
64. win.data <- list( C = data1$C, A = as.numeric(data1$A), M = as.numeric(data1$M), d = as.numeric(scale(data1$dcol)), dcos = as.numeric(scale(data1$dcos)), cebo = as.numeric(data1$cebo), descart = as.numeric(data1$descart), cat = as.numeric(data1$cat), logh = log(data1$h), n = n, n0 = n0 )
    
65. 
    
66. # Initial values
    
67. inits <- function(){ list( r = runif(1, 1, 5), a0 = rnorm(1), a = rnorm(1), b.A = rnorm(1), b0.A = rnorm(1), b.M = rnorm(1), b0.M = rnorm(1), b.d = rnorm(1), b0.d = rnorm(1), b.descart = rnorm(1), b0.descart = rnorm(1), b0.cat = rnorm(1), b.dcos = rnorm(1), b.cebo = rnorm(1) ) }
    
68. 
    
69. # Parameters
    
70. params <- c( "r", "mu", "lambda", "a", "a0", "b.A", "b0.A", "b.M", "b0.M", "b.d", "b0.d", "b.descart", "b0.descart", "b0.cat", "b.dcos", "b.cebo", "res", "pred", "var" )
    
71. 
    
72. # MCMC settings
    
73. nc <- 3
    
74. ni <- 30000
    
75. nb <- 3000
    
76. nt <- 3
    
77. 
    
78. # WinBUGS execution
    
79. out <- bugs( data = win.data, inits = inits, parameters.to.save = params, model.file = "ZINB.txt", n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = TRUE )
    
80. 
    
81. print(out, dig = 2)

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Undefined winbugs Defined DEFINE WINBUG undefined equivalent different problems message

Ok, I found a missing zero in one of the priors:

b.descart ~ dnorm(.0,0.001)  
Apparently:

BUGS doesn't take .0 as a 0.0, and
WinBUGS doesn't take .0 as a syntax error
Another point: the model neither runs with more than three covariates nor more than two categorical covariates. Any clue?

As an aside, the value I gave to the constant K makes the deviance and DIC to get wildly high values.