[Case Study]Hierarchical Multinomial Models using WinBUGS

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

1. Dirichlet / Gamma model...
   
2. model {
   
3. for (i in 1:npupil) {
   
4. Goals ~ dcat(p[School,])
   
5. }
   
6. for (j in 1:nschool) {
   
7. for (k in 1:3) {
   
8. p[j,k] <- q[j,k]/sum(q[j,])
   
9. q[j,k] ~ dgamma(a[k], 1)
   
10. }
    
11. }
    
12. for (k in 1:3) {
    
13. a[k] ~ dgamma(1, 0.001)
    
14. p.pop[k] <- a[k]/sum(a[]) # population mean of p[,k]
    
15. }
    
16. dummy <- Gender[1] # stop WinBUGS complaining about unused variable
    
17. }

复制代码

Data:
Click arrow for data - same data used for all models in this example

Inits:
list(a=c(1,1,1))

Poor convergence. Click arrow for trace plots of a[], which are highly cross-correlated.

No difference between schools...

1. model {
   
2. for (i in 1:npupil) {
   
3. Goals ~ dcat(p[i,])
   
4. for (k in 1:3) {
   
5. p[i,k] <- q[i,k]/sum(q[i,])
   
6. log(q[i,k]) <- a[i,k]
   
7. }
   
8. a[i,1] <- b[1] + b.boy*Gender
   
9. a[i,2] <- b[2]
   
10. a[i,3] <- 0
    
11. }
    
12. b[1] ~ dnorm(0, 0.0001)
    
13. b[2] ~ dnorm(0, 0.0001)
    
14. b.boy ~ dnorm(0, 0.0001)
    
15. or.boy <- exp(b.boy)
    
16. 
    
17. qboy[1] <- exp(b[1] + b.boy); qboy[2] <- exp(b[2]); qboy[3] <- 1
    
18. qgirl[1] <- exp(b[1]); qgirl[2] <- exp(b[2]); qgirl[3] <- 1
    
19. 
    
20. # Probabilities of preferring 1) sports 2) popularity 3) grades for boys and girls separately
    
21. for (k in 1:3) {
    
22.    p.boy[k] <- qboy[k]/sum(qboy[])
    
23.    p.girl[k] <- qgirl[k]/sum(qgirl[])
    
24. }
    
25. dummy <- School[1] + nschool
    
26. }

复制代码

Inits:
list(b.boy=0, b=c(0,0))

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   b.boy   0.9855   0.2476   0.005047   0.5125   0.9823   1.473   4003   50000
   or.boy   2.763   0.6999   0.01422   1.669   2.671   4.364   4003   50000
   
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
   Dbar   Dhat   pD   DIC   
Goals   957.489   954.473   3.016   960.506   
total   957.489   954.473   3.016   960.506   
   

1. Independent school effects...
   
2. model {
   
3. for (i in 1:npupil) {
   
4. Goals ~ dcat(p[i,])
   
5. for (k in 1:3) {
   
6. p[i,k] <- q[i,k]/sum(q[i,])
   
7. log(q[i,k]) <- a[i,k]
   
8. }
   
9. a[i,1] <- b[School, 1] + b.boy*Gender
   
10. a[i,2] <- b[School, 2]
    
11. a[i,3] <- 0
    
12. }
    
13. for (j in 1:nschool) {
    
14. b[j,1] ~ dnorm(0, 0.0001)
    
15. b[j,2] ~ dnorm(0, 0.0001)
    
16. b[j,3] <- 0
    
17. for (k in 1:3) {   
    
18. egirl[j,k] <- exp(b[j,k])
    
19. }
    
20. eboy[j,1] <- exp(b[j,1] + b.boy); for (k in 2:3) {eboy[j,k] <- exp(b[j,k])}
    
21. for (k in 1:3) {
    
22. p.girl[j,k] <- egirl[j,k] / sum(egirl[j,])
    
23. p.boy[j,k] <- eboy[j,k] / sum(eboy[j,])
    
24. }
    
25. }
    
26. b.boy ~ dnorm(0, 0.0001)
    
27. or.boy <- exp(b.boy)
    
28. }

复制代码

Inits:
list(b.boy=0, b=structure(.Data=c(0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA),.Dim=c(9,3)))

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   b.boy   1.088   0.2639   0.005174   0.5781   1.084   1.61   4002   50000
   or.boy   3.074   0.8356   0.01662   1.783   2.956   5.005   4002   50000

   Dbar   Dhat   pD   DIC   
Goals   936.567   917.409   19.158   955.725   
total   936.567   917.409   19.158   955.725   

1. Hierarchical school effects...
   
2. model {
   
3. for (i in 1:npupil) {
   
4. Goals ~ dcat(p[i,])
   
5. for (k in 1:3) {
   
6. p[i,k] <- q[i,k]/sum(q[i,])
   
7. log(q[i,k]) <- a[i,k]
   
8. }
   
9. a[i,1] <- b[School, 1] + b.boy*Gender
   
10. a[i,2] <- b[School, 2]
    
11. a[i,3] <- 0
    
12. }
    
13. for (j in 1:nschool) {
    
14. b[j,1] ~ dnorm(mu[1], tau[1])
    
15. mub2[j] <- mu[2] + s[2]/s[1]*cor*(b[j,1] - mu[1])
    
16. b[j,2] ~ dnorm(mub2[j], taub2)
    
17. b[j,3] <- 0
    
18. for (k in 1:3) {   
    
19. egirl[j,k] <- exp(b[j,k])
    
20. }
    
21. eboy[j,1] <- exp(b[j,1] + b.boy); for (k in 2:3) { eboy[j,k] <- exp(b[j,k]) }
    
22. for (k in 1:3) {
    
23. p.girl[j,k] <- egirl[j,k] / sum(egirl[j,])
    
24. p.boy[j,k] <- eboy[j,k] / sum(eboy[j,])
    
25. }
    
26. }
    
27. vb2 <- (1 - cor*cor)*v[2]
    
28. tau[1] <- 1/v[1]
    
29. taub2 <- 1/vb2
    
30. for (k in 1:2) {
    
31. mu[k] ~ dnorm(0, 0.0001)
    
32. v[k] <- s[k]*s[k]
    
33. s[k] ~ dunif(0, 100)
    
34. }
    
35. cor ~ dunif(0, 1)
    
36. b.boy ~ dnorm(0, 0.0001)
    
37. or.boy <- exp(b.boy)
    
38. }

复制代码

Inits:
list(mu=c(0,0), s=c(1,1), cor=0.5)

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   b.boy   1.053   0.2549   0.003429   0.5596   1.052   1.56   4002   200000
   or.boy   2.961   0.7715   0.01031   1.75   2.862   4.758   4002   200000

   Dbar   Dhat   pD   DIC   
Goals   937.938   924.479   13.459   951.397   
total   937.938   924.479   13.459   951.397

Monitor p.girl and p.boy to obtain the data for the graph in Figure 10.5, which was produced in R.

扫码加我 拉你入群

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

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

分享0 收藏2 回帖

关键词：Hierarchical Multinomial Case study winbugs models

辛苦了，楼主