Skip to content

Commit 244b19e

Browse files
martin-smiramartin-smira
committed
add MPT models
1 parent 3a6887d commit 244b19e

File tree

2 files changed

+333
-0
lines changed

2 files changed

+333
-0
lines changed

‎Bayesian_Cognitive_Modeling/CaseStudies/MPT/MPT_1_Stan.R‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
# clears workspace:
2+
rm(list=ls())
3+
4+
library(rstan)
5+
6+
model <- "
7+
// Multinomial Processing Tree
8+
data {
9+
int<lower=1> n;
10+
int<lower=0,upper=n> k[4];
11+
}
12+
parameters {
13+
real<lower=0,upper=1> c;
14+
real<lower=0,upper=1> r;
15+
real<lower=0,upper=1> u;
16+
}
17+
transformed parameters {
18+
simplex[4] theta;
19+
20+
// MPT Category Probabilities for Word Pairs
21+
theta[1] <- c * r;
22+
theta[2] <- (1 - c) * u ^ 2;
23+
theta[3] <- (1 - c) * 2 * u * (1 - u);
24+
theta[4] <- c * (1 - r) + (1 - c) * (1 - u) ^ 2;
25+
}
26+
model {
27+
// Priors
28+
c ~ beta(1, 1); // can be removed
29+
r ~ beta(1, 1); // can be removed
30+
u ~ beta(1, 1); // can be removed
31+
32+
// Data
33+
k ~ multinomial(theta);
34+
}"
35+
36+
k <- c(45, 24, 97, 254)
37+
n <- sum(k)
38+
data <- list(k=k, n=n) # To be passed on to Stan
39+
40+
myinits <- list(
41+
list(c=.5, r=.5, u=.5))
42+
43+
parameters <- c("c", "r", "u") # Parameters to be monitored
44+
45+
# The following command calls Stan with specific options.
46+
# For a detailed description type "?stan".
47+
samples_1 <- stan(model_code=model,
48+
data=data,
49+
init=myinits, # If not specified, gives random inits
50+
pars=parameters,
51+
iter=11000,
52+
chains=1,
53+
thin=1,
54+
warmup=1000, # Stands for burn-in; Default = iter/2
55+
# seed=123 # Setting seed; Default is random seed
56+
)
57+
58+
k <- c(106, 41, 107, 166)
59+
n <- sum(k)
60+
data <- list(k=k, n=n) # To be passed on to Stan
61+
62+
samples_2 <- stan(fit=samples_1,
63+
data=data,
64+
init=myinits, # If not specified, gives random inits
65+
pars=parameters,
66+
iter=11000,
67+
chains=1,
68+
thin=1,
69+
warmup=1000, # Stands for burn-in; Default = iter/2
70+
# seed=123 # Setting seed; Default is random seed
71+
)
72+
73+
k <- c(243, 64, 65, 48)
74+
n <- sum(k)
75+
data <- list(k=k, n=n) # To be passed on to Stan
76+
77+
samples_3 <- stan(fit=samples_1,
78+
data=data,
79+
init=myinits, # If not specified, gives random inits
80+
pars=parameters,
81+
iter=11000,
82+
chains=1,
83+
thin=1,
84+
warmup=1000, # Stands for burn-in; Default = iter/2
85+
# seed=123 # Setting seed; Default is random seed
86+
)
87+
# Now the values for the monitored parameters are in the "samples" object,
88+
# ready for inspection.
89+
90+
c1 <- extract(samples_1)$c
91+
r1 <- extract(samples_1)$r
92+
u1 <- extract(samples_1)$u
93+
c2 <- extract(samples_2)$c
94+
r2 <- extract(samples_2)$r
95+
u2 <- extract(samples_2)$u
96+
c3 <- extract(samples_3)$c
97+
r3 <- extract(samples_3)$r
98+
u3 <- extract(samples_3)$u
99+
100+
windows(14, 7)
101+
layout(matrix(1:3, 1, 3, byrow=T))
102+
plot(density(c3), xlim=c(0, 1), lty="dotted")
103+
lines(density(c2), lty="dashed")
104+
lines(density(c1))
105+
plot(density(r3), xlim=c(0, 1), lty="dotted")
106+
lines(density(r2), lty="dashed")
107+
lines(density(r1))
108+
plot(density(u3), xlim=c(0, 1), lty="dotted")
109+
lines(density(u2), lty="dashed")
110+
lines(density(u1))
111+

‎Bayesian_Cognitive_Modeling/CaseStudies/MPT/MPT_2_Stan.R‎

