*! version 1.1 8december2005
*! Jean-Benoit Hardouin
************************************************************************************************************
* Simirt: Simulation of dichotomous data following an IRT model (Rasch model,
* OPLM, Birnbaum model, 3-PLM, 4-PLM, 5-PAM
*
* Version 1 : May 9, 2005 (Jean-Benoit Hardouin)
* Version 1.1 : December 8, 2005 (Jean-Benoit Hardouin) /*group and deltagroup options*/
*
* Jean-benoit Hardouin, Regional Health Observatory of Orléans - France
* jean-benoit.hardouin@orscentre.org
*
* News about this program : http://anaqol.free.fr
* FreeIRT Project : http://freeirt.free.fr
*
* Copyright 2005 Jean-Benoit Hardouin
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
************************************************************************************************************/
program define simirt , rclass
version 8.0
syntax  [, NBObs(integer 2000) Dim(string) MU(string) COV(string) DISc(string) DIFf(string) PMIN(string) PMAX(string) ACC(string) clear STOre(string) REPlace PREFix(string) DRAW GRoup(real 0) DELtagroup(real 0)]


/********************************************************************************
TESTS
********************************************************************************/


if `group'<0|`group'>1 {
   di in red "{p}The {hi:group} option defines a probability. The values defined by this option must be greater (or equal) to 0 and lesser (or equal) to 1.{p_end}"
   error 198
   exit
}

if "`clear'"==""&"`store'"=="" {
   di in red "You must use at least one of these two options: clear and/or store."
   error 198
   exit
}

if "`dim'"!="" {
   local nbdim:word count `dim'
   if `nbdim'>2 {
      di in red "You can simulate data with one or two dimensions, and you have indicated `nbdim' dimensions in the {hi:dim} option. Please correct it."
      error 198
      exit
   }
   if `nbdim'==1 {
      local dim1=`dim'
      local nbitems=`dim1'
      local dim=1
   }
   else if `nbdim'==2 {
      local dim1:word 1 of `dim'
      local dim2:word 2 of `dim'
      local nbitems=`dim1'+`dim2'
      local dim=2
   }
   else {
      di in red "Your {hi:dim} option is uncorrect, please correct it. This option must indicate one or two integer(s)."
      error 198
      exit
   }
   if "`diff'"!="" {
      local nbdiff:word count `diff'
      if (`nbdiff'==3&`dim'==1)|(`nbdiff'==5&`dim'==2)|`nbdiff'==1 {
         local typediff:word 1 of `diff'
         if "`typediff'"!="gauss"&"`typediff'"!="uniform" {
            local typediff values
         }
      }
      if "`typediff'"=="values"&`nbdiff'!=`nbitems' {
          di in red "You have indicated a number of difficulty parameters ({hi:diff} option) different of the number of items to simulate ({hi:dim} option). Please correct these options."
          error 198
          exit
      }
   }
   if "`diff'"=="" {
      if `dim'==1 {
         local diff gauss 0 1
      }
      else if `dim'==2 {
         local diff gauss 0 1 0 1
      }
      local typediff gauss
      local nbdiff:word count `diff'
   }
}
else if "`dim'"==""{
   if "`diff'"=="" {
      di in red "You must indicate the number of items to simulate with the {hi:dim} or the {hi:diff} option(s)."
      error 198
      exit
   }
   else {
      local nbdiff:word count `diff'
      if `nbdiff'!=`nbitems' {
          di in red "You have indicated a number of difficulty parameters ({hi:diff} option) different of the number of items to simulate ({hi:dim} option). Please correct these options."
          error 198
          exit
      }
   }
}
if "`draw'"!=""&`dim'!=1 {
   di in red "The {hi:draw} option is available only with unidimensional simulated data."
   error 198
   exit
}
if (`group'!=0|`deltagroup'!=0)&`dim'!=1 {
   di in red "The {hi:group} and the {hi:deltagroup} options are available only with unidimensional simulated data."
   error 198
   exit
}

if "`prefix'"=="" {
   local prefix item
}

