teratail header banner
teratail header banner
質問するログイン新規登録

質問編集履歴

1

内容

2016/06/26 04:30

投稿

kkoorroo
kkoorroo

スコア8

title CHANGED
File without changes
body CHANGED
@@ -12,9 +12,7 @@
12
12
 
13
13
  ###該当のソースコード
14
14
  ```R
15
- #### R2Winbugs Code for a Bayesian hierarchical test of the Geurts data.
15
+
16
- #### Do children with ADHD perform worse on the Wisconsin Card Sorting Task
17
- #### than children who develop typically (i.e., normal controls)?
18
16
 
19
17
  rm(list=ls()) # clears workspace
20
18
 
@@ -23,108 +21,7 @@
23
21
  library(R2WinBUGS)
24
22
  bugsdir <- "C:/Program Files/WinBUGS14"
25
23
 
26
- ### Geurts data:
27
- # Normal Controls:
28
- num.errors <- c(15,10, 61,11, 60, 44, 63, 70, 57,11, 67, 21, 89,12, 63, 11, 96,10, 37,19, 44,18, 78, 27, 60,14)
29
- nc <- c(89,74,128,87,128,121,128,128,128,78,128,106,128,83,128,100,128,73,128,86,128,86,128,100,128,79)
30
- kc <- nc - num.errors
31
- nsc <- length(kc)
32
- # ADHD:
33
- num.errors <- c(88, 50, 58,17, 40, 18,21, 50, 21, 69,19, 29,11, 76, 46, 36, 37, 72,27, 92,13, 39, 53, 31, 49,
34
- 57,17,10,12,21, 39, 43, 49,17, 39,13, 68, 24, 21,27, 48, 54, 41, 75, 38, 76,21, 41, 61,24, 28,21)
35
- na <- c(128,128,128,86,128,117,89,128,110,128,93,107,87,128,128,113,128,128,98,128,93,116,128,116,128,
36
- 128,93,86,86,96,128,128,128,86,128,78,128,111,100,95,128,128,128,128,128,128,98,127,128,93,110,96)
37
- ka <- na - num.errors
38
- nsa <- length(ka)
39
-
40
- # two-sided p-value = .72
41
- t.test(kc/nc, ka/na, alternative = c("two.sided"), paired=F)
42
24
 
43
- data <- list("nc","kc","nsc","na","ka","nsa") # to be passed on to WinBUGS
44
-
45
- myinits <- list(
46
- list(mu = 0, sigma = 1, delta = 0),
47
- list(mu = -.8, sigma = 2, delta = -.5),
48
- list(mu = .8, sigma = 1.5, delta = .5))
49
-
50
- parameters <- c("delta")
51
-
52
- # The following command calls WinBUGS with specific options.
53
- # For a detailed description see Sturtz, Ligges, & Gelman (2005).
54
- samples <- bugs(data, inits=myinits, parameters,
55
- model.file ="Geurts.txt",
56
- n.chains=3, n.iter=100000, n.burnin=1000, n.thin=1,
57
- DIC=T, bugs.directory=bugsdir, #DIC=T or else error message
58
- codaPkg=F, debug=T)
59
- # Now the values for delta are in the "samples" object, ready for inspection.
60
-
61
- # Collect posterior samples across all chains:
62
- delta.posterior <- samples$sims.list$delta
63
-
64
- #============ Unrestricted analysis ===========================
65
- #============ BFs based on logspline fit ===========================
66
- library(polspline) # this package can be installed from within R
67
- fit.posterior <- logspline(delta.posterior)
68
-
69
- # 95% confidence interval:
70
- x0 <- qlogspline(0.025,fit.posterior)
71
- x1 <- qlogspline(0.975,fit.posterior)
72
-
73
- posterior <- dlogspline(0, fit.posterior) # this gives the pdf at point delta = 0
74
- prior <- dnorm(0) # height of order-restricted prior at delta = 0
75
- BF01 <- posterior/prior
76
- # BF01 = 3.96, BF10 = 0.25
77
-
78
- #####################################################################
79
- ### Plot Prior and Posterior
80
- #####################################################################
81
- par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
82
- font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
83
- Nbreaks <- 80
84
- y <- hist(delta.posterior, Nbreaks, plot=F)
85
- plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="S", lwd=2, lty=2,
86
- xlim=c(-3,3), ylim=c(0,2), xlab=" ", ylab="Density", axes=F)
87
- axis(1, at = c(-3,-2,-1,0,1,2,3), lab=c("-3","-2", "-1", "0", "1", "2", "3"))
88
- axis(2)
89
- mtext(expression(delta), side=1, line = 2.8, cex=2)
90
- #now bring in log spline density estimation:
91
- par(new=T)
92
- plot(fit.posterior, ylim=c(0,2), xlim=c(-3,3), lty=1, lwd=1, axes=F)
93
- points(0, dlogspline(0, fit.posterior),pch=19, cex=2)
94
- # plot the prior:
95
- par(new=T)
96
- plot ( function( x ) dnorm( x, 0, 1 ), -3, 3, xlim=c(-3,3), ylim=c(0,2), lwd=2, lty=1, ylab=" ", xlab = " ")
97
- points(0, dnorm(0), pch=19, cex=2)
98
-
99
- text(0.3, 1.5, labels = "Posterior", cex = 1.5, pos=4)
100
- text(1, 0.3, labels = "Prior", cex=1.5, pos=4)
101
-
102
-
103
- Geurts.txtは以下の内容です(これは自分で書きました).
104
- # Geurts
105
- model{
106
- for (i in 1:nsc){
107
- kc[i] ~ dbin(thetac[i], nc[i])
108
- thetac[i] <- phi(phic[i])
109
- phic[i] ~ dnorm(muc, lambda)
110
- }
111
- for (j in 1:nsa){
112
- ka[j] ~ dbin(thetaa[j], na[j])
113
- thetaa[j] <- phi(phia[j])
114
- phia[j] ~ dnorm(mua, lambda)
115
- }
116
- muc <- mu+alpha/2
117
- mua <- mu-alpha/2
118
- # Priors
119
- mu ~ dnorm(0,1)
120
- sigma ~ dunif(0,10)
121
- alpha <- delta*sigma
122
- lambda <- pow(sigma,-2)
123
- delta ~dnorm(0,1)
124
- # sampling from Prior Distribution for Delta
125
- deltaprior ~ dnorm(0,1)
126
- }
127
-
128
25
  ```
129
26
 
130
27
  ###試したこと