*! version 2.4 june2020
*! Myriam Blanchin - Priscilla Brisson
************************************************************************************************************
* ROSALI: RespOnse-Shift ALgorithm at Item-level
* Response-shift detection based on Rasch models family
*
* Version 1 : December 21, 2016 (Myriam Blanchin) /*rspcm122016*/
* Version 1.1 : October 13, 2017 (Myriam Blanchin) /*option: MODA, automatic recoding of unused response categories*/
* Version 2 : April, 2018 (Myriam Blanchin - Priscilla Brisson) /*option: GROUP, dichotomous group variable*/
* Version 2.1 : October, 2018 (Myriam Blanchin - Priscilla Brisson) /* Version 1.1 + Version 2 */
* Version 2.2 : February, 2019 (Priscilla Brisson) /* option nodif, optimization */
* Version 2.3 : December, 2019 (Priscilla Brisson) /* option detail, + petites corrections */
* Version 2.4 : June, 2020 (Myriam Blanchin) /* debug option detail + step C, modifs sorties et help */
*
* Myriam Blanchin, SPHERE, Faculty of Pharmaceutical Sciences - University of Nantes - France
* myriam.blanchin@univ-nantes.fr
*
* Priscilla Brisson, SPHERE, Faculty of Pharmaceutical Sciences - University of Nantes - France
* priscilla.brisson@univ-nantes.fr
*
* 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 rosali_original, rclass

timer clear 1
timer on 1

syntax varlist(min=2 numeric) [if] [,GROUP(varlist) NODIF PRO DETail]

preserve
version 15
tempfile saverspcm
capture qui save `saverspcm',replace
local save1=_rc

if "`if'"!="" {
    qui keep `if'
}

if "`pro'" != "" {
	di "START"
}

/**************************************************************************/
set more off
set matsize 5000
constraint drop _all

local gp "`group'"

tokenize `varlist'
local nbitems:word count `varlist'

	/* Vérif nb d'items pair */
local mod=mod(`nbitems',2)
if `mod'!=0 {
  di as error "You must enter an even number of items : the first half of the items represents the items at time 1 and the second half the items at time 2"
  error 198
  exit
}

local nbitems=`nbitems'/2


if "`group'"=="" & "`nodif'"!="" {
	di as error "nodif can only be used with the group option ({hi:nodif} option). Please correct this option."
	error 198
	exit
}

local nbc: word count `group'
if `nbc' >= 2 {
	di as error "Only one variable can be used for group option ({hi:group} option). Please correct this option."
	error 198
	exit
}

	/* Vérif qu'il y a 2 groupes si l'option groupe est choisie */
if "`group'"!="" {
  qui tab `group'
  local nbgrp = r(r)
  if `nbgrp' != 2 {
	di as error "Only 2 groups are possible for the group option ({hi:group} option). Please correct this option."
	error 420
	exit
  }
}
/* recoder la variable de groupe en 0, 1*/