local nbprefix:word count `prefix'
local prefix1:word 1 of `prefix'
if `dim'==2&`nbprefix'>=2 {
   local prefix2:word 2 of `prefix'
}
else if `dim'==2 {
   local prefix2:word 1 of `prefix'
}
if `dim'==2&"`prefix1'"=="`prefix2'" {
   local prefix1 `prefix1'A
   local prefix2 `prefix2'B
}

*di in ye "dim : `dim' ; diff : `diff' ; nbdiff : `nbdiff'"

local nbcov:word count `cov'

if `dim'==1&`nbcov'>1 {
   di in red "You simulate one dimension. You must indicate only the variance of the simulated latent trait in the {hi:cov} option."
   error 198
   exit
}
else if `dim'==2&`nbcov'!=3&`nbcov'>0 {
   di in red "You simulate two dimensions. You must indicate exactly 3 values in the {hi:cov} option (Variance of the first simulated latent trait, Variance of the second simulated latent trait, Covariance between the two simulated latent traits)."
   error 198
   exit
}
else if `nbcov'==0 {
   if `dim'==1 {
      local cov "1"
   }
   else if `dim'==2 {
      local cov "1 1 0"
   }
   local nbcov:word count `cov'
}

if `nbcov'==1 {
   if `cov'<0 {
      di in red "The variance of your latent trait can not be negative. Please correct your {hi:cov} option."
      error 198
   }
}
else if `nbcov'==3 {
   local cov1:word 1 of `cov'
   local cov2:word 2 of `cov'
   local cov3:word 3 of `cov'
   local rho=`cov3'/sqrt(`cov1'*`cov2')
   if `cov1'<0|`cov2'<0|`rho'<-1|`rho'>1 {
      di in red "Your covariance matrix defined by the {hi:cov} option is not correct. Please correct it."
      error 198
      exit
   }
}

local nbmu:word count `mu'
if `nbmu'!=`dim'&`nbmu'!=0 {
   di in red "You must indicate as many values in the {hi:mu} option as the number of dimension(s) (`dim')"
   error 198
   exit
}
local nbdisc:word count `disc'
if `nbdisc'!=`nbitems'&`nbdisc'!=0 {
   di in red "You must indicate as many values in the {hi:disc} option as items defined by the {hi:dim} and the {hi:diff} options (`nbitems')"
   error 198
   exit
}
local nbpmin:word count `pmin'
if `nbpmin'!=`nbitems'&`nbpmin'!=0 {
   di in red "You must indicate as many values in the {hi:pmin} option as items defined by the {hi:dim} and the {hi:diff} options (`nbitems')"
   error 198
   exit
}
local nbpmax:word count `pmax'
if `nbpmax'!=`nbitems'&`nbpmax'!=0 {
   di in red "You must indicate as many values in the {hi:pmax} option as items defined by the {hi:dim} and the {hi:diff} options (`nbitems')"
   error 198
   exit
}
local nbacc:word count `acc'
if `nbacc'!=`nbitems'&`nbacc'!=0 {
   di in red "You must indicate as many values in the {hi:acc} option as items defined by the {hi:dim} and the {hi:diff} options (`nbitems')"
   error 198
   exit
}

preserve

tempfile saveraschbin
capture qui save `saveraschbin'




/********************************************************************************
PARAMETERS
********************************************************************************/


/*
scalar define hour=real(substr("$S_TIME",1,2))
scalar define min=real(substr("$S_TIME",4,2))
scalar define sec=real(substr("$S_TIME",7,2))
scalar define jour=real(substr("$S_DATE",1,2))
*/
local hour=real(substr("$S_TIME",1,2))
local min=real(substr("$S_TIME",4,2))
local sec=real(substr("$S_TIME",7,2))
local jour=real(substr("$S_DATE",1,2))




if "$seed"!="" {
   global seed2=int($seed)
}
else {
   global seed2=0
}
global seed=$seed2+256484+`sec'*1000000+`min'*10000+`hour'*100+`jour'
if $seed>2^31-1 {
   global seed=int($seed/10)
}
qui set seed $seed

