if (FALSE) {
## Euro 2016 example
data(Euro2016)
posterior1 <- MCMCpaircompare(pwc.data=Euro2016,
theta.constraints=list(Ukraine="-",
Portugal="+"),
alpha.fixed=TRUE,
verbose=10000,
burnin=10000, mcmc=500000, thin=100,
store.theta=TRUE, store.alpha=FALSE)
## alternative identification constraints
posterior2 <- MCMCpaircompare(pwc.data=Euro2016,
theta.constraints=list(Ukraine="-",
Portugal=1),
alpha.fixed=TRUE,
verbose=10000,
burnin=10000, mcmc=500000, thin=100,
store.theta=TRUE, store.alpha=FALSE)
## a synthetic data example with estimated rater-specific parameters
set.seed(123)
I <- 65 ## number of raters
J <- 50 ## number of items to be compared
## raters 1 to 5 have less sensitivity to stimuli than raters 6 through I
alpha.true <- c(rnorm(5, m=0.2, s=0.05), rnorm(I - 5, m=1, s=0.1))
theta.true <- sort(rnorm(J, m=0, s=1))
n.comparisons <- 125 ## number of pairwise comparisons for each rater
## generate synthetic data according to the assumed model
rater.id <- NULL
item.1.id <- NULL
item.2.id <- NULL
choice.id <- NULL
for (i in 1:I){
for (c in 1:n.comparisons){
rater.id <- c(rater.id, i+100)
item.numbers <- sample(1:J, size=2, replace=FALSE)
item.1 <- item.numbers[1]
item.2 <- item.numbers[2]
item.1.id <- c(item.1.id, item.1)
item.2.id <- c(item.2.id, item.2)
eta <- alpha.true[i] * (theta.true[item.1] - theta.true[item.2])
prob.item.1.chosen <- pnorm(eta)
u <- runif(1)
if (u <= prob.item.1.chosen){
choice.id <- c(choice.id, item.1)
}
else{
choice.id <- c(choice.id, item.2)
}
}
}
item.1.id <- paste("item", item.1.id+100, sep=".")
item.2.id <- paste("item", item.2.id+100, sep=".")
choice.id <- paste("item", choice.id+100, sep=".")
sim.data <- data.frame(rater.id, item.1.id, item.2.id, choice.id)
## fit the model
posterior <- MCMCpaircompare(pwc.data=sim.data,
theta.constraints=list(item.101=-2,
item.150=2),
alpha.fixed=FALSE,
verbose=10000,
a=0, A=0.5,
burnin=10000, mcmc=200000, thin=100,
store.theta=TRUE, store.alpha=TRUE)
theta.draws <- posterior[, grep("theta", colnames(posterior))]
alpha.draws <- posterior[, grep("alpha", colnames(posterior))]
theta.post.med <- apply(theta.draws, 2, median)
alpha.post.med <- apply(alpha.draws, 2, median)
theta.post.025 <- apply(theta.draws, 2, quantile, prob=0.025)
theta.post.975 <- apply(theta.draws, 2, quantile, prob=0.975)
alpha.post.025 <- apply(alpha.draws, 2, quantile, prob=0.025)
alpha.post.975 <- apply(alpha.draws, 2, quantile, prob=0.975)
## compare estimates to truth
par(mfrow=c(1,2))
plot(theta.true, theta.post.med, xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5),
col=rgb(0,0,0,0.3))
segments(x0=theta.true, x1=theta.true,
y0=theta.post.025, y1=theta.post.975,
col=rgb(0,0,0,0.3))
abline(0, 1, col=rgb(1,0,0,0.5))
plot(alpha.true, alpha.post.med, xlim=c(0, 1.2), ylim=c(0, 3),
col=rgb(0,0,0,0.3))
segments(x0=alpha.true, x1=alpha.true,
y0=alpha.post.025, y1=alpha.post.975,
col=rgb(0,0,0,0.3))
abline(0, 1, col=rgb(1,0,0,0.5))
}
Run the code above in your browser using DataLab