if "`group'"!="" {
	qui tab `gp', matrow(rep)
	qui matrix list rep
	if rep[1,1]+rep[2,1] != 1 & rep[1,1]*rep[2,1] != 0 {
		forvalues i=1/`=rowsof(rep)'{
			qui replace `gp'=`i'-1 if `gp'==rep[`i',1]
			di "WARNING : `gp' `=rep[`i',1]' is now `gp' `=`i'-1' "
		}
	}
	forvalues g = 0/1 {
		qui tab `gp' if `gp' == `g'
		local nbp_gp`g' = r(N)
	}
}



/*item rename*/
/*
Items au temps 1 : 1 à nbitems ``j''
Items au temps 2 : nbitems à 2*nbitems ``=`j'+`nbitems'''

Si t varie, puis num item : ``=(`t'-1)*`nbitems'+`j'''
*/


local com_z = 0 // Indicatrice de recodage
	/*verif modalités répondues*/
if "`gp'" == "" {							// Si pas d'option groupe
	forvalues j = 1 / `nbitems' {
		local recoda_`j' = 0
		qui tab ``j'', matrow(rect1_`j')				// Récupération des infos moda du temps 1
		local minm`j'_t1 = rect1_`j'[1,1]
		local maxm`j'_t1 = rect1_`j'[r(r),1]

		qui tab ``=`j'+`nbitems''', matrow(rect2_`j')			// Récupération des infos moda du temps 2
		local minm`j'_t2 = rect2_`j'[1,1]
		local maxm`j'_t2 = rect2_`j'[r(r),1]

		local minm_`j' = min(`minm`j'_t1',`minm`j'_t2')		// Info moda pour l'item j
		local maxm_`j' = max(`maxm`j'_t1',`maxm`j'_t2')
		local nbm_`j' = `=`maxm_`j''-`minm_`j'''

		if `minm_`j'' != 0 & `com_z' == 0 {
			local com_z = 1
		}


		//Recodage des réponses en 0, 1, 2, etc...
		forvalues r = 0/`=`maxm_`j''-1'  {
				qui replace ``j'' = `r' if ``j'' == `=`r'+`minm_`j'''
				qui replace ``=`j'+`nbitems''' = `r' if ``=`j'+`nbitems''' == `=`r'+`minm_`j'''
		}

		// Vérif. Que toutes les modas sont utilisées & concordance entre temps
		forvalues m = 0/`nbm_`j'' {
			qui count if ``j'' == `m'
			local nb_rn1 = r(N)
			qui count if ``=`j'+`nbitems''' == `m'
			local nb_rn2 = r(N)
			local nb_rn = min(`nb_rn1',`nb_rn2')

			if `nb_rn' == 0 {		// Une moda n'est pas utilisée
				local recoda_`j' = 1
				if   `m' == 0 | `m' <= `minm`j'_t1' | `m' <= `minm`j'_t2' { // La moda 0 ou les moda min ne sont pas utilisées
					local stop = 1
					forvalues k = 1/`=`nbm_`j''-`m'' {
						qui count if ``j'' == `=`m' + `k''
						local v`k'1 = r(N)
						qui count if ``=`j'+`nbitems''' == `=`m' + `k''
						local v`k'2 = r(N)
						if (`v`k'1' != 0 | `v`k'2' != 0) & `stop' != 0 {
							qui replace ``j''= `=`m'+`k'' if ``j''==`m'
							qui replace ``=`j'+`nbitems'''=`=`m'+`k'' if ``=`j'+`nbitems'''==`m'
							di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `=`m'+`k'' merged "
							local stop = 0
						}
					}
				}
				else if `m' >= `maxm`j'_t1' | `m' >= `maxm`j'_t2' | `m' == `maxm_`j'' {	// La (ou les) moda max ne sont pas utilisée(s)
					local stop = 1
					forvalues k = 1/`m' {
						qui count if ``j'' == `=`m' - `k''
						local v`k'1 = r(N)
						qui count if ``=`j'+`nbitems''' == `=`m' - `k''
						local v`k'2 = r(N)
						if (`v`k'1' != 0 | `v`k'2' != 0) & `stop' != 0 {
							qui replace ``j''=`=`m' - `k'' if ``j''==`m'
							qui replace ``=`j'+`nbitems'''=`=`m' - `k'' if ``=`j'+`nbitems'''==`m'
							di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `=`m'-`k'' merged"
							local stop = 0
						}
					}
				}
				else {
					if runiform()>0.5{											// Tirage au sort pour regrouper
						local stop = 1
						forvalues k = 1/`m' {
							qui count if ``j'' == `=`m' - `k''
							local v`k'1 = r(N)
							qui count if ``=`j'+`nbitems''' == `=`m' - `k''
							local v`k'2 = r(N)
							if (`v`k'1' != 0 | `v`k'2' != 0) & `stop' != 0 {
								qui replace ``j''= `=`m'-`k'' if ``j''==`m'
								qui replace ``=`j'+`nbitems''' =`=`m'-`k'' if ``=`j'+`nbitems''' ==`m'
								di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `=`m'-`k'' merged"
								local stop = 0
							}
						}
					}
					else {
						local stop = 1
						forvalues k = 1/`=`nbm_`j''-`m'' {
							qui count if ``j'' == `=`m' + `k''
							local v`k'1 = r(N)
							qui count if ``=`j'+`nbitems''' == `=`m' + `k''
							local v`k'2 = r(N)
							if (`v`k'1' != 0 | `v`k'2' != 0) & `stop' != 0 {
								qui replace ``j''=`=`m' + `k'' if ``j''==`m'
								qui replace ``=`j'+`nbitems'''=`=`m' + `k'' if ``=`j'+`nbitems'''==`m'
								di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `=`m'+`k'' merged"
								local stop = 0
							}
							else {
								if `stop' != 0 {
									qui replace ``j''= `nbm_`j'' if ``j''==`m'
									qui replace ``=`j'+`nbitems'''= `nbm_`j'' if ``=`j'+`nbitems'''==`m'
									di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `nbm_`j'' merged"
									local stop = 0
								}
							}
						}
					}
				}
			}
		}
	}
}
else {									// Cas où l'option groupe est utilisée
	forvalues j = 1 / `nbitems' {
		local recoda_`j' = 0
		qui tab ``j'' if `gp' == 0, matrow(rect1_g0_`j') matcell(nbrt1_g0_`j')					// Récupération des infos moda du temps 1pour chaque groupe
		local minm`j'_t1_g0 = rect1_g0_`j'[1,1]
		local maxm`j'_t1_g0 = rect1_g0_`j'[r(r),1]

		qui tab ``j'' if `gp' == 1, matrow(rect1_g1_`j') matcell(nbrt1_g1_`j')
		local minm`j'_t1_g1 = rect1_g1_`j'[1,1]
		local maxm`j'_t1_g1 = rect1_g1_`j'[r(r),1]

		qui tab ``=`j'+`nbitems''' if `gp' == 0, matrow(rect2_g0_`j') matcell(nbrt2_g0_`j')				// Récupération des infos moda du temps 2 pour chaque groupe
		local minm`j'_t2_g0 = rect2_g0_`j'[1,1]
		local maxm`j'_t2_g0 = rect2_g0_`j'[r(r),1]

		qui tab ``=`j'+`nbitems''' if `gp' == 1	, matrow(rect2_g1_`j') matcell(nbrt2_g1_`j')
		local minm`j'_t2_g1 = rect2_g0_`j'[1,1]
		local maxm`j'_t2_g1 = rect2_g0_`j'[r(r),1]

		local minm_`j' = min(`minm`j'_t1_g0',`minm`j'_t2_g0',`minm`j'_t1_g1',`minm`j'_t2_g1')		// Info moda pour l'item j
		local maxm_`j' = max(`maxm`j'_t1_g0',`maxm`j'_t2_g0',`maxm`j'_t1_g1',`maxm`j'_t2_g1')
		local nbm_`j' = `=`maxm_`j''-`minm_`j''+1'

		if `minm_`j'' != 0 & `com_z' == 0 {
			local com_z = 1
		}
		//Recodage des réponses en 0, 1, 2, etc...
		forvalues r = 0/`=`maxm_`j''-1'  {
			qui replace ``j'' = `r' if ``j'' == `=`r'+`minm_`j'''
			qui replace ``=`j'+`nbitems''' = `r' if ``=`j'+`nbitems''' == `=`r'+`minm_`j'''
		}

		// Vérif. Que toutes les modas sont utilisées & concordance entre temps
		forvalues m = 0/`=`nbm_`j''-1' {
			qui count if ``j'' == `m' & `gp' == 0
			local nb_rn1_g0 = r(N)
			qui count if ``j'' == `m' & `gp' == 1
			local nb_rn1_g1 = r(N)
			qui count if ``=`j'+`nbitems''' == `m' & `gp' == 0
			local nb_rn2_g0 = r(N)
			qui count if ``=`j'+`nbitems''' == `m' & `gp' == 1
			local nb_rn2_g1 = r(N)
			local nb_rn = min(`nb_rn1_g0',`nb_rn2_g0',`nb_rn1_g1',`nb_rn2_g1')

			if `nb_rn' == 0 {		// Une moda n'est pas utilisée
				local recoda_`j' = 1
				if  `m' == 0 | `m' < `minm`j'_t1_g0' | `m' < `minm`j'_t2_g0' | `m' < `minm`j'_t1_g1' | `m' < `minm`j'_t2_g1' { // La moda 0 n'est pas utilisée
					local stop = 1
					forvalues k = 1/`=`nbm_`j''-`m'' {
						qui count if ``j'' == `=`m' + `k'' & `gp' == 0
						local v`k'1_0 = r(N)
						qui count if ``j'' == `=`m' + `k'' & `gp' == 1
						local v`k'1_1 = r(N)
						qui count if ``=`j'+`nbitems''' == `=`m' + `k'' & `gp' == 0
						local v`k'2_0 = r(N)
						qui count if ``=`j'+`nbitems''' == `=`m' + `k'' & `gp' == 1
						local v`k'2_1 = r(N)
						if (`v`k'1_0' != 0 | `v`k'2_0' != 0 | `v`k'1_1' != 0 | `v`k'2_1' != 0) & `stop' != 0 {
							qui replace ``j''= `=`m'+`k'' if ``j''==`m'
							qui replace ``=`j'+`nbitems'''=`=`m'+`k'' if ``=`j'+`nbitems'''==`m'
							di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `=`m'+`k'' merged"
							local stop = 0
						}
					}
				}
				else if `m' == `=`nbm_`j''-1' | `m' >= `maxm`j'_t2_g0' | `m' >= `maxm`j'_t1_g1' | `m' >= `maxm`j'_t2_g1' {	// La moda max n'est pas utilisée
					local stop = 1
					forvalues k = 1/`=`m'' {
						qui count if ``j'' == `=`m' - `k'' & `gp' == 0
						local v`k'1_0 = r(N)
						qui count if ``j'' == `=`m' - `k'' & `gp' == 1
						local v`k'1_1 = r(N)
						qui count if ``=`j'+`nbitems''' == `=`m' - `k'' & `gp' == 0
						local v`k'2_0 = r(N)
						qui count if ``=`j'+`nbitems''' == `=`m' - `k'' & `gp' == 1
						local v`k'2_1 = r(N)
						if (`v`k'1_0' != 0 | `v`k'2_0' != 0 | `v`k'1_1' != 0 | `v`k'2_1' != 0 ) & `stop' != 0 {
							qui replace ``j''= `=`m' - `k'' if ``j''==`m'
							qui replace ``=`j'+`nbitems'''= `=`m' - `k'' if ``=`j'+`nbitems'''==`m'
							di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `=`m'-`k'' merged"
							local stop = 0
						}
					}
				}
				else {														// Moda central non utilisée
					if runiform()>0.5{											// Tirage au sort pour regrouper
						local stop = 1
						forvalues k = 1/`m' {
							qui count if ``j'' == `=`m' - `k'' & `gp' == 0
							local v`k'1_0 = r(N)
							qui count if ``j'' == `=`m' - `k'' & `gp' == 1
							local v`k'1_1 = r(N)
							qui count if ``=`j'+`nbitems''' == `=`m' - `k'' & `gp' == 0
							local v`k'2_0 = r(N)
							qui count if ``=`j'+`nbitems''' == `=`m' - `k'' & `gp' == 1
							local v`k'2_1 = r(N)
							if (`v`k'1_0' != 0 | `v`k'2_0' != 0 | `v`k'1_1' != 0 | `v`k'2_1' != 0) & `stop' != 0 {
								qui replace ``j''= `=`m'-`k'' if ``j''==`m'
								qui replace ``=`j'+`nbitems'''=`=`m'-`k'' if ``=`j'+`nbitems'''==`m'
								di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `=`m'-`k'' merged"
								local stop = 0
							}
						}
					}
					else {
						local stop = 1
						forvalues k = 1/`=`nbm_`j''-`m'' {
							qui count if ``j'' == `=`m' + `k'' & `gp' == 0
							local v`k'1_0 = r(N)
							qui count if ``j'' == `=`m' + `k'' & `gp' == 1
							local v`k'1_1 = r(N)
							qui count if ``=`j'+`nbitems''' == `=`m' + `k'' & `gp' == 0
							local v`k'2_0 = r(N)
							qui count if ``=`j'+`nbitems''' == `=`m' + `k'' & `gp' == 1
							local v`k'2_1 = r(N)
							if (`v`k'1_0' != 0 | `v`k'2_0' != 0 | `v`k'1_1' != 0 | `v`k'2_1' != 0) & `stop' != 0{
								qui replace ``j''=`=`m' + `k'' if ``j''==`m'
								qui replace ``=`j'+`nbitems'''=`=`m' + `k'' if ``=`j'+`nbitems'''==`m'
								di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `=`m'+`k'' merged"
								local stop = 0
							}
							else {
								if `stop' != 0 {
									qui replace ``j''= `nbm_`j'' if ``j''==`m'
									qui replace ``=`j'+`nbitems'''= `nbm_`j'' if ``=`j'+`nbitems'''==`m'
									di  "WARNING: items ``j'' & ``=`j'+`nbitems''': answers `m' and `nbm_`j'' merged"
									local stop = 0
								}
							}
						}
					}
				}
			}
		}
	}
}

if `com_z' == 1 {
	di
	di "WARNING : Automatic recoding, the first response category is 0. see {help rosali:help rosali}."
	di
}

forvalues j =1/`nbitems' {
		qui tab ``j'', matrow(rec)					// Récupération des infos moda du temps 1
		local nbm`j'_t1 = r(r)

		qui tab ``=`j'+`nbitems''' 					// Récupération des infos moda du temps 2
		local nbm`j'_t2 = r(r)

	local nbm_`j' = max(`nbm`j'_t1', `nbm`j'_t2')
	//Recodage des réponses en 0, 1, 2, etc...
	forvalues r = 0/`=`nbm_`j''-1'  {
		qui replace ``j'' = `r' if ``j'' == `=rec[`=`r'+1',1]'
		qui replace ``=`j'+`nbitems''' = `r' if ``=`j'+`nbitems''' == `=rec[`=`r'+1',1]'
	}
}

/* Calcul de nbmoda & nbdif */
forvalues j = 1/`nbitems' {
	qui tab ``j''
	local nbmoda_`j' = r(r)
	local nbdif_`j' = r(r) - 1
}