if "`typediff'"=="uniform" {
   if `dim'==1 {
      if `nbdiff'==3 {
         local min:word 2 of `diff'
         local max:word 3 of `diff'
      }
      else if `nbdiff'==1 {
         local min=-2
         local max=2
      }
      else {
         di in red "Your {hi:diff} option is uncorrect. Please correct it."
      }
      local diff
      forvalues i=1/`nbitems' {
          local diff `diff' `=`min'+(`max'-`min')*`i'/(`nbitems'+1)'
      }
   }
   if `dim'==2 {
      if `nbdiff'==5 {
         local min1:word 2 of `diff'
         local max1:word 3 of `diff'
         local min2:word 4 of `diff'
         local max2:word 5 of `diff'
      }
      else if `nbdiff'==1 {
         local min1=-2
         local max1=2
         local min2=-2
         local max2=2
      }
      else {
         di in red "Your {hi:diff} option is uncorrect. Please correct it."
         exit
      }
      local diff
      forvalues i=1/`dim1' {
          local diff `diff' `=`min1'+(`max1'-`min1')*`i'/(`dim1'+1)'
      }
      forvalues i=1/`dim2' {
          local diff `diff' `=`min2'+(`max2'-`min2')*`i'/(`dim2'+1)'
      }
   }
}
if "`typediff'"=="gauss" {
   if `dim'==1 {
      if `nbdiff'==3 {
         local mean:word 2 of `diff'
         local var:word 3 of `diff'
      }
      else if `nbdiff'==1 {
         local mean=0
         local var=1
      }
      else {
         di in red "Your {hi:diff} option is uncorrect. Please correct it."
         error 198
         exit
      }
      local diff
      forvalues i=1/`nbitems' {
          local tmp=invnorm(`i'/(`nbitems'+1))*sqrt(`var')+`mean'
          local diff `diff' `tmp'
      }
   }
   if `dim'==2 {
      if `nbdiff'==5 {
         local mean1:word 2 of `diff'
         local var1:word 3 of `diff'
         local mean2:word 4 of `diff'
         local var2:word 5 of `diff'
      }
      else if `nbdiff'==1 {
         local mean1=0
         local var1=1
         local mean2=0
         local var2=1
      }
      else {
        di in red "Your {hi:diff} option is uncorrect. Please correct it."
        error 198
        exit
      }
      local diff
      forvalues i=1/`dim1' {
          local tmp=invnorm(`i'/(`dim1'+1))*sqrt(`var1')+`mean1'
          local diff `diff' `tmp'
      }
      forvalues i=1/`dim2' {
          local tmp=invnorm(`i'/(`dim2'+1))*sqrt(`var2')+`mean2'
          local diff `diff' `tmp'
      }
   }
}


tempname matmu matcov matdiff matdisc matpmin matpmax matacc
matrix define `matdiff'=J(`nbitems',1,0)
matrix define `matdisc'=J(`nbitems',1,0)
matrix define `matpmin'=J(`nbitems',1,0)
matrix define `matpmax'=J(`nbitems',1,0)
matrix define `matacc'=J(`nbitems',1,0)
matrix define `matmu'=J(`dim',1,0)
matrix define `matcov'=J(`=(`dim'+1)*`dim'/2',1,0)