Lines changed: 222 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,222 @@
1+
# clears workspace:
2+
rm(list=ls())
3+
4+
library(rstan)
5+
6+
model <- "
7+
// Multinomial Processing Tree with Latent Traits
8+
data {
9+
int<lower=1> nsubjs;
10+
int<lower=1> nparams;
11+
int<lower=0,upper=20> k[nsubjs,4];
12+
cov_matrix[nparams] I;
13+
}
14+
transformed data {
15+
int df;
16+
vector[3] mudeltahat;
17+
18+
df <- nparams + 1;
19+
mudeltahat[1] <- 0;
20+
mudeltahat[2] <- 0;
21+
mudeltahat[3] <- 0;
22+
}
23+
parameters {
24+
vector[nparams] deltahat[nsubjs];
25+
cov_matrix[nparams] Sigma;
26+
real muchat;
27+
real murhat;
28+
real muuhat;
29+
real<lower=0,upper=100> xichat;
30+
real<lower=0,upper=100> xirhat;
31+
real<lower=0,upper=100> xiuhat;
32+
}
33+
transformed parameters {
34+
simplex[4] theta[nsubjs];
35+
vector<lower=0,upper=1>[nsubjs] c;
36+
vector<lower=0,upper=1>[nsubjs] r;
37+
vector<lower=0,upper=1>[nsubjs] u;
38+
matrix[nparams,nparams] rho;
39+
40+
for (i in 1:nsubjs) {
41+
vector[nsubjs] deltachat;
42+
vector[nsubjs] deltarhat;
43+
vector[nsubjs] deltauhat;
44+
45+
deltachat[i] <- deltahat[i,1];
46+
deltarhat[i] <- deltahat[i,2];
47+
deltauhat[i] <- deltahat[i,3];
48+
49+
// Probitize Parameters c, r, and u
50+
c[i] <- Phi(muchat + xichat * deltachat[i]);
51+
r[i] <- Phi(murhat + xirhat * deltarhat[i]);
52+
u[i] <- Phi(muuhat + xiuhat * deltauhat[i]);
53+
54+
// MPT Category Probabilities for Word Pairs
55+
theta[i,1] <- c[i] * r[i];
56+
theta[i,2] <- (1 - c[i]) * (u[i]) ^ 2;
57+
theta[i,3] <- (1 - c[i]) * 2 * u[i] * (1 - u[i]);
58+
theta[i,4] <- c[i] * (1 - r[i]) + (1 - c[i]) * (1 - u[i]) ^ 2;
59+
}
60+
for (i1 in 1:nparams)
61+
for (i2 in 1:nparams)
62+
rho[i1,i2] <- Sigma[i1,i2] / sqrt(Sigma[i1,i1] * Sigma[i2,i2]);
63+
}
64+
model {
65+
// Priors
66+
muchat ~ normal(0, 1);
67+
murhat ~ normal(0, 1);
68+
muuhat ~ normal(0, 1);
69+
70+
Sigma ~ wishart(df, I);
71+
// Individual Effects
72+
deltahat ~ multi_normal(mudeltahat, Sigma);
73+
// Data
74+
for (i in 1:nsubjs)
75+
k[i] ~ multinomial(theta[i]);
76+
}
77+
generated quantities {
78+
real<lower=0,upper=1> muc;
79+
real<lower=0,upper=1> mur;
80+
real<lower=0,upper=1> muu;
81+
real sigmac;
82+
real sigmar;
83+
real sigmau;
84+
85+
// Post-Processing Means, Standard Deviations, Correlations
86+
muc <- Phi(muchat);
87+
mur <- Phi(murhat);
88+
muu <- Phi(muuhat);
89+
sigmac <- xichat * sqrt(Sigma[1,1]);
90+
sigmar <- xirhat * sqrt(Sigma[2,2]);
91+
sigmau <- xiuhat * sqrt(Sigma[3,3]);
92+
}"
93+
94+
trial_1 <- list(subjs=21, items=20, response=structure(
95+
.Data=c(2,4,4,10,2,1,3,14,2,2,5,11,6,0,4,10,1,0,4,15,1,0,2,17,1,2,4,13,4,1,
96+
6,9,5,1,4,10,1,0,9,10,5,0,3,12,0,1,6,13,1,5,7,7,1,1,4,14,2,2,3,13,2,
97+
1,5,12,2,0,6,12,1,0,5,14,2,1,8,9,3,0,2,15,1,2,3,14),
98+
.Dim = c(4, 21)))
99+
100+
trial_2 <- list(subjs=21, items=20, response=structure(
101+
.Data = c(7,5,3,5,5,2,3,10,6,2,7,5,9,4,2,5,2,2,7,9,1,3,3,13,5,0,5,10,7,3,4,
102+
6,7,3,6,4,4,1,10,5,9,1,2,8,3,1,6,10,3,5,9,3,2,0,6,12,8,0,3,9,3,2,
103+
7,8,7,1,5,7,2,1,6,11,5,3,5,7,5,0,6,9,6,2,2,10),
104+
.Dim = c(4, 21)))
105+
106+
trial_6 <- list(subjs=21, items=20, response=structure(
107+
.Data = c(14,3,1,2,12,3,1,4,18,0,1,1,15,3,0,2,7,1,10,2,3,6,11,0,8,4,3,5,17,
108+
1,1,1,13,4,3,0,11,6,1,2,16,1,2,1,10,1,3,6,7,13,0,0,8,4,3,5,16,1,
109+
1,2,5,4,7,4,15,0,5,0,6,3,6,5,17,2,0,1,17,1,0,2,8,3,6,3),
110+
.Dim = c(4, 21)))
111+
112+
response_1 <- t(trial_1$response)
113+
response_2 <- t(trial_2$response)
114+
response_6 <- t(trial_6$response)
115+
116+
nparams <- 3 # Number of free parameters per participant: c_i, r_i, u_i
117+
I <- diag(3) # Identity matrix for Wishart
118+
119+
k <- response_1
120+
nsubjs <- nrow(k) # Number of word pairs per participant
121+
122+
# To be passed on to Stan
123+
data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I)
124+
125+
myinits <- list(
126+
list(deltahat=matrix(0, 21, 3), Sigma=diag(3),
127+
muchat=0, murhat=0, muuhat=0, xichat=1, xirhat=1, xiuhat=1))
128+
129+
# Parameters to be monitored
130+
parameters <- c("muc", "mur", "muu", "sigmac", "sigmar", "sigmau", "rho")
131+
132+
# Run higher iterations for better estimate
133+
myiterations <- 3300
134+
mywarmup <- 300
135+
136+
# The following command calls Stan with specific options.
137+
# For a detailed description type "?stan".
138+
samples_1 <- stan(model_code=model,
139+
data=data,
140+
init=myinits, # If not specified, gives random inits
141+
pars=parameters,
142+
iter=myiterations,
143+
chains=1,
144+
thin=1,
145+
warmup=mywarmup, # Stands for burn-in; Default = iter/2
146+
# seed=123 # Setting seed; Default is random seed
147+
)
148+
149+
k <- response_2
150+
nsubjs <- nrow(k) # Number of word pairs per participant
151+
# To be passed on to Stan
152+
data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I)
153+
samples_2 <- stan(fit=samples_1,
154+
data=data,
155+
init=myinits, # If not specified, gives random inits
156+
pars=parameters,
157+
iter=myiterations,
158+
chains=1,
159+
thin=1,
160+
warmup=mywarmup, # Stands for burn-in; Default = iter/2
161+
# seed=123 # Setting seed; Default is random seed
162+
)
163+
164+
k <- response_6
165+
nsubjs <- nrow(k) # Number of word pairs per participant
166+
# To be passed on to Stan
167+
data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I)
168+
samples_6 <- stan(fit=samples_1,
169+
data=data,
170+
init=myinits, # If not specified, gives random inits
171+
pars=parameters,
172+
iter=myiterations,
173+
chains=1,
174+
thin=1,
175+
warmup=mywarmup, # Stands for burn-in; Default = iter/2
176+
# seed=123 # Setting seed; Default is random seed
177+
)
178+
# Now the values for the monitored parameters are in the "samples" object,
179+
# ready for inspection.
180+
181+
muc1 <- extract(samples_1)$muc
182+
mur1 <- extract(samples_1)$mur
183+
muu1 <- extract(samples_1)$muu
184+
muc2 <- extract(samples_2)$muc
185+
mur2 <- extract(samples_2)$mur
186+
muu2 <- extract(samples_2)$muu
187+
muc6 <- extract(samples_6)$muc
188+
mur6 <- extract(samples_6)$mur
189+
muu6 <- extract(samples_6)$muu
190+
191+
#### Plots posteriors of the group--level c, r, and u parameters
192+
windows(14,7)
193+
layout(matrix(1:3,1,3,byrow=T))
194+
par(mar = c(4.25, 3, 1, 1))
195+
par(cex.axis=0.9)
196+
plot(density((muc1)),ylim = c(0,6),xlim=c(0,1), axes=F,xlab="c",ylab="",
197+
main="")
198+
axis(1,at=seq(0,1,0.2))
199+
axis(2,at = c(0,12),labels=F,lwd.ticks=0)
200+
mtext("Density",2,1,cex=0.55,las=0)
201+
par(cex=0.60)
202+
legend(0.3, 10, c("Trial1","Trial 2","Trial 3"),lty = c(1,2,3),col=c("black"),
203+
text.col = "black")
204+
par(cex=0.65)
205+
lines(density((muc2)),lty=2)
206+
lines(density((muc6)),lty=3)
207+
208+
plot(density((mur1),to=0.75),ylim = c(0,13),xlim=c(0,1), axes=F,xlab="r",
209+
ylab="",main="")
210+
axis(1,at=seq(0,1,0.2))
211+
axis(2,at = c(0,13),labels=F,lwd.ticks=0)
212+
mtext("Density",2,1,cex=0.55,las=0)
213+
lines(density((mur2)),lty=2)
214+
lines(density((mur6)),lty=3)
215+
216+
plot(density((muu1)),ylim = c(0,9),xlim=c(0,1), axes=F,xlab="u",ylab="",
217+
main="")
218+
axis(1,at=seq(0,1,0.2))
219+
axis(2,at = c(0,10),labels=F,lwd.ticks=0)
220+
mtext("Density",2,1,cex=0.55,las=0)
221+
lines(density((muu2)),lty=2)
222+
lines(density((muu6)),lty=3)

0 commit comments

Comments
 (0)