#----------------------------------
#Section 1 - Flexible parameters
#----------------------------------

library(ltm)

NSim <- 10 #Number of simulations run

#Number of items in a test
N<-10

#Number of categories in each item (it is best it k be an odd number)
k<-5

#Step. The change in difficulty between each of the item's possible categories.
#recall that we have k categories whithin each item. The difficulties of each categrory differ in h
h <- 0.4

#Number of respondents that take the test
S <- 150

#In our simulations, we consider that subjects have two latent traits: one which the test was designed to measure and another which, albeit distinct, is correlated to the first.
#The squared correlation between both latent trait is the following r2:
r2<-0.36    

#Ogive D parameter (1 = logistic, 1.7 = normal)
D = 1

#This function generates a response given by a subject  to an item
resp<-function(subject,item,trait)
{
  #Calculate the probability of each category
  
  p.vec = NULL
  
  for(category in 1:k) #For each category, calculate:
  {
    p = ogive(theta[subject,trait],a[item],b[item,category],D)
    p.vec = c(p.vec,p)
  }
  
  #Generates a random number
  
  seed<-runif(1,0.2,0.8)
  
  #Selects an answer
  
  temp = which(p.vec>seed) #we use this temporary vector in order to avoid warnings should its length be zero
  if(length(temp)>0){resp=max(temp)}else{resp=1}
  resp
}

for(count in 1:NSim)
{
  '*************************************************************************************************************'
  
  'Items average bs are uniformly distributed between -2 and 1.
For now, we ignore other possible distributions for b.
  If everything works, then we may venture into considering other distributions for b'

  bmed <- runif(N,-2.5,2.5)
  
  #Sorts items acording to their difficulty
  bmed<-sort(bmed)
  
  #Creates a discrimination parameter for each of the N items in the test.
  a <- rnorm(N,1.25,0.07)
  
  #Creates a main latent trait which the subjects posess and which the test is designed to measure.
  theta1<-rnorm(S,0,1)
  #Creates the secondary trait whose squared correlation with theta1 was defined by the parameter r2
  sd.epsilon<-sd(theta1)*sqrt((1/r2)-1)   #Defines the standard deviation of the error term to be added to theta in order to define theta2
  theta2<-theta1+rnorm(S,0,sd.epsilon)
  
  
  #=======================================================
  #The user need enter information only up to this point.
  #Here we generate the test
  #=======================================================
  
  #Creates a matrix of bs.
  #Each row of this matrix designates one of the N items in the test.
  #Each column designates one of the k categories in its respective item (line)
  b<- matrix(rep(NA,  N * k), ncol = k)
  centro<-(k+1)/2
  for (i in 1:N)
  {
    for (j in 1:k)
    {
      b[i,j]<-bmed[i]+(j-(k+1)/2)*h #constant stem with bmed being the average b. The (k+1)/2 works if k is an odd number 
    }
  }
  #Here, we had a piece of code that sorted the bs whithin each line.
  #However, as we have constructed the b's now (by adding of subtracting multiples of h), there is no need for sorting.
  #If, however, we wish to consider randomly distributed values of h in the future, then we will have to sort these values before adding them to bmed
  
  
  ogive<-function(x,a,b,D)
  {
    1/(1+exp(-D*a*(x-b)))
  }
  
  #----------------------------------------------------------------------------------------------------------
  #I have made this piece of code from scratch substituting a previous translation of an Excel simulation.
  #Arguably, this portion of the code may need some correction.
  #----------------------------------------------------------------------------------------------------------
  
  
  #theta is a matrix where each column is a theta, so we can run everything for one column and then for the other one
  theta <- as.matrix(cbind(theta1,theta2))
  
  SOR <- NULL
  SOR <- array(rep(NA,S*N*2),dim=c(S,N,2))
  #The Stock_of_Responses (SOR) is an array composed of 2 matrices, one referring to latent trait theta1 and the other to latent trait theta2.
  #Each of this matrices is an SxN matrix, with one row for each respondent and one column for each of the test N items.
  #Thus, each cell designates the category that the row's subject selected when responding to the column's item using the latent trait of the respective matrix.
  #i.e. SOR[3,2,1] is the category selected by the 3rd subject to the 2nd question using latent trait theta1.
    
  for(trait in 1:2)
  {
    #Here we generate the responses for each individual to each item and stock them in the SOR
    for(i in 1:S){for(j in 1:N) {SOR[i,j,trait]=resp(i,j,trait)}}
  }
  
  
  #Choose the items to be included in each category
  E<-c(1,2,3,7)                    #Easy items
  I<-c(4,5,6,7)                    #Intermediate items
  D<-c(3,8,9,10)                   #Difficult items
  Sp<-c(2,4,6,8)                    #Sparced items
  All<-c(1,2,3,4,5,6,7,8,9,10)     #All items
  
  
  #The followinf function creates a matrix to input the L responses of interest for each of the S subjects
  IM<-function(items,responses,trait) 
  {
    
    L<-length(items)
    A<-matrix(rep(NA,S*L),nrow=S)  
    
    for(i in 1:S) #for each person...
    {  
      for(j in 1:L)
      {
        A[i,j]<-responses[i,items[j],trait]  #...gather the responses s/he gave to each of the items of interest.
      }
    }
    
    A
  }
  
  IME1<-IM(E,SOR,1)      #Input matrix for easy items on latent trait 1. This contains all relevant data from this items.   
  IME2<-IM(E,SOR,2)
  IMI1<-IM(I,SOR,1)   
  IMI2<-IM(I,SOR,2)
  IMD1<-IM(D,SOR,1)   
  IMD2<-IM(D,SOR,2)   
  IMS1<-IM(Sp,SOR,1)   
  IMS2<-IM(Sp,SOR,2)
  IMA1<-IM(All,SOR,1)   
  IMA2<-IM(All,SOR,2)
  
  IRTE1 <- grm(IME1, constrained=FALSE,start.val="random")
  IRTE2 <- grm(IME2, constrained=FALSE,start.val="random")
  IRTI1 <- grm(IMI1, constrained=FALSE,start.val="random")
  IRTI2 <- grm(IMI2, constrained=FALSE,start.val="random")
  IRTD1 <- grm(IMD1, constrained=FALSE,start.val="random")
  IRTD2 <- grm(IMD2, constrained=FALSE,start.val="random")
  IRTS1 <- grm(IMS1, constrained=FALSE,start.val="random")
  IRTS2 <- grm(IMS2, constrained=FALSE,start.val="random")
  IRTA1 <- grm(IMA1, constrained=FALSE,start.val="random")
  IRTA2 <- grm(IMA2, constrained=FALSE,start.val="random")
}