forvalues i=1/`nbitems' {
   local tmp:word `i' of `diff'
   matrix `matdiff'[`i',1]=`tmp'
   if `nbdisc'!=0 {
      local tmp:word `i' of `disc'
      matrix `matdisc'[`i',1]=`tmp'
   }
   else {
      matrix `matdisc'[`i',1]=1
   }
   if `nbpmin'!=0 {
      local tmp:word `i' of `pmin'
      matrix `matpmin'[`i',1]=`tmp'
   }
   else {
      matrix `matpmin'[`i',1]=0
   }
   if `nbpmax'!=0 {
      local tmp:word `i' of `pmax'
      matrix `matpmax'[`i',1]=`tmp'
   }
   else {
      matrix `matpmax'[`i',1]=1
   }
   if `nbacc'!=0 {
      local tmp:word `i' of `acc'
      matrix `matacc'[`i',1]=`tmp'
   }
   else {
      matrix `matacc'[`i',1]=1
   }
}

if `nbcov'==3|`nbcov'==1 {
   local tmp:word 1 of `cov'
   matrix `matcov'[1,1]=`tmp'
   if `nbcov'==3&`dim'==2 {
      local tmp:word 2 of `cov'
      matrix `matcov'[2,1]=`tmp'
      local tmp:word 3 of `cov'
      matrix `matcov'[3,1]=`tmp'
   }
}
else {
   matrix `matcov'[1,1]=1
   if `dim'==2 {
      matrix `matcov'[2,1]=1
      matrix `matcov'[3,1]=0
   }
}
if (`nbmu'==`dim') {
   local tmp:word 1 of `mu'
   matrix `matmu'[1,1]=`tmp'
   if `dim'==2 {
      local tmp:word 2 of `mu'
      matrix `matmu'[2,1]=`tmp'
   }
}
if `dim'==2 {
   local corr=`matcov'[3,1]/sqrt(`matcov'[1,1]*`matcov'[2,1])
}

/********************************************************************************
ITEMS CHARACTERISTIC CURVES
********************************************************************************/


if "`draw'"!="" {
   *set trace on
   drop _all
   qui set obs 2001
   qui gen lt1=_n
   qui replace lt1=(lt1-1001)/1000*4*sqrt(`cov')+`matmu'[1,1]
   label variable lt1 "Latent trait"
   local dess

   forvalues i=1/`dim1' {
      qui gen `prefix1'`i'=`matpmin'[`i',1]+(`matpmax'[`i',1]-`matpmin'[`i',1])*(1/(1+exp(-`matdisc'[`i',1]*(lt1-`matdiff'[`i',1]))))^(`matacc'[`i',1])
      local dess `dess' (line `prefix'`i' lt1)
   }
   graph twoway `dess' , ylabel(0(.25)1) legend(off) ytitle("Probability of a positive response")
}


/********************************************************************************
SIMULATION
********************************************************************************/

drop _all
qui set obs `nbobs'
qui gen lt1=invnorm(uniform())
if `dim'==2 {
   qui gen lt2=invnorm(uniform())
   qui replace lt2=`corr'*lt1+sqrt(1-(`corr')^2)*lt2
   qui replace lt2=lt2*sqrt(`matcov'[2,1])+`matmu'[2,1]
}
qui replace lt1=lt1*sqrt(`matcov'[1,1])+`matmu'[1,1]

if `dim'==1&`group'!=0 {
   qui gen group=uniform()<`group'
   qui replace lt1=lt1+`deltagroup'*(1-`group') if group==1
   qui replace lt1=lt1-`deltagroup'*`group' if group==0
}

