This is a questionnaire generator based on a Multilevel Partial Credit model. This code is partially based on the Okan Bulut’s work. You can find the Simulating Response Data With Partial Credit Model webpage here

You can generate questionnaire responses by school based on both item and person parameters.

set.seed(1234)

items<-4 ## number of items
cat<-5### number of categories by item
n.col<-1000 ## number of schools

n<-rbinom(n.col,30,0.7) ## random number of students by school.
# In this case, it was simulated very close to 30 students (with probability 0.7)
# you should change the probability parameter in rbinom function if
# you want to change the the number of students

eta.school<-rnorm(n.col) ## eta by school (school ability parameter)
theta<-lapply(1:n.col, function(a) rnorm(n[a],eta.school[a],1)) ### theta=eta+theta_p

### item parameters
beta <- rnorm(items,0,1)
tau <- t(replicate(items, diff(c(0, sort(runif((cat-2),-1,1)), 0))))
for(i in 1:items)
{tau[i,]=sort(tau[i,])}

names.tau<-c(paste("tau",1:(items-1)))
item.par <- cbind(beta, tau)
colnames(item.par)<-c("beta",paste("tau",1:(cat-1)))

q.school<-list() ## list that contains the questionnaires
for(m in 1:n.col)
{
responses <- matrix(NA, nrow = length(theta[[m]]), ncol = items) ## response matrix for each questionnaire

## This function generates the questionnaries
for(k in 1:length(theta[[m]])) {
for(i in 1:items) {
measure<-0
cumexp <- vector()
cumexp <- 1
for(l in 2:cat) {
measure <- measure + theta[[m]][k]-item.par[i,1]-item.par[i,l]
cumexp[l] <- cumexp[(l-1)]+exp(measure)
}

U <- runif(1, 0, 1)
U <- U * cumexp[cat]

for(j in 1:cat) {
if(U <= cumexp[j]) {
responses[k,i] <- (j-1)
break
}
}
}
}
responses <- as.data.frame(responses) # each questionnarie is kept in a data frame
q.school[[m]]<-data.frame(school=rep(m,dim(responses)),responses,theta=theta[[m]]) ## keep the id school, questionnaire and theta
}

q.school.all<-do.call(rbind,q.school)
q.school.all <- as.data.frame(q.school.all)