local maxdif = 0
local nbmoda_sum = 0
forvalues j = 1/`nbitems' {
	if `maxdif' < `nbdif_`j'' {
		local maxdif = `nbdif_`j''
	}
	local nbmoda_sum = `nbmoda_sum' + `nbdif_`j''
}

/* Au moins 2 moda par item */
forvalues j=1/`nbitems' {
	if `nbmoda_`j'' == 1 {
	di as error "``j'' have only one response category, the analysis can be performed only if each item has at least 2 response categories"
		error 198
		exit
	}
}

local coln ""
forvalues j =1 /`nbitems' {
	local coln "`coln' ``j''"
}

matrix nbmod = J(2,`nbitems',.)

matrix colnames nbmod = `coln'
matrix rownames nbmod = NbModa Recoding

forvalues j = 1/`nbitems' {
	matrix nbmod[1,`j'] = `nbmoda_`j''
	matrix nbmod[2,`j'] = `recoda_`j''
}

*Erreur si plus de 200 difficultés
local nb_test = 0
forvalues j=1/`nbitems' {
	local nb_test = `nb_test'+`nbmoda_`j'' -1
}

if `nb_test' >= 200 {
	di as error "The total number of items difficulties to be estimated must be less than 200 ({hi:moda} option option)."
	error 198
	exit
}

local nbitp = 0

forvalues j = 1/`nbitems' {
	if `nbmoda_`j'' >= 2 {
		local nbitp = `nbitp' + 1
	}
}

qui count
local nbpat = r(N)


/*********************************
* AFFICHAGE INITIAL
*********************************/
di
di _col(5) "{hline 78}"
di  _col(15) "Time 1" _col(42) "Time 2" _col(65) "Nb of Answer Cat."
di _col(5) "{hline 78}"
forvalues j=1/`nbitems' {
    di as text _col(15) abbrev("``j''",20)  _col(42) abbrev("``=`j'+`nbitems'''",20) _col(65) `nbmoda_`j''
}
di _col(5) "{hline 78}"
if "`group'" != "" {
	di _col(10) "Nb of patients: "  abbrev("`gp'",20) " 0 = `nbp_gp0' ;", abbrev("`gp'",20) " 1 = `nbp_gp1'"
	di _col(5) "{hline 78}"
}
else {
	di _col(10) "Nb. of patients: `nbpat'"
	di _col(5) "{hline 78}"
}
di
if `nbitems' == 1 {
	di as error "The analysis can only be performed with at least 2 items."
	error 198
	exit
}
forvalues j = 1/`nbitems' {
	if `nbmoda_`j'' == 2 {
		di "WARNING: ``j'' has only 2 response categories, no distinction can be made between uniform or non-uniform recalibration."
	}
	if `nbmoda_`j'' == 1 {
		di as error "Only `nbmoda_`j'' response categories of item ``j'' were used by the sample, the analysis cannot be performed."
		error 198
		exit
	}
	if `nbmoda_`j'' == 0 {
		di as error "No response categories of item ``j'' were used by the sample, the analysis cannot be performed."
		error 198
		exit
	}
}
di
if "`group'" != "" {
	di _col(2) as text "For all models : - mean of the latent trait in `gp' 0 at time 1 is constrained at 0"
	di _col(19) "- equality of variances between groups"
	di
}
else {
	di _col(2) as text "For all models : mean of the latent trait at time 1 is constrained at 0"
	di
}



/*********************************
* DEFINITION DES CONTRAINTES
*********************************/

if "`group'"!="" { // Contraintes si option groupe
	*EGALITE ENTRE GROUPES A T1 (1-200)
	forvalues j=1/`nbitems'{
		forvalues p=1/`nbdif_`j''{
			constraint `=0+`maxdif'*(`j'-1)+`p'' [`p'.``j'']0bn.`gp'=[`p'.``j'']1.`gp'
		}
	}

	*DIF UNIFORME A T1 (201-400)
	forvalues j=1/`nbitems'{
		forvalues p=2/`nbdif_`j''{
			constraint `=200+`maxdif'*(`j'-1)+`p'' [`p'.``j'']1.`gp'-[`p'.``j'']0bn.`gp'=`p'*[1.``j'']1.`gp'-`p'*[1.``j'']0bn.`gp'
		}
	}

	*EGALITES ENTRE T1 et T2, groupe 0 (401-600)
	forvalues j=1/`nbitems'{
		forvalues p=1/`nbdif_`j''{
			constraint `=400+`maxdif'*(`j'-1)+`p'' [`p'.``j'']0bn.`gp'=[`p'.``=`j'+`nbitems''']0bn.`gp'
		}
	}

	*EGALITES ENTRE T1 et T2, groupe 1 (601-800)
	forvalues j=1/`nbitems'{
		forvalues p=1/`nbdif_`j''{
			constraint `=600+`maxdif'*(`j'-1)+`p'' [`p'.``j'']1.`gp'=[`p'.``=`j'+`nbitems''']1.`gp'
		}
	}

	* RC COMMUNE (801-1000)
	forvalues j=1/`nbitems'{
		forvalues p=1/`nbdif_`j''{
			constraint `=800+`maxdif'*(`j'-1)+`p'' [`p'.``=`j'+`nbitems''']0bn.`gp'-[`p'.``j'']0bn.`gp'=[`p'.``=`j'+`nbitems''']1.`gp'-[`p'.``j'']1.`gp'
		}
	}

	* RC UNIFORME, groupe 0 (1001-1200)
	forvalues j=1/`nbitems'{
		forvalues p=2/`nbdif_`j''{
			constraint `=1000+`maxdif'*(`j'-1)+`p'' `p'*([1.``=`j'+`nbitems''']0bn.`gp'-[1.``j'']0bn.`gp')=[`p'.``=`j'+`nbitems''']0bn.`gp'-[`p'.``j'']0bn.`gp'
		}
	}

	* RC UNIFORME, groupe 1 (1201-1400)
	forvalues j=1/`nbitems'{
		forvalues p=2/`nbdif_`j''{
			constraint `=1200+`maxdif'*(`j'-1)+`p'' `p'*([1.``=`j'+`nbitems''']1.`gp'-[1.``j'']1.`gp')=[`p'.``=`j'+`nbitems''']1.`gp'-[`p'.``j'']1.`gp'
		}
	}

	*Sans interaction temps x groupe
	constraint 1999 [/]:mean(THETA2)#1.`gp'-[/]:mean(THETA2)#0bn.`gp'=[/]:mean(THETA1)#1.`gp'-[/]:mean(THETA1)#0bn.`gp'
}
else { //Contraintes si pas d'option groupe
	*EGALITE ENTRE T1 et T2 (401-600)
	forvalues j=1/`nbitems'{
		forvalues p=1/`nbdif_`j''{
			constraint `=400+`maxdif'*(`j'-1)+`p'' [`p'.``j'']:_cons = [`p'.``=`j'+`nbitems''']:_cons
		}
	}
	*RC UNIFORME (1001-1200)
	forvalues j=1/`nbitems'{
		forvalues p=2/`nbdif_`j''{
			constraint `=1000+`maxdif'*(`j'-1)+`p'' `p'*([1.``=`j'+`nbitems''']:_cons - [1.``j'']:_cons)=[`p'.``=`j'+`nbitems''']:_cons -[`p'.``j'']:_cons
		}
	}
}

/*********************************
* MATRICE DES RESULTATS
*********************************/
matrix dif_rc=J(`nbitems',8,.)
matrix colnames dif_rc=DIFT1 DIFU RC RC_DIF RCG0 RCUG0 RCG1 RCUG1
local rown ""

forvalues j =1 /`nbitems' {
	local rown "`rown' ``j''"
}
matrix rownames dif_rc = `rown'

*Nb modalité max
local nbdif_max = 0
forvalues j=1/`nbitems' {
	if `nbdif_max' < `nbdif_`j'' {
		local nbdif_max = `nbdif_`j''
	}
}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////					PARTIE 1 : DIF A T1 ?						////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

if "`group'"!="" & "`nodif'"=="" { // PARTIE 1 = Slmt si option group & pas de "nodif"

		di _dup(49) "_ "
		di
		*di  as input "PART 1: DETECTION OF DIFFERENCE IN ITEM DIFFICULTIES BETWEEN GROUPS AT TIME 1"

    *********************************
  	**			MODEL B			   **
  	*********************************

  	local model ""
  	forvalues j=1/`nbitems'{
  		forvalues p=1/`nbdif_`j''{
  			local model "`model' (`p'.``j''<-THETA@`p')"
  		}
  	}

  	qui gsem `model', mlogit tol(0.01) iterate(100) group(`gp') ginvariant(coef loading cons) var(0: THETA@v) var(1:THETA@v) latent(THETA) nocapslatent
  	/* Stockage des estimations du modèle */
  	estimates store modeldifB
  	matrix val_mB = r(table)
  	matrix esti_B = e(b)


  	/* Calcul des difficultés d'item (delta_j) */
  	matrix delta_mB=J(`nbitems',`=`nbdif_max'*2',.)
  	local name_partOneC ""

  	forvalues p=1/`nbdif_max' {
  		forvalues g=0/1 {
  			local name_partOneC "`name_partOneC' delta_`p'_gp`g'"
  		}
  	}
  	local name_partOneL ""

  	forvalues j=1/`nbitems' {
  			local name_partOneL "`name_partOneL' ``j''"
  	}

  	matrix colnames delta_mB = `name_partOneC'
  	matrix rownames delta_mB = `name_partOneL'
  		matrix delta_mB_se=J(`nbitems',`=`nbdif_max'*2',.)
  	local name_partOneC_se ""

  	forvalues p=1/`nbdif_max' {
  		forvalues g=0/1 {
  			local name_partOneC_se "`name_partOneC_se' delta_`p'_gp`g'_se"
  		}
  	}

  	matrix colnames delta_mB_se = `name_partOneC_se'
  	matrix rownames delta_mB_se = `name_partOneL'

  	forvalues j=1/`nbitems'{
  		forvalues p=1/`nbdif_`j''{
  			forvalues g=0/1{
  				qui lincom -[`p'.``j'']:`g'.`gp'
  				local delta`j'_`p'g`g'mB=r(estimate)
  				local delta`j'_`p'g`g'mB_se=r(se)
  				if `p'>1{
  					qui lincom [`=`p'-1'.``j'']:`g'.`gp' - [`p'.``j'']:`g'.`gp'
  					local delta`j'_`p'g`g'mB = r(estimate)
  					local delta`j'_`p'g`g'mB_se = r(se)
  				}
  				matrix delta_mB[`j',`=2*`p'-1+`g'']=`delta`j'_`p'g`g'mB'
  				matrix delta_mB_se[`j',`=2*`p'-1+`g'']=`delta`j'_`p'g`g'mB_se'
  			}
  		}
  	}

  	matrix var_mB = (val_mB[1,"/var(THETA)#0bn.`gp'"]\val_mB[2,"/var(THETA)#0bn.`gp'"])

  	/*group effect*/
  	qui lincom [/]:mean(THETA)#1.`gp'-[/]:mean(THETA)#0bn.`gp'
  	local geffmB=r(estimate)
  	local segeffmB=r(se)
  	qui test [/]:mean(THETA)#1.`gp'-[/]:mean(THETA)#0bn.`gp'=0
  	local gcmBp=r(p)
  	local gcmBchi=r(chi2)
  	local gcmBdf=r(df)


  	*********************************
  	**			MODEL A			   **
  	*********************************

  	local model ""
  	forvalues j=1/`nbitems'{
  		forvalues p=1/`nbdif_`j''{
  			local model "`model' (`p'.``j''<-THETA@`p')"
  		}
  	}

  	 qui gsem `model', mlogit tol(0.01) iterate(100) group(`gp') ginvariant(coef loading means) var(0: THETA@v) var(1:THETA@v) from(esti_B, skip) latent(THETA) nocapslatent

  	/* Stockage des estimations du modèle */
  	estimates store modeldifA
  	matrix val_mA = r(table)
  	matrix esti_A = e(b)

  	/* Calcul des difficultés d'item (delta_j) */
  	matrix delta_mA=J(`nbitems',`=`nbdif_max'*2',.)
  	local name_partOneC ""

  	forvalues p=1/`nbdif_max' {
  		forvalues g=0/1 {
  			local name_partOneC "`name_partOneC' delta_`p'_gp`g'"
  		}
  	}
  	local name_partOneL ""

  	forvalues j=1/`nbitems' {
  			local name_partOneL "`name_partOneL' ``j''"
  	}
  	matrix colnames delta_mA = `name_partOneC'
  	matrix rownames delta_mA = `name_partOneL'
  	matrix delta_mA_se=J(`nbitems',`=`nbdif_max'*2',.)
  	local name_partOneC_se ""

  	forvalues p=1/`nbdif_max' {
  		forvalues g=0/1 {
  			local name_partOneC_se "`name_partOneC_se' delta_`p'_gp`g'_se"
  		}
  	}

  	matrix colnames delta_mA_se = `name_partOneC_se'
  	matrix rownames delta_mA_se = `name_partOneL'

  	forvalues j=1/`nbitems'{
  		forvalues p=1/`nbdif_`j''{
  			forvalues g=0/1{
  				qui lincom -[`p'.``j'']:`g'.`gp'
  				local delta`j'_`p'g`g'mA=r(estimate)
  				local delta`j'_`p'g`g'mA_se=r(se)
  				if `p'>1{
  					qui lincom [`=`p'-1'.``j'']:`g'.`gp' - [`p'.``j'']:`g'.`gp'
  					local delta`j'_`p'g`g'mA = r(estimate)
  					local delta`j'_`p'g`g'mA_se = r(se)
  				}
  				matrix delta_mA[`j',`=2*`p'-1+`g'']=`delta`j'_`p'g`g'mA'
  				matrix delta_mA_se[`j',`=2*`p'-1+`g'']=`delta`j'_`p'g`g'mA_se'
  			}
  		}
  	}
  	//Variance et se mA
  	matrix var_mA = (val_mA[1,"/var(THETA)#0bn.`gp'"]\val_mA[2,"/var(THETA)#0bn.`gp'"])

	*********************************
	*************MODEL C*************
	*********************************
	// Etape itérative si lrtest significatif
	local nb_stepC = 0
  qui lrtest modeldifA modeldifB
  local diftestp=r(p)
	if `diftestp'<0.05{  /*If pvalue(LRtest)<0.05 then step C*/
		di
		*di as input "PROCESSING STEP C"
		di

		/*test DIF pour chaque item*/
		local boucle = 1
		local stop = 0
		while `boucle'<=`=`nbitp'-1' & `stop'==0{ /*on s'arrête quand on a libéré du DIF sur (tous les items-1) ou lorsqu'il n'y a plus de tests significatifs*/
			local nb_stepC = `boucle'
			local pajust=0.05/`=`nbitp'+1-`boucle''
			/*réinitialisation de la matrice de test*/
			matrix test_difu_`boucle'=J(`nbitems',3,.)
			matrix colnames test_difu_`boucle'=chi_DIFU df_DIFU pvalueDIFU
			matrix test_dif_`boucle'=J(`nbitems',3,.)
			matrix colnames test_dif_`boucle'=chi_DIF df_DIF pvalueDIF
			local nbsig=0
			local minpval=1
			local itemdif=0
			if "`detail'" != ""{

				di as text "Loop `boucle'"
				di as text _col(5) "Adjusted  alpha: " %6.4f  `pajust'
				di
				di as text _col(10) "{hline 65}"
				di as text _col(10) "Freed item" _col(31) "Chi-Square"  _col(48) "DF" _col(57) "P-Value"
				di as text _col(10) "{hline 65}"
			}
			/*boucle de test*/
			forvalues j=1/`nbitems'{
				//if `nbdif_`j'' > 2 {
					local model ""
					local listconst ""
					if dif_rc[`j',1]==. | dif_rc[`j',1]==0 { /*si pas de DIF déjà détecté sur l'item j*/
						/*on libère le DIF de l'item i: pas de contraintes*/
						forvalues k=1/`nbitems'{ /*contraintes pour les autres items (si DIF NU sur item k, pas de contraintes*/
							if `k'!=`j' & `nbmoda_`j'' >= 2 {
								if dif_rc[`k',1]==. | dif_rc[`k',1]==0 {/*pas de DIF sur item k: contraintes 1-200*/
									forvalues p=1/`nbdif_`k''{
										qui local listconst "`listconst' `=0+`maxdif'*(`k'-1)+`p''"
										qui constraint list `=0+`maxdif'*(`k'-1)+`p''
									}
								}
								else{
									if dif_rc[`k',2]!=. & dif_rc[`k',2]!= 0 & `nbmoda_`k'' > 2 { /*DIF U: contraintes 201-400*/
										forvalues p=2/`nbdif_`k''{
											qui local listconst "`listconst' `=200+`maxdif'*(`k'-1)+`p''"
											qui constraint list `=200+`maxdif'*(`k'-1)+`p''
										}
									}
								}
							}
						}
						forvalues jj=1/`nbitems'{
							forvalues p=1/`nbdif_`jj''{
								local model "`model' (`p'.``jj''<-THETA@`p')"
							}
						}
						qui gsem `model', mlogit tol(0.01) iterate(100) group(`gp') ginvariant(coef loading) var(0: THETA@v) var(1:THETA@v) constraint(`listconst') from(esti_B) latent(THETA) nocapslatent
						estimates store modeldif3b`boucle'it`i'

						*************************
						*****test DIF item i*****
						*************************
						qui test [1.``j'']1.`gp'=[1.``j'']0bn.`gp'
						if `nbmoda_`j'' > 2 {
							forvalues p=2/`nbdif_`j''{
								qui test [`p'.``j'']1.`gp'=[`p'.``j'']0bn.`gp', acc
							}
						}
						matrix test_dif_`boucle'[`j',1]=(r(chi2),r(df),r(p))

						/* Test DIF Uniforme */
						if `nbmoda_`j'' > 2 {
							qui test 2*([1.``j'']1.`gp'-[1.``j'']0bn.`gp')=[2.``j'']1.`gp'-[2.``j'']0bn.`gp'
							forvalues p=3/`nbdif_`j''{
								qui test `p'*([1.``j'']1.`gp'-[1.``j'']0bn.`gp')=[`p'.``j'']1.`gp'-[`p'.``j'']0bn.`gp', acc
							}
							matrix test_difu_`boucle'[`j',1]=(r(chi2), r(df), r(p))
						}

						if test_dif_`boucle'[`j',3]<`pajust'{/*si DIF sur item i*/
							local ++nbsig
							if test_dif_`boucle'[`j',3]<`minpval'{
								local minpval=test_dif_`boucle'[`j',3]
								local itemdif=`j'
							}
						}
						if "`detail'" != "" {
							di as text  _col(10) abbrev("``j''",15)  as result _col(31) %6.3f test_dif_`boucle'[`j',1] _col(48) test_dif_`boucle'[`j',2] _col(57)  %6.4f test_dif_`boucle'[`j',3]
						}
					}
			}
			/*si nb de tests significatifs=0, on arrête*/
			if `nbsig'==0{
				local stop=1
				if `boucle' == 1 {
					if "`detail'" != "" {
						di as text _col(10) "{hline 65}"
						di
						di as result  "No significant test: no difference between groups detected, no DIF detected"
						di
					}
				}
				else {
					if "`detail'" != ""{
						di as text _col(10) "{hline 65}"
						di
						di as result "No other significant tests"
						di
					}
				}
			}
			else{/*si nb de tests significatifs>0, mise à jour de la matrice de résultats*/
				matrix dif_rc[`itemdif',1]=`boucle'
				if "`detail'" != ""{
					di as text _col(10) "{hline 65}"
					di
					di as result "Difference between groups on ``itemdif'' at time 1"
				}
				if `nbmoda_`itemdif'' > 2 {
					if "`detail'" != "" {

						di
						di %~60s as text "Test of uniform difference"
						di _col(10) "{hline 40}"
						di _col(10) as text "Chi-square" _col(28) "DF" _col(40) "P-value"
						di _col(10) as result %4.2f `=test_difu_`boucle'[`itemdif',1]' _col(28) `=test_difu_`boucle'[`itemdif',2]' _col(40) %4.2f `=test_difu_`boucle'[`itemdif',3]'
						di _col(10) as text "{hline 40}"
					}
					if test_difu_`boucle'[`itemdif',3]<0.05{	/*DIF NU détectée*/
						matrix dif_rc[`itemdif',2]=0
							di
							di as result "``itemdif'' : Non-uniform differences of item difficulties between groups at T1"
							di
					}
					else{/*DIF U détectée*/
						matrix dif_rc[`itemdif',2]=`boucle'
							di
							di as result "``itemdif'' : Uniform differences of item difficulties between groups at T1"
							di
					}
				}
				else {
					// Différence entre groupes au temps 1 mais slmt 2 moda. donc pas de U ou NU
					di _col(15) _dup(60) "-"
				}
			}
			local ++boucle
		}
	}

	/* MODELE FINAL DE LA PARTIE 1. Si DIFT1 détecté (=Au moins 2 boucles dans l'étape C)*/
	if `nb_stepC' > 1 {
		forvalues j=1/`nbitems'{
			local model ""
			local listconst ""
			if dif_rc[`j',1]==. | dif_rc[`j',1]==0 { /*si pas de DIF: contraintes 1-200*/
				forvalues p=1/`nbdif_`j''{
					qui local listconst "`listconst' `=0+`maxdif'*(`j'-1)+`p''"
					qui constraint list `=0+`maxdif'*(`j'-1)+`p''
				}
			}
			else {
				if dif_rc[`j',2]!=. & dif_rc[`j',2]!=0 { /*DIF U: contraintes 201-400*/
					forvalues p=2/`nbdif_`j''{
						qui local listconst "`listconst' `=200+`maxdif'*(`j'-1)+`p''"
						qui constraint list `=200+`maxdif'*(`j'-1)+`p''
					}
				}
			}
		}
		forvalues j=1/`nbitems'{
			forvalues p=1/`nbdif_`j''{
				local model "`model' (`p'.``j''<-THETA@`p')"
			}
		}

		qui gsem `model', mlogit tol(0.01) iterate(100) group(`gp') ginvariant(coef loading) var(0: THETA@v) var(1:THETA@v) constraint(`listconst') from(esti_B) latent(THETA) nocapslatent
		/* Stockage des estimations du modèle */
		estimates store modeldifCFin
		matrix val_mC = r(table)

		/* Calcul des difficultés d'item (delta_j) */
		matrix delta_mCFin=J(`nbitems',`=`nbdif_max'*2',.)
		local name_partOneC ""
		forvalues p=1/`nbdif_max' {
			forvalues g=0/1 {
				local name_partOneC "`name_partOneC' delta_`p'_gp`g'"
			}
		}
		local name_partOneL ""
		forvalues j=1/`nbitems' {
			local name_partOneL "`name_partOneL' ``j''"
		}
		matrix colnames delta_mCFin = `name_partOneC'
		matrix rownames delta_mCFin = `name_partOneL'

		matrix delta_mCFin_se=J(`nbitems',`=`nbdif_max'*2',.)
		local name_partOneC_se ""

		forvalues p=1/`nbdif_max' {
			forvalues g=0/1 {
				local name_partOneC_se "`name_partOneC_se' delta_`p'_gp`g'_se"
			}
		}
		matrix colnames delta_mCFin_se = `name_partOneC_se'
		matrix rownames delta_mCFin_se = `name_partOneL'

		forvalues j=1/`nbitems'{
			forvalues p=1/`nbdif_`j''{
				forvalues g=0/1{
					qui lincom -[`p'.``j'']:`g'.`gp'
					local delta`j'_`p'g`g'mCFin=r(estimate)
					local delta`j'_`p'g`g'mCFin_se=r(se)
					if `p'>1{
						qui lincom [`=`p'-1'.``j'']:`g'.`gp' - [`p'.``j'']:`g'.`gp'
						local delta`j'_`p'g`g'mCFin = r(estimate)
						local delta`j'_`p'g`g'mCFin_se = r(se)
					}
					matrix delta_mCFin[`j',`=2*`p'-1+`g'']=`delta`j'_`p'g`g'mCFin'
					matrix delta_mCFin_se[`j',`=2*`p'-1+`g'']=`delta`j'_`p'g`g'mCFin_se'
				}
			}
		}
		if "`group'" != "" { //Variance et se mA
			matrix var_mC = (val_mC[1,"/var(THETA)#0bn.`gp'"]\val_mC[2,"/var(THETA)#0bn.`gp'"])
		}
			/*group effect*/
		qui lincom [/]:mean(THETA)#1.`gp'-[/]:mean(THETA)#0bn.`gp'
		local geffmCFin=r(estimate)
		local segeffmCFin=r(se)
		qui	test [/]:mean(THETA)#1.`gp'-[/]:mean(THETA)#0bn.`gp'=0
		local gcmCFinp=r(p)
		local gcmCFinchi=r(chi2)
		local gcmCFindf=r(df)
	}
}

	*********************************
	 ***		BILAN  			***
	*********************************

	if "`group'" != "" & "`nodif'" == "" {
	di
	di %~84s as result "SUMMARY"
	di as result _col(2)  "{hline 80}"
	di as result 		 _col(18) "Difference in"
	di as result _col(2) "Item" _col(18) "groups at T1" _col(36) "Recalibration"  _col(54) "RC " abbrev("`gp'",10) " 0"  _col(72) "RC " abbrev("`gp'",10) " 1"
	di as result _col(2) "{hline 80}"
	forvalues j=1/`nbitems' {
		local RC
		local RCg0
		local RCg1
		local difft1
		if (dif_rc[`j',3] != . & dif_rc[`j',3] != 0 & dif_rc[`j',4] == 0) {
			local RC "Common"
		}
		if (dif_rc[`j',3] != . & dif_rc[`j',3] != 0 & dif_rc[`j',4] != 0) {
			local RC "Differential"
		}
		if `nbmoda_`j'' > 2 {
			if (dif_rc[`j',6]!=. & dif_rc[`j',6] != 0) {
				local RCg0 "Uniform"
			}
			if (dif_rc[`j',6] == 0) {
				local RCg0 "Non-uniform"
			}
				if (dif_rc[`j',8]!=. & dif_rc[`j',8] != 0) {
				local RCg1 "Uniform"
			}
			if ( dif_rc[`j',8] == 0) {
				local RCg1 "Non-uniform"
			}
			if (dif_rc[`j',1] != . ) {
				if (dif_rc[`j',2]!=0) {
					local difft1 "Uniform"
				}
				else {
					local difft1 "Non-uniform"
				}
			}
		}
		else {
			if dif_rc[`j',6] != . {
				local RCg0 " X "
			}
			if dif_rc[`j',8] != . {
				local RCg1 " X "
			}
			if dif_rc[`j',1] != . {
				local difft1 " X "
			}
		}
		di as result _col(2) abbrev("``j''",15) as text _col(18) "`difft1'"  _col(36) "`RC'"  _col(54) "`RCg0'"  _col(72) "`RCg1'"
	}
	di as result _col(2) "{hline 80}"
	di
}
else if "`group'" != "" & "`nodif'" != "" {
	di
	di %~90s as result "SUMMARY"
	di as result _col(10)  "{hline 70}"
	di as result _col(10) "Item" _col(26) "Recalibration"  _col(46) "RC `gp' 0"  _col(62) "RC `gp' 1"
	di _col(10) "{hline 70}"
	forvalues j=1/`nbitems' {
		local RC
		local RCg0
		local RCg1
		if (dif_rc[`j',3] != . & dif_rc[`j',3] != 0 & dif_rc[`j',4] == 0) {
			local RC "Common"
		}
		if (dif_rc[`j',3] != . & dif_rc[`j',3] != 0 & dif_rc[`j',4] != 0) {
			local RC "Differential"
		}
		if `nbmoda_`j'' > 2 {
			if (dif_rc[`j',6]!=. & dif_rc[`j',6] != 0) {
				local RCg0 "Uniform"
			}
			if (dif_rc[`j',6] == 0) {
				local RCg0 "Non-uniform"
			}
			if (dif_rc[`j',8]!=. & dif_rc[`j',8] != 0) {
				local RCg1 "Uniform"
			}
			if ( dif_rc[`j',8] == 0) {
				local RCg1 "Non-uniform"
			}
		}
		else {
			if dif_rc[`j',6] != . {
				local RCg0 " X "
			}
			if dif_rc[`j',8] != . {
				local RCg1 " X "
			}
		}
		di as result _col(10) "``j''" as text  _col(26) "`RC'"  _col(44) "`RCg0'"  _col(62) "`RCg1'"
	}
	di as result _col(10) "{hline 70}"
}
else if "`group'" == "" {
	di
	di %~60s as result "SUMMARY"
	di as result _col(10)  "{hline 40}"
	di _col(10) "Item" _col(36) "Recalibration"
	di _col(10) "{hline 40}"
	forvalues j=1/`nbitems' {
		local RC
		if dif_rc[`j',3] != . {
			if `nbmoda_`j'' > 2 {
				if (dif_rc[`j',6]!=. & dif_rc[`j',6] != 0) {
					local RC "Uniform"
				}
				if (dif_rc[`j',6] == 0) {
					local RC "Non-uniform"
				}
			}
			else {
				local RC " X "
			}
		}
		di as result _col(10) "``j''" as text  _col(38) "`RC'"
	}
	di as result _col(10) "{hline 40}"
	di
}

matrix dif_detect = J(1,`nbitems',.)
local numdif=1
forvalues j=1/`nbitems' {
  if dif_rc[`j',1] != . {
    matrix dif_detect[1,`numdif']=`j'
    local numdif = `numdif'+1
  }
}
return matrix difitems = dif_detect

capture qui use `saverspcm', clear

end