di in gr "{hline 75}"
di _col(1) in gr "Items" _col(18) "Difficulty" _col(34) "Discr." _col(45) "Pmin" _col(58) "Pmax" _col(73) "Acc"
di in gr "{hline 75}"
forvalues i=1/`dim1' {
   qui gen `prefix1'`i'=`matpmin'[`i',1]+(`matpmax'[`i',1]-`matpmin'[`i',1])*(1/(1+exp(-`matdisc'[`i',1]*(lt1-`matdiff'[`i',1]))))^(`matacc'[`i',1])
   qui replace `prefix1'`i'=uniform()<=`prefix1'`i'
   di _col(1) in gr "`prefix1'`i'" _col(20) in ye %8.4f `matdiff'[`i',1] _col(32) %8.4f `matdisc'[`i',1] _col(44) %6.4f `matpmin'[`i',1] _col(56) %6.4f `matpmax'[`i',1] _col(68) %8.4f `matacc'[`i',1]
}
if `dim'==2 {
   di in gr "{hline 75}"
   forvalues i=`=`dim1'+1'/`=`dim1'+`dim2'' {
      qui gen `prefix2'`=`i'-`dim1''=`matpmin'[`i',1]+(`matpmax'[`i',1]-`matpmin'[`i',1])*(1/(1+exp(-`matdisc'[`i',1]*(lt2-`matdiff'[`i',1]))))^(`matacc'[`i',1])
      qui replace `prefix2'`=`i'-`dim1''=uniform()<=`prefix2'`=`i'-`dim1''
      di _col(1) in gr "`prefix2'`=`i'-`dim1''" _col(20) in ye %8.4f `matdiff'[`i',1] _col(32) %8.4f `matdisc'[`i',1] _col(44) %6.4f `matpmin'[`i',1] _col(56) %6.4f `matpmax'[`i',1] _col(68) %8.4f `matacc'[`i',1]
   }
}
di in gr "{hline 75}"
di



qui corr lt*,cov
local var_1=r(Var_1)
return scalar var_1=r(Var_1)
if `dim'==2 {
   local cov_12=r(cov_12)
   local var_2=r(Var_2)
   return scalar cov_12=r(cov_12)
   return scalar var_2=r(Var_2)
   qui corr lt*
   local rho=r(rho)
   return scalar rho=r(rho)
}
qui su lt1
local mean1=r(mean)
return scalar mean_1=r(mean)
if `dim'==2 {
   qui su lt2
   local mean2=r(mean)
   return scalar mean_2=r(mean)
}
if `dim'==1&`group'!=0 {
   qui su lt1 if group==0
   local mean_0=r(mean)
   qui su lt1 if group==1
   local mean_1=r(mean)
   local delta=`mean_1'-`mean_0'
}

return scalar nbobs=`nbobs'

di in gr "{hline 50}"
di _col(1) in gr  "Latent trait" _c
if `dim'==2 {
   di in gr "s"_c
}
di in gr _col(29) "Theorical" _col(42) "Observed"
di in gr "{hline 50}"
forvalues i=1/`dim' {
   di _col(1) in gr "Mean(lt`i')" _col(30) in ye %8.4f `matmu'[`i',1] _col(42) %8.4f `mean`i''
   di _col(1) in gr "Variance(lt`i')" _col(30) in ye %8.4f `matcov'[`i',1] _col(42) %8.4f `var_`i''
}
if `dim'==1&`group'!=0 {
   di _col(1) in gr "Mean(lt1) group 0" _col(30) in ye %8.4f `matmu'[1,1]-`deltagroup'*`group' _col(42) %8.4f `mean_0'
   di _col(1) in gr "Mean(lt1) group 1" _col(30) in ye %8.4f `matmu'[1,1]+`deltagroup'*(1-`group') _col(42) %8.4f `mean_1'
   qui count if group==1
   local prop=r(N)/`nbobs'
   di _col(1) in gr "Proportion group 1" _col(30) in ye %8.4f `group' _col(42) %8.4f `prop'
}
if `dim'==2 {
   di _col(1) in gr "Covariance" _col(30) in ye %8.4f `matcov'[3,1] _col(42) %8.4f `cov_12'
   di _col(1) in gr "Correlation" _col(31) in ye %7.4f `corr' _col(43) %7.4f `rho'
}
di in gr "{hline 50}"

/********************************************************************************
CLEAR AND/OR STORE
********************************************************************************/


if "`clear'"!="" {
   restore, not
   preserve
}
if "`store'"!="" {
   save "`store'",`replace'
}
if "`clear'"=="" {
   use "`saveraschbin'",replace
}
end