質問編集履歴
1
内容
title
CHANGED
File without changes
|
body
CHANGED
@@ -12,9 +12,7 @@
|
|
12
12
|
|
13
13
|
###該当のソースコード
|
14
14
|
```R
|
15
|
-
|
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
|
###試したこと
|