*! Version 4.2 21September2019
*! Jean-Benoit Hardouin, Myriam Blanchin
************************************************************************************************************
* Stata program : pcm
* Estimate the parameters of the Partial Credit Model
* Version 1 : December 17, 2007 [Jean-Benoit Hardouin]
* Version 2 : July 15, 2011 [Jean-Benoit Hardouin]
* Version 2.1 : October 18th, 2011 [Jean-Benoit Hardouin] : -fixedvar- option, new presentation
* Version 2.2 : October 23rd, 2013 [Jean-Benoit Hardouin] : correction of -fixedvar- option
* Version 2.3 : April 10th, 2014 [Jean-Benoit Hardouin] : correction of -fixedvar- option
* Version 2.3 : April 10th, 2014 [Jean-Benoit Hardouin] : correction of -fixedvar- option
* Version 3 : July 6th, 2019 [Jean-Benoit Hardouin] : New version using gsem
* Version 3.1 : July 9th, 2019 [Jean-Benoit Hardouin] : Small corrections
* Version 3.2 : July 17th, 2019 [Jean-Benoit Hardouin] : Small corrections
* Version 3.3 : July 25th, 2019 [Jean-Benoit Hardouin] : -pce- option
* Version 3.4 : August 23th, 2019 [Jean-Benoit Hardouin] : Correction of a bug
* Version 3.5 : August 29th, 2019 [Jean-Benoit Hardouin] : Correction of a bug with modamax``i''
* Version 4: September 13th, 2019 [Myriam Blanchin]: addition of longitudinal pcm
* Version 4.1: September 15th, 2019 [Jean-Benoit Hardouin]: correction of a small bug in the outputs
* Version 4.2: September 27th, 2019 [Jean-Benoit Hardouin]
*
*
* Jean-benoit Hardouin, Myriam Blanchin - University of Nantes - France
* INSERM UMR 1246-SPHERE "Methods in Patient Centered Outcomes and Health Research", Nantes University, University of Tours
* jean-benoit.hardouin@univ-nantes.fr myriam.blanchin@univ-nantes.fr
*
* News about this program : http://www.anaqol.org
*
* Copyright 2007, 2011, 2013, 2014, 2019 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 pcm, rclass 
syntax varlist(min=2 numeric) [iweight] [, CONTinuous(varlist) CATegorical(varlist) ITerate(int 100) TOLerance(real 0.01) model DIFFiculties(name) VARiance(real -1) rsm Graphs noGRAPHItems filesave dirsave(string) docx(string) extension(string) PCE WMLiterate(int 1) GENLT(string) GENINF(string) REPlace postpce visit(varname) id(varname) eqset1(varlist) eqset2(varlist) EQGraph minsize(int 30)]

preserve

/*************************************************************************************************************
QUELQUES TESTS
*************************************************************************************************************/

if `variance'!=-1&`variance'<=0 {
	di in red "The -variance- option cannot be negative"
	exit 198
}
if `variance'!=-1&"`visit'"!="" {
	di in red "The -variance- and -visit- options cannot be used simultaneously."*
	exit 198
}
if "`genlt'"!=""|"`geninf'"!="" {
    capture confirm new variable `genlt' se`genlt' `geninf'
	if _rc!=0&"`replace'"=="" {
	    di in red "The variables `genlt', se`genlt' and/or `geninf' alreday exist. Please modify the -genlt- and/or -geninf- option"
		exit 198
	}
	if _rc!=0&"`replace'"!="" {
	    qui drop `genlt' se`genlt' `geninf' 
	}
}
if ("`eqset1'"!=""&"`eqset2'"=="")|("`eqset1'"==""&"`eqset2'"!="") {
    di in red "The two options -eqset1- and -eqset2- must be used simultaneously"
	exit 198
}
if ("`eqset1'"!=""&"`graphs'"!="") {
    di in red "The two options -eqset1- and -graph- cannot be used simultaneously"
	exit 198
}




/*************************************************************************************************************
GESTION DES VARIABLES CONTINUES ET CATEGORIELLES
*************************************************************************************************************/
if "`visit'"!=""{
	if "`id'"==""{
		di in red "Option -visit- must be combined with option -id-. Please fill in the -id- option"
		exit 198
	}
	qui levelsof `visit'
	local levelsofv `r(levels)'
	local nbvisits=r(r)
	local timemin: word 1 of `levelsofv'
	local timemax: word `nbvisits' of `levelsofv'
	if `timemax'>5{
		di as error "You must use a discrete time variable (-visit- option) with less than 5 measurement occasions"
		error 198
	}
	if `timemin'!=1{
		di as error "You must use a -visit- variable coded at 1 for the first visit"
		error 198
	}
	qui reshape wide `varlist', i(`id') j(`visit')
	local multivist=1
}
else {
	local timemax=1
	foreach i in `varlist' {
	   *rename `i' `i'1
	}
	local multivisit
}
qui count
local nbobs=r(N)

local timelist
forvalues t=1/`timemax'{
	local timelist `timelist' T`t'
}

local modcont
local premodcont
local nbpar=0
local nbcont=0
local nbcat=0
if "`continuous'"!="" {
	tokenize `continuous'
	local nbcont : word count `continuous'
	local continuous
	forvalues i=1/`nbcont' {
		local cont`i' ``i''
		local continuous `continuous' ``i''
		local modcont `modcont' ``i''
	    local ++nbpar
	}
	local premodcont (`modcont'->T1)
	local modcont (`modcont'->`timelist')
}

local modcat
local premodcat
if "`categorical'"!="" {
	tokenize `categorical'
	local nbcat : word count `categorical'
	local categorical
	forvalues i=1/`nbcat' {
		local cat`i' ``i''
		local categorical `categorical' ``i''
		local modcat `modcat' i.``i''
		qui levelsof ``i''
		local levelsof``i'' `r(levels)'
		local nbpar=`nbpar'+`r(r)'-1
	    *di "categorical : ``i'' levels : `levelsof``i'''"
	}
	local premodcat (`modcat'->T1)
	local modcat (`modcat'->`timelist')
}


if "`dirsave'"=="" {
   local dirsave `c(pwd)'
}

/*************************************************************************************************************
GESTION DES ITEMS ET TESTS
*************************************************************************************************************/

tokenize `varlist'
local nbitems : word count `varlist'
marksample touse ,novarlist
*preserve

local modamax=1
local modamin=0
local pbmin
local nbdiff=0
forvalues i=1/`nbitems' {
	local modamax`i'=1
	local modamax``i''=1
	if `timemax'>1 {
		forvalues t=1/`timemax'{
			 *qui replace ``i'`t''=``i'`t''-`min' 
			 qui su ``i''`t'
			 if `r(min)'!=`modamin' {
				 local modamin=r(min)
				 local pbmin `pbmin' ``i'`t''
			 }
			 if `r(max)'>`modamax' {
				 local modamax=r(max)
			 }
			 if `r(max)'>`modamax`i'' {
				local modamax`i'=r(max)
			}
		}
	}
	else {
		qui su ``i''
		if `r(min)'!=`modamin' {
		    local modamin=r(min)
			local pbmin `pbmin' ``i''
		}
		if `r(max)'>`modamax' {
			local modamax=r(max)
		}
		if `r(max)'>`modamax`i'' {
			local modamax`i'=r(max)
			local modamax``i''=r(max)
		}
	}
	if "`rsm'"=="" {
		local nbdiff=`nbdiff'+`modamax`i''
	}
}
if "`rsm'"!="" {
   local nbdiff=`nbitems'+`modamax'-1
}
if `modamin'!=0 {
   di as error "The minimal answer category of each item must be coded by 0. This is not the case for the following items: `pbmin' (`modamin') "
   error 198
}
qui count
local nbind=r(N)
*set trace on
local code
local precode
if `timemax'>1 {
	forvalues k=1/`modamax' {
	   forvalues t=1/`timemax'{
			local code`k'
		   forvalues i=1/`nbitems' {
			   if `k'<=`modamax`i'' {
				  local code`k' `code`k'' `k'.``i''`t'
			   }
		   }
		   local code`k' (`code`k''<-T`t'@`k')
		   if `t'==1{
				local precode `precode' `code`k''
		   }
		   local code `code' `code`k''
		}
	}
}
else {
	forvalues k=1/`modamax' {
		local code`k'
		forvalues i=1/`nbitems' {
			if `k'<=`modamax`i'' {
				local code`k' `code`k'' `k'.``i''
			}
		}
		local code`k' (`code`k''<-T@`k')
		local precode `precode' `code`k''
		local code `code' `code`k''
	}
}

/*************************************************************************************************************
OPTION PCE
*************************************************************************************************************/

if "`pce'"!=""&"`difficulties'"==""&"`visit'"=="" {
    tempname sedelta b
	qui raschpce `varlist'
	local ll=r(ll)
	matrix `sedelta'=r(sedelta)
	matrix `sedelta'=`sedelta''
	matrix `b'=r(b)
	*matrix `b'=`b''
	local difficulties `b'
	*matrix list `b'
	matrix loulou=`b'
	return matrix diff_parm=`b'
	`qui' pcm `varlist', diff(loulou)  geninf(TInf_0) genlt(lt_0) /*postpce*/
	*exit
}


/*************************************************************************************************************
RECUPERATION DES PARAMETRES DE DIFFICULTES ET DEFINITION DES CONTRAINTES
*************************************************************************************************************/
if "`difficulties'"!=""&"`rsm'"!="" {
   di as error "You can not defined in the same time the difficulties and the rsm options"
   error 198
}
local t=1
local constraints
local codemean
local codevar
local codecov
forvalues j=2/`timemax'{
	forvalues i=1/`nbitems' {
		forvalues k=1/`modamax`i'' {
			qui constraint `t' [`k'.``i''`j']_cons=[`k'.``i''`multivisit']_cons
			local constraints `constraints' `t'
			local ++t
			
		}
	}	
	if "`continuous'"=="" &  "`categorical'"==""{
		local codemean `codemean' T`j'@m`j'
		local codevar `codevar' T`j'@v`j'
		forvalues l=1/`=`j'-1'{
			local codecov `codecov' T`l'*T`j'@cov`l'`j'
		}
	}
	else{
		local codevar `codevar' e.T`j'@v`j'
		forvalues l=1/`=`j'-1'{
			local codecov `codecov' e.T`l'*e.T`j'@cov`l'`j'
		}	
	}
}
if `timemax'>1{
	if "`continuous'"=="" &  "`categorical'"==""{	
		local codelg means(T1@0 `codemean') var(T1@v1 `codevar') cov(`codecov')
	}
	else{
		local codelg var(e.T1@v1 `codevar') cov(`codecov')
	}
}
else {
    local constrvar
    if `variance'>0 {
		if "`continuous'"=="" &  "`categorical'"==""{	
			local constrvar var(T1@`variance') 
		}
		else{
			local constrvar var(e.T1@`variance') 
		}
	}
}

if "`difficulties'"!="" {
    tempname beta
	matrix `beta'=J(`nbitems',`modamax',.)
	matrix list `difficulties'
	forvalues i=1/`nbitems' {
		forvalues k=1/`modamax`i'' {
			if `difficulties'[`i',`k']==. {
				 di as error "The kth difficulty parameter of the item ``i'' is not correctly defined in the difficulties matrix"
				 error 198
			}
			else {
				if `k'==1 {
					matrix `beta'[`i',1]=-`difficulties'[`i',1]
  			    }
				else {
				    matrix `beta'[`i',`k']=`beta'[`i',`=`k'-1']-`difficulties'[`i',`k']
				}
				qui constraint `t' [`k'.``i''`multivisit']_cons=`beta'[`i',`k']
				local constraints `constraints' `t'
				local ++t
			}
		}
	}
}



/*************************************************************************************************************
DEFINITION DES CONTRAINTES POUR UN RSM
*************************************************************************************************************/
if "`rsm'"!="" {
   local constraints
   forvalues k=2/`modamax' {
       forvalues i=2/`nbitems' {
	       qui constraint `t'   [`=`k'-1'.``i''`multivisit']_cons-[`k'.``i''`multivisit']_cons+[1.``i''`multivisit']_cons=[`=`k'-1'.`1'1]_cons-[`k'.`1'1]_cons+[1.`1'1]_cons
		   local constraints `constraints' `t'
		   local ++t
	   }
   }
}

/*************************************************************************************************************
MODELE
*************************************************************************************************************/


discard
*di "`qui' gsem `code' `modcont' `modcat' ,iterate(`iterate') tol(`tolerance') constraint(`constraints') latent(`timelist') `codelg' "
if "`model'"!="" {
    local qui
}
else {
	local qui qui
}
if `timemax'==1{
   `qui' gsem `code' `modcont' `modcat' ,iterate(`iterate') tol(`tolerance') constraint(`constraints') latent(`timelist') `constrvar'
}
else{
    di "`qui' gsem `precode' `premodcont' `premodcat',iterate(`iterate') tol(`tolerance') "
    `qui' gsem `precode' `premodcont' `premodcat',iterate(`iterate') tol(`tolerance') constraint(`constraints') 
     matrix esti_B = e(b)
    di "`qui' gsem `code' `modcont' `modcat' ,iterate(`iterate') tol(`tolerance') constraint(`constraints') latent(`timelist') `codelg' from(esti_B,skip)"
    `qui' gsem `code' `modcont' `modcat' ,iterate(`iterate') tol(`tolerance') constraint(`constraints') latent(`timelist') `codelg' from(esti_B,skip)
}
local ll=e(ll)

*set trace on
tempvar latent score group selatent latent2 miss
tempname groups
*capture qui predict mu, mu
*su mu 
qui predict `latent'*,latent se(`selatent'*)

if "`genlt'"!="" {
	if `timemax'==1 {
		qui gen `genlt'=`latent'`i'
		qui gen se`genlt'=`selatent'`i'
	}
	forvalues t=2/`timemax' {
		qui gen `genlt'`t'=`latent'`t'
		qui gen se`genlt'`t'=`selatent'`t'
	}
}

set seed 123456
if `timemax'>1 {
	forvalues t=1/`timemax'{
		qui gen `latent2'`t'=`latent'`t'+invnorm(uniform())*`selatent'`t'
		local listit 
		forvalues i=1/`nbitems' {
			local listit `listit' ``i''`t'
		}
		qui genscore `listit',score(`score'`t')
		qui gengroup `latent'`t',newvariable(`group'`t') continuous minsize(`minsize')
	}
}
else {
	qui gen `latent2'=`latent'+invnorm(uniform())*`selatent'
	local listit 
	forvalues i=1/`nbitems' {
		local listit `listit' ``i''
	}
	qui genscore `listit',score(`score')
	qui gengroup `latent',newvariable(`group') continuous minsize(`minsize')
}

/*time 1 only*/
qui levelsof `group'`multivisit'
local nbgroups=r(r)
matrix `groups'=J(`nbgroups',`=`nbitems'+3',.)
forvalues g=1/`nbgroups' {
	matrix `groups'[`g',`=`nbitems'+3']=0
	forvalues i=1/`nbitems' {
			qui count if ``i''`multivisit'!=.&`group'`multivisit'==`g'
			local n=r(N)
			if `n'>0 {
				qui su ``i''`multivisit' if `group'`multivisit'==`g'
				matrix `groups'[`g',`i']=r(mean)
				matrix `groups'[`g',`=`nbitems'+3']=`groups'[`g',`=`nbitems'+3']+`r(mean)'
			}
			else {
				matrix `groups'[`g',`i']=.
				matrix `groups'[`g',`=`nbitems'+3']=.
			}		
	}
	qui su `latent'`multivisit' if `group'`multivisit'==`g'
	matrix `groups'[`g',`=`nbitems'+1']=r(mean)
	qui count if `group'`multivisit'==`g'
	matrix `groups'[`g',`=`nbitems'+2']=r(N)
}

/*number of non-missing on all time points*/
egen `miss'=rowmiss(`score'*)
qui count if `miss'==0
local nbobsssmd=r(N)
drop `miss'	

di
di as text "Number of individuals:" %6.0f as result `nbobs'
di as text "Number of complete individuals:" %6.0f as result `nbobsssmd'
di as text "Number of items:" %6.0f as result `nbitems'

di as text "Marginal log-likelihood:" %12.4f as result `ll'
di 
return scalar ll=`ll'



*set trace on
/*************************************************************************************************************
RECUPERATION DES ESTIMATIONS DES PARAMETRES DE DIFFICULTE
*************************************************************************************************************/

tempname diff diffmat vardiff diffmat2
*set trace on
qui matrix `diffmat'=J(`nbitems',`modamax',.)
qui matrix `diffmat2'=J(`nbitems',`modamax',.)
qui matrix `diff'=J(`nbdiff',6,.)
local rn
*qui matrix `vardiff'=J(`nbdiff',`nbdiff',.)
*matrix list `diff'
*set trace on
local t=1
forvalues i=1/`nbitems' {
    qui matrix `diffmat'[`i',1]=-_b[1.``i''`multivisit':_cons]
    qui matrix `diffmat2'[`i',1]=-_b[1.``i''`multivisit':_cons]
	qui lincom -_b[1.``i''`multivisit':_cons]
    qui matrix `diff'[`t',1]=`r(estimate)'
    qui matrix `diff'[`t',2]=`r(se)'
    qui matrix `diff'[`t',3]=`r(z)'
    qui matrix `diff'[`t',4]=`r(p)'
    qui matrix `diff'[`t',5]=`r(lb)'
    qui matrix `diff'[`t',6]=`r(ub)'
	local rn `rn' 1.``i''`multivisit'
	local ++t
	local sum _b[1.``i''`multivisit':_cons]
	if "`rsm'"=="" {
	    forvalues k=2/`modamax`i'' {
			local sum "_b[`k'.``i''`multivisit':_cons]-(`sum')"
			*di "``i''`multivisit' `k' `sum'"
			local sum2 "_b[`=`k'-1'.``i''`multivisit':_cons]-_b[`k'.``i''`multivisit':_cons]"
			qui lincom (`sum2')
			*set trace on
			qui matrix `diffmat'[`i',`k']=`r(estimate)'
			qui matrix `diffmat2'[`i',`k']=`diffmat2'[`i',`=`k'-1']+`diffmat'[`i',`k']
			qui matrix `diff'[`t',1]=`r(estimate)'
			qui matrix `diff'[`t',2]=`r(se)'
			qui matrix `diff'[`t',3]=`r(z)'
			qui matrix `diff'[`t',4]=`r(p)'
			qui matrix `diff'[`t',5]=`r(lb)'
			qui matrix `diff'[`t',6]=`r(ub)'
			*qui matrix `vardiff'[`t',`t']=`r(se)'^2
			*set trace off
			local rn `rn' `k'.``i''`multivisit'
			local ++t
		}
    }
}
if "`rsm'"!="" {
    forvalues k=2/`modamax' {
		qui lincom _b[`=`k'-1'.`1'`multivisit':_cons]-_b[`k'.`1'`multivisit':_cons]+_b[1.`1'`multivisit':_cons] /*``i'' instead of `i'?*/
		qui matrix `diff'[`t',1]=`r(estimate)'
		qui matrix `diff'[`t',2]=`r(se)'
		qui matrix `diff'[`t',3]=`r(z)'
		qui matrix `diff'[`t',4]=`r(p)'
		qui matrix `diff'[`t',5]=`r(lb)'
		qui matrix `diff'[`t',6]=`r(ub)'
		forvalues i=1/`nbitems' {
		    qui matrix `diffmat'[`i',`k']=`diff'[`t',1]+`diffmat'[`i',1]
		    qui matrix `diffmat2'[`i',`k']=`diffmat'[`i',`k']+`diffmat2'[`i',`=`k'-1']
		}
		local rn `rn' tau`k'
		local ++t
	}
}
local cn Estimate S.e. z p "Lower bound" "Upper Bound"
matrix colnames `diff'=`cn'
matrix rownames `diff'=`rn'
*matrix list `diff'
*matrix list `diffmat'
*matrix list `diffmat2'
*matrix list `vardiff'

/*************************************************************************************************************
RECUPERATION DES ESTIMATIONS DES PARAMETRES POUR LES COVARIABLES, MOYENNES ET VARIANCES
*************************************************************************************************************/
tempname covariates
local nbcov=0
forvalues j=2/`timemax'{
	local nbcov=`nbcov'+`j'-1
}
qui matrix `covariates'=J(`=`nbpar'+`timemax'+2*`nbcov'',6,.)

*set trace on
local t=1


forvalues j=1/`=`timemax'-1'{
	forvalues k=`=`j'+1'/`timemax'{
		if "`categorical'"=="" & "`continuous'"=="" {
			if `j'==1{
				qui lincom [/]mean(T`k')
			}
			else{
				qui lincom [/]mean(T`k')-[/]mean(T`j')
			}
			qui matrix `covariates'[`t',1]=`r(estimate)'
			qui matrix `covariates'[`t',2]=`r(se)'
			qui matrix `covariates'[`t',3]=`r(z)'
			qui matrix `covariates'[`t',4]=`r(p)'
			qui matrix `covariates'[`t',5]=`r(lb)'
			qui matrix `covariates'[`t',6]=`r(ub)'
			local ++t
		}
		else{
			if "`categorical'"!=""{
				local first=0 
				foreach l in `levelsof`cat1'' {
					if `first'==0 {
						local ++first
					}
					else{
						if `first'==1 {
							qui lincom [T`k']`l'.`cat1'-[T`j']`l'.`cat1'
							qui matrix `covariates'[`t',1]=`r(estimate)'
							qui matrix `covariates'[`t',2]=`r(se)'
							qui matrix `covariates'[`t',3]=`r(z)'
							qui matrix `covariates'[`t',4]=`r(p)'
							qui matrix `covariates'[`t',5]=`r(lb)'
							qui matrix `covariates'[`t',6]=`r(ub)'
							local ++t
							local ++first
						}
					}
				}
			}
			else{
				qui lincom [T`k']`cont1'-[T`j']`cont1'
				qui matrix `covariates'[`t',1]=`r(estimate)'
				qui matrix `covariates'[`t',2]=`r(se)'
				qui matrix `covariates'[`t',3]=`r(z)'
				qui matrix `covariates'[`t',4]=`r(p)'
				qui matrix `covariates'[`t',5]=`r(lb)'
				qui matrix `covariates'[`t',6]=`r(ub)'
				local ++t
			}		
		}
	}
}	
		
forvalues j=1/`timemax'{
	if "`continuous'"!=""|"`categorical'"!="" {
		qui lincom _b[/var(e.T`j')]
	}
	else {
		qui lincom _b[/var(T`j')]
	}
	qui matrix `covariates'[`t',1]=`r(estimate)'
	qui matrix `covariates'[`t',2]=`r(se)'
	qui matrix `covariates'[`t',3]=`r(z)'
	qui matrix `covariates'[`t',4]=`r(p)'
	qui matrix `covariates'[`t',5]=`r(lb)'
	qui matrix `covariates'[`t',6]=`r(ub)'
	local ++t
}
forvalues j=1/`=`timemax'-1'{
	if "`continuous'"!=""|"`categorical'"!="" {
		forvalues k=`=`j'+1'/`timemax'{
			qui lincom _b[/cov(e.T`j',e.T`k')]
			qui matrix `covariates'[`t',1]=`r(estimate)'
			qui matrix `covariates'[`t',2]=`r(se)'
			qui matrix `covariates'[`t',3]=`r(z)'
			qui matrix `covariates'[`t',4]=`r(p)'
			qui matrix `covariates'[`t',5]=`r(lb)'
			qui matrix `covariates'[`t',6]=`r(ub)'
			local ++t
		}
	}
	else{
		forvalues k=`=`j'+1'/`timemax'{
			qui lincom _b[/cov(T`j',T`k')]
			qui matrix `covariates'[`t',1]=`r(estimate)'
			qui matrix `covariates'[`t',2]=`r(se)'
			qui matrix `covariates'[`t',3]=`r(z)'
			qui matrix `covariates'[`t',4]=`r(p)'
			qui matrix `covariates'[`t',5]=`r(lb)'
			qui matrix `covariates'[`t',6]=`r(ub)'
			local ++t
		}
	}
}
forvalues i=1/ `nbcont' {
   qui lincom `cont`i''
   qui matrix `covariates'[`t',1]=`r(estimate)'
   qui matrix `covariates'[`t',2]=`r(se)'
   qui matrix `covariates'[`t',3]=`r(z)'
   qui matrix `covariates'[`t',4]=`r(p)'
   qui matrix `covariates'[`t',5]=`r(lb)'
   qui matrix `covariates'[`t',6]=`r(ub)'
   local ++t
}
forvalues i=1/ `nbcat' {
   local first=0
   foreach j in `levelsof`cat`i''' {
	   if `first'==0 {
	      local ++first
	   }
	   else {
		   qui	lincom `j'.`cat`i''
		   qui matrix `covariates'[`t',1]=`r(estimate)'
		   qui matrix `covariates'[`t',2]=`r(se)'
		   qui matrix `covariates'[`t',3]=`r(z)'
		   qui matrix `covariates'[`t',4]=`r(p)'
		   qui matrix `covariates'[`t',5]=`r(lb)'
		   qui matrix `covariates'[`t',6]=`r(ub)'
		   local ++t
	   }
   }
}
*matrix list `covariates'


/*************************************************************************************************************
OUTPUTS
*************************************************************************************************************/

if "`postpce'"=="" {
	local t=1
	local diffname
	*set trace on
	di "{hline 73}"
	di  as text _col(60) "<--95% IC -->"
	di  "Items" _col(60) "Lower" _col(68) "Upper"
	di "parameters" _col(13) "category" _col(25) "Estimate" _col(37) "s.e." _col(48) "z" _col(56) "p" _col(59) " Bound" _col(68) "Bound"
	di "{hline 73}"
	*set trace on
	forvalues i=1/`nbitems' {
	   *local l=1
	   forvalues j=1/`modamax`i'' {
		  if "`rsm'"==""|`j'==1 {
			  if `j'==1 {
				 di as text abbrev("``i''",19) _c
			  }
			  di as text _col(20) %5.2f "`j'" as result _col(28) %5.2f `diff'[`t',1]  _col(36) %5.2f `diff'[`t',2] _col(44) %5.2f `diff'[`t',3] _col(52) %5.2f `diff'[`t',4] _col(60) %5.2f `diff'[`t',5] _col(68) %5.2f `diff'[`t',6] 
			  local ++t
			  *local ++l
			  local diffname `diffname' `j'.``i''
		  }
	   }
	}
	if "`rsm'"!="" {
		forvalues k=2/`modamax' {
			di as text "tau`k'"  as result _col(28) %5.2f `diff'[`t',1]  _col(36) %5.2f `diff'[`t',2] _col(44) %5.2f `diff'[`t',3] _col(52) %5.2f `diff'[`t',4] _col(60) %5.2f `diff'[`t',5] _col(68) %5.2f `diff'[`t',6] 
			local diffname `diffname' tau`k'
			local ++t
		}
	}
	di as text "{hline 73}"
	local t=1
	local listmoy
	local listvar
	local listcov
	forvalues j=1/`timemax'{
		local listvar `listvar' Variance_T`j'
		forvalues k=`=`j'+1'/`timemax'{
			local listcov `listcov' Cov_T`j'_T`k'
		}
		forvalues k=`=`j'+1'/`timemax'{
			local listmoy `listmoy' Mean_diff_T`j'_T`k'
		}
	}
	local n: word count `listmoy' `listvar' `listcov' `continuous' 
	forvalues i=1/`n' {
		local v: word `i' of `listmoy' `listvar' `listcov' `continuous' 
		di as text _col(1) %5.2f "`v'" as result _col(28) %5.2f `covariates'[`t',1]  _col(36) %5.2f `covariates'[`t',2] _col(44) %5.2f `covariates'[`t',3] _col(52) %5.2f `covariates'[`t',4] _col(60) %5.2f `covariates'[`t',5] _col(68) %5.2f `covariates'[`t',6] 
		local ++t
	}

	local rn Variance `continuous'

	local n: word count of  `categorical' 
	local catname
	forvalues i=1/`n' {
		local v: word `i' of `categorical'
		local first=1
		local saute=1
		foreach j in `levelsof`cat`i''' {
			if `saute'==0 {
				if `first'==1 {
					di as text _col(1) abbrev("`v'",19) _c
				}
				di  as text _col(20) %5.2f "`j'" as result _col(28) %5.2f `covariates'[`t',1]  _col(36) %5.2f `covariates'[`t',2] _col(44) %5.2f `covariates'[`t',3] _col(52) %5.2f `covariates'[`t',4] _col(60) %5.2f `covariates'[`t',5] _col(68) %5.2f `covariates'[`t',6] 
				local ++first
				local rn `rn' `j'.`n'
				local ++t
				local catname `catname' `j'.`v'
			}
			else {
			   local saute=0
			}
		}
		*local ++t
	}
	di as text "{hline 73}"
	if "`visit'"==""{
		di
		qui su `latent'
		*qui local PSI=1-(`r(sd)')^2/((`covariates'[1,1])+(`r(sd)')^2)
		*di as text "Variance of the estimated latent variable: " as result %4.2f `=(`r(sd)')^2'
        tempvar se2latent
		qui gen `se2latent'=(`selatent')^2
		qui su `se2latent'
		di as text "Mean squared std error of the latent variable: " as result %4.2f `=(`r(mean)')'
		di as text "Global variance of the latent variable: " as result %4.2f `=((`covariates'[1,1])+(`r(mean)'))'
		local PSI=(`covariates'[1,1])/((`covariates'[1,1])+(`r(mean)'))
		di as text "PSI: " as result %4.2f `PSI' _c
		if "`continuous'"!=""|"`categorical'"!="" {
		    di as text " (without adjustment on covariates)"
	    }
		else {
		    di
		}
		di
	}



	matrix colnames `covariates'=`cn'
	matrix rownames `covariates'=`rn'
}

/*************************************************************************************************************
FIT TESTS
*************************************************************************************************************/
if "`visit'"==""{
	tempname fit
	qui matrix `fit'=J(`nbitems',4,.)
	matrix colnames `fit'=OUTFIT INFIT "Standardized OUTFIT" "Standardized INFIT"
	matrix rownames `fit'=`varlist'
	*matrix list `fit'

	tempvar Tcum TInf cum
	qui gen `Tcum'=0
	qui gen `TInf'=0
	if "`postpce'"=="" {
		di as text "{hline 90}"
		di as text _col(60) "<---  Standardized   --->"
		di as text "Items" _col(34) "OUTFIT" _col(50) "INFIT" _col(64) "OUTFIT" _col(80) "INFIT"
		di as text "{hline 90}"
		di as text "Referenced values*" _col(29) "[" %4.2f `=1-6/sqrt(`nbobs')' ";" %4.2f `=1+6/sqrt(`nbobs')' "]" _col(44) "[" %4.2f `=1-2/sqrt(`nbobs')' ";" %4.2f `=1+2/sqrt(`nbobs')' "]" _col(60) "[-2.6;2.6]" _col(75) "[-2.6;2.6]"
		di as text "Referenced values**" _col(29) "[0.75;1.30]" _col(44) "[0.75;1.30]" _col(60) "[-2.6;2.6]" _col(75) "[-2.6;2.6]"
		di as text "{hline 90}"
	}
	*set trace on
	local chi2=0
	local chi2_old=0
	forvalues g=1/`nbgroups' {
	   local chi2_g`g'=0
	   local chi2_old_g`g'=0
	}
	forvalues i=1/`nbitems' {
		if "`rsm'"=="" {
			local mm=`modamax`i''
		}
		else {
			local mm `modamax'
		}	
		tempvar cum_old``i'' c_old0_``i'' Inf_old``i'' y_old``i'' y2_old``i''  	
		tempvar cum``i'' c0_``i'' Inf``i'' C``i'' C2``i'' C3``i'' y``i'' y2``i'' z``i'' z2``i'' i``i'' 

		local d=1
		local d_old=1
		qui gen `cum``i'''=0
		qui gen `cum_old``i'''=0
		forvalues k=1/`mm' {
			local d `d'+exp(`k'*`latent2'-`diffmat2'[`i',`k'])
			local d_old `d_old'+exp(`k'*`latent'-`diffmat2'[`i',`k'])
		}
		qui gen `c0_``i'''=1/(`d')
		qui gen `c_old0_``i'''=1/(`d_old')
		forvalues k=1/`mm' {
			tempvar c`k'_``i'' c_old`k'_``i''
			qui gen `c`k'_``i'''=exp(`k'*`latent2'-`diffmat2'[`i',`k'])/(`d')
			qui gen `c_old`k'_``i'''=exp(`k'*`latent'-`diffmat2'[`i',`k'])/(`d')
			qui replace `cum``i'''=`cum``i'''+`c`k'_``i'''*`k'
			qui replace `cum_old``i'''=`cum_old``i'''+`c_old`k'_``i'''*`k'
		}
		qui gen `Inf``i'''=0
		qui gen `Inf_old``i'''=0
		qui gen `C``i'''=0
		forvalues k=0/`mm' {
			qui replace `Inf``i'''=`Inf``i'''+(`k'-`cum``i''')^2*`c`k'_``i'''
			qui replace `Inf_old``i'''=`Inf_old``i'''+(`k'-`cum_old``i''')^2*`c_old`k'_``i'''
			qui replace `C``i'''=`C``i'''+(`k'-`cum``i''')^4*`c`k'_``i'''
		}
		qui count if ``i''!=.
		local n``i''=r(N)
		
		qui gen `C2``i'''=`C``i'''/((`Inf``i''')^2)
		qui su `C2``i'''
		local q2o``i''=(`r(mean)'-1)/((`n``i'''))

		qui gen `C3``i'''=`C``i'''-(`Inf``i''')^2
		qui su `C3``i'''
		local n=r(sum)
		qui su `Inf``i'''
		local d=r(sum)
		local q2i``i''=`n'/((`d')^2)
		
		//di "``i'' qo = `=sqrt(`q2o``i''')' qi = `=sqrt(`q2i``i''')'"
		
		qui replace `Tcum'=`Tcum'+`cum``i'''
		qui replace `TInf'=`TInf'+`Inf``i'''
		qui gen `y``i'''=``i''-`cum``i'''
		qui gen `y_old``i'''=``i''-`cum_old``i'''
		qui gen `y2``i'''=(`y``i''')^2
		qui gen `y2_old``i'''=(`y_old``i''')^2
		qui gen `z``i'''=(`y``i'''/sqrt(`Inf``i'''))
		local chi2_``i''=0
		local chi2_old_``i''=0
		forvalues g=1/`nbgroups' {
			qui su `y2``i''' if `group'==`g'
			local n=r(sum)
			qui su ``i'' if `group'==`g'
			local n1=r(sum)
			qui su `cum``i''' if `group'==`g'
			local n2=r(sum)
			qui su `Inf``i''' if `group'==`g'
			local d=r(sum)
			*qui count if `group'==`g'
			*local eff=r(N)
			*di "chi2_`g'_``i''=`chi2'+/*`eff'**/(`n1'-`n2')^2/(`d')"
			local chi2=`chi2'+/*`eff'**/(`n1'-`n2')^2/(`d')
			local chi2_``i''=`chi2_``i'''+/*`eff'**/(`n1'-`n2')^2/(`d')
			local chi2_g`g'=`chi2_g`g''+/*`eff'**/(`n1'-`n2')^2/(`d')
			qui su `y2_old``i''' if `group'==`g'
			local n_old=r(sum)
			qui su ``i'' if `group'==`g'
			local n1_old=r(sum)
			qui su `cum_old``i''' if `group'==`g'
			local n2_old=r(sum)
			qui su `Inf_old``i''' if `group'==`g'
			local d_old=r(sum)
			local chi2_old=`chi2_old'+(`n1_old'-`n2_old')^2/(`d_old')
			local chi2_old_``i''=`chi2_old_``i'''+(`n_old')/(`d_old')
			local chi2_old_g`g'=`chi2_old_g`g''+(`n_old')/(`d_old')
		}
		di "Item ``i'' Chi2``i''=`chi2_``i''' et chi2=`chi2'  Chi2_old=`chi2_old_``i''' et chi2_old=`chi2_old' "
		*su `z``i'''
		label variable `z``i''' "Standardized residuals associated to ``i''"
		label variable `latent' "Latent trait"
		*set trace on
		if  "`graphs'"!=""&"`graphitems'"=="" {
			if "`filesave'"!="" {
				local fs saving("`dirsave'//residuals_``i''",replace)
			}
			qui graph twoway scatter `z``i''' `latent', name(residuals``i'',replace) title("Standardized residuals associated to ``i''") `fs'
		}
		*set trace off
		qui gen `z2``i'''=(`z``i''')^2
		qui su `z2``i'''
		local OUTFIT``i''=`r(mean)'
		qui matrix `fit'[`i',1]=`OUTFIT``i'''
		local OUTFITs``i''=((`r(mean)')^(1/3)-1)*(3/sqrt(`q2o``i'''))+sqrt(`q2o``i''')/3
		qui matrix `fit'[`i',3]=`OUTFITs``i'''
		qui su `Inf``i''' if ``i''!=.
		local sumw``i''=r(sum)
		qui gen `i``i'''=`Inf``i'''*`z2``i''' 
		qui su `i``i''' if ``i''!=.
		local INFIT``i'' = `=`r(sum)'/`sumw``i''''
		qui matrix `fit'[`i',2]=`INFIT``i'''
		local INFITs``i''=(`=`r(sum)'/`sumw``i''''^(1/3)-1)*(3/sqrt(`q2i``i'''))+sqrt(`q2i``i''')/3
		qui matrix `fit'[`i',4]=`INFITs``i'''
		if "`postpce'"=="" {
			di abbrev("``i''",19) _col(35) %5.3f `OUTFIT``i''' _col(50) %5.3f `INFIT``i''' _col(64) %6.3f `OUTFITs``i''' _col(79) %6.3f `INFITs``i''' 
		}
	}
	if "`postpce'"=="" {
		di as text "{hline 90}"
		di as text "*: As suggested by Wright (Smith, 1998)
		di as text "**: As suggested by Bond and Fox (2007)
	}
	if "`geninf'"!="" {
	    gen `geninf'=`TInf'
	}
}
*set trace off
/*************************************************************************************************************
ESTIMATION OF THE WEIGHTED ML ESTIMATORS
**************************************************************************************************************/
*set trace on
*di "estimation `wmliterate'"
if "`postpce'"!="" {
    local conv=10
	local it=`wmliterate'
	di "Iteration : `it'"
	while(`conv'>=1) {
		di "Itération `it' : conv=`conv'"
		tempvar   sinf
		qui gen `sinf'=sqrt(TInf_`=`it'-1')
		`qui' pcm `varlist' [iweight=`sinf'],diff(loulou) wmliterate(`it') geninf(TInf_`it') genlt(lt_`it') 
		tempvar  ecart_`it' 
		qui gen `ecart_`it''=abs(lt_`it'-lt_`=`it'-1')
		qui su `ecart_`it''
		local conv =r(mean)
		local ++it
	}
    exit	
}


/*************************************************************************************************************
Categories/Items/Test Characteristics Curves and Information graphs
*************************************************************************************************************/
if "`visit'"==""{
	if "`graphs'"!="" {

		tempfile savefile
		qui save `savefile'

		*qui clear
		qui drop _all
		
		qui set obs 2000
		qui gen u=(_n-1000)/250
		qui gen Tcum=0
		qui gen TInf=0
		forvalues i=1/`nbitems' {
		   local scatteri`i' 
		   forvalues g=1/`nbgroups' {
			  local x=`groups'[`g',`=`nbitems'+1']
			  local y=`groups'[`g',`i']
			  local s1=`groups'[`g',`=`nbitems'+2']
			  local seuil=30
			  local s vtiny
			  *set trace on
			  foreach lab in  /*tiny*/ vsmall small medsmall medium medlarge  large   vlarge  huge vhuge /*ehuge*/ {
				  if `s1'>`seuil' {
					 local s `lab'
				  }
				  local seuil=`seuil'+10
			  }
			  local scatteri`i' `scatteri`i'' || scatteri `y' `x' , mcolor(black) msize(`s') legend(off)
			*set trace off
			}
			local d=1
			qui gen cum``i''=0
			*set trace on
			if "`rsm'"=="" {
				local mm=`modamax`i''
			}
			else {
				local mm `modamax'
			}	
			forvalues k=1/`mm' {
				local d `d'+exp(`k'*u-`diffmat2'[`i',`k'])
			}
			qui gen c0_``i''=1/(`d')
			label variable c0_``i'' "Pr(``i''=0)"
			forvalues k=1/`mm' {
				qui gen c`k'_``i''=exp(`k'*u-`diffmat2'[`i',`k'])/(`d')
				qui replace cum``i''=cum``i''+c`k'_``i''*`k'
				label variable c`k'_``i'' "Pr(``i''=`k')"
			}
			qui gen Inf``i''=0
			forvalues k=1/`mm' {
				qui replace Inf``i''=Inf``i''+(`k'-cum``i'')^2*c`k'_``i''
			}
			if "`graphitems'"=="" {
				if "`filesave'"!="" {
					local fsc saving("`dirsave'//CCC_``i''",replace)
					local fsi saving("`dirsave'//ICC_``i''",replace)
				}
				qui graph twoway line c*_``i'' u , name(CCC``i'', replace) title(Categories Characteristic Curve (CCC) of ``i'') ytitle("Probability") xtitle("Latent trait") `fsc'
				qui graph twoway line cum``i'' u, name(ICC``i'',replace) title("Item Characteristic Curve (ICC) of ``i''") ytitle("Score to the item") xtitle("Latent trait") `scatteri`i'' `fsi'
			}
			qui replace Tcum=Tcum+cum``i''
			qui replace TInf=TInf+Inf``i''	
			label variable Inf``i'' "``i''"
		}
		if "`filesave'"!="" {
			local fst saving("`dirsave'//TCC",replace)
			local fsteo saving("`dirsave'//TCCeo",replace)
			local fsi saving("`dirsave'//ICC",replace)
			local fsti saving("`dirsave'//TIC",replace)
			local fsm saving("`dirsave'//map",replace)
		}
		qui graph twoway line Tcum u, name(TCC,replace) title("Test Characteristic Curve (TCC)") ytitle("Score to the test") xtitle("Latent trait") `fst'
		qui graph twoway line Inf* u, name(IIC,replace) title("Item Information Curves") ytitle("Information") xtitle("Latent trait")  `fsi'
		qui graph twoway line TInf u, name(TIC,replace) title("Test Information Curve") ytitle("Information") xtitle("Latent trait")  `fsti'

		local scatteri
		forvalues g=1/`nbgroups' {
			local x=`groups'[`g',`=`nbitems'+1']
			local y=`groups'[`g',`=`nbitems'+3']
			  local s1=`groups'[`g',`=`nbitems'+2']
			  local seuil=30
			  local s vtiny
			  *set trace on
			  foreach lab in  tiny vsmall small medsmall medium medlarge  large   vlarge  huge vhuge /*ehuge*/ {
				  if `s1'>`seuil' {
					 local s `lab'
				  }
				  local seuil=`seuil'+10
			  }
			local scatteri `scatteri' || scatteri `y' `x' , mcolor(black) msize(`s') legend(off)
		}
		qui graph twoway line Tcum u , name(TCCeo,replace) title("Test Characteristic Curve (TCC)") ytitle("Score to the test") xtitle("Latent trait") `scatteri' `fsteo'
	
	
	
	
	
/*************************************************************************************************************
MAP
*************************************************************************************************************/
	*set trace on
		gen eff=0
		local effmax=0
		forvalues g=1/`nbgroups' {
			local eff=`groups'[`g',`=`nbitems'+2']
			if `groups'[`g',`=`nbitems'+2']>`effmax' {
			   local effmax=`groups'[`g',`=`nbitems'+2']
			}
			local lat=round(`groups'[`g',`=`nbitems'+1'],0.004)
			*di "replace eff=`eff' if round(u,0.004)==`lat'"
			qui replace eff=`eff' if round(u,0.004)==`lat'
		}
		gen density=normalden(u)*sqrt(`covariates'[1,1])
		label variable eff "Frequencies"
		label variable u "Latent trait"
		label variable TInf "Information curve"
		label variable density "Density function of the latent trait"
		local scatteri
		qui su u if eff!=0
		*set trace on
		local floor=floor(`r(min)')
		local ceil=ceil(`r(max)')
		forvalues i=1/`nbitems' {
		   local color`i':word `i' of `color'
		   local y=-`i'*`effmax'/`nbitems'
		   forvalues l=1/`modamax' {
			   local x=`diffmat'[`i',`l']
			   local scatteri `scatteri' || scatteri `y' `x' "`l'" ,mcolor(black) mlabcolor(black)
			   if `x'<`floor' {
				  local floor=floor(`x')
			   }
			   if `x'>`ceil'&`x'!=. {
				  local ceil=ceil(`x')
			   }
		   }
		   local scatteri `scatteri' || scatteri `y' `floor' "``i''",mcolor(black) mlabcolor(black) msize(vtiny)
		}
		qui su eff
		local maxe=ceil(`=(floor(`r(max)'/10)+1)*10')
		qui su TInf
		local maxi=ceil(`r(max)')
		qui su density
		local maxd=ceil(`r(max)')
		qui drop if u<`floor'|u>`ceil'
		qui graph twoway (bar eff u, barwidth(.2) yaxis(1) legend(off) xlabel(`floor'(1)`ceil')) (line TInf u,yaxis(2)) (line density u,yaxis(3))    `scatteri'  , name(map,replace) ytitle("Frequencies")  ylabel(-`maxi'(`=`maxi'/5')`maxi' ,axis(2)) ylabel(-`maxd'(`=`maxd'/5')`maxd' ,axis(3)) ylabel(-`maxe'(`=`maxe'/5')`maxe' ,axis(1)) title("Individuals/items representations (Map)") xsize(12) ysize(9) note("Red line: Information curve - Green line : Density of the latent trait") xtitle("Latent trait") `fsm'


		qui use `savefile', clear

	}
}



/*************************************************************************************************************
RESULTS BY GROUP
*************************************************************************************************************/
if "`visit'"==""{
	*set trace on
    tempname matscorelt
	qui matrix `matscorelt'=J(`=`nbitems'*`modamax'+1',3,.)
	di
	di as text "{hline 57}"
	di _col(31) "Latent Variable" _col(50) "Expected"
	di "Group" _col(10) "Score" _col(20) "Freq" _col(32) "Mean" _col(42) "s.e." _col(53) "Score"
	di as text "{hline 57}"
	forvalues g=1/`nbgroups' {
		qui count if `group'`multivisit'==`g' 
		local eff`g'=r(N)
		qui count if `score'`multivisit'!=.&`group'`multivisit'==`g'
		local n=r(N)
		di as text "`g' (n=" as result `eff`g'' as text ")" _c
		if `n'>0 {
			qui su `score'`multivisit' if `group'`multivisit'==`g'
			local scoremin`g'=`r(min)'
			local scoremax`g'=`r(max)'
			forvalues s=`scoremin`g''/`scoremax`g'' {
				qui count if `group'`multivisit'==`g'&`score'`multivisit'==`s'
				local eff=r(N)
				qui su `latent' if `group'`multivisit'==`g'&`score'`multivisit'==`s'
				local mean=r(mean)
				qui su `selatent' if `group'`multivisit'==`g'&`score'`multivisit'==`s'
				local se=r(mean)
				qui su `Tcum' if `group'`multivisit'==`g'&`score'`multivisit'==`s'
				local exp=r(mean)
				if `eff'>0 {
				   di as text _col(10) %5.0f `s' as result _col(20) %4.0f `eff' _col(30) %6.3f `mean' _col(40) %6.3f `se' _col(53) %5.2f `exp'
				}
				qui matrix `matscorelt'[`=`s'+1',1]=`s'
				qui matrix `matscorelt'[`=`s'+1',2]=`mean'
				qui matrix `matscorelt'[`=`s'+1',3]=`se'
			}
		}
		qui count if `group'`multivisit'==`g'&`score'`multivisit'==.
		local eff=r(N)
		qui su `latent' if `group'`multivisit'==`g'&`score'`multivisit'==.
		local mean=r(mean)
		qui su `selatent' if `group'`multivisit'==`g'&`score'`multivisit'==.
		local se=r(mean)
		qui su `Tcum' if `group'`multivisit'==`g'&`score'`multivisit'==.
		local exp=r(mean)
		if `eff'>0 {
			di  as text _col(10) "    ." as result _col(20) %4.0f `eff' _col(30) %6.3f `mean' _col(40) %6.3f `se' _col(53) %5.2f `exp'
		}
		di as text "{hline 57}"
	}
	*matrix list `matscorelt'
}	

/*************************************************************************************************************
EQUATING
*************************************************************************************************************/
*set trace on
if "`eqset1'"!="" {

    tokenize `eqset1'
	local nbset1: word count `eqset1'
	forvalues i=1/`nbset1' {
	   local eq1_`i':word `i' of `eqset1'
	   *di "set1 (`nbset1') : `eq1_`i''"
	}
	*di "`eqset1'"
    tokenize `eqset2'
	local nbset2: word count  `eqset2'
	*di "`eqset2'"
	forvalues i=1/`nbset2' {
	   local eq2_`i':word `i' of `eqset2'
	   *di "set2  (`nbset2') : `eq2_`i''"
	}	



	tempfile fileeq
	qui save `fileeq',replace
	forvalues t=1/2{
		local scoremaxset`t'=0
		forvalues i=1/`nbset`t'' {
		    *di "`eq`t'_`i''"
			*local tmp=""
			local scoremaxset`t'=`scoremaxset`t''+`modamax`eq`t'_`i'''
		}
	}
	drop _all
	*set trace on
	qui set obs `=(`scoremaxset1'+`scoremaxset2'+2)*3'
	forvalues t=1/2 {
		qui gen scoreset`t'=.
		qui gen scoreset`t'm=.
		qui gen scoreset`t'p=.	
	}
	forvalues i=0/`scoremaxset1' {
	    qui replace scoreset1=`i' in `=`i'+1'
	    qui replace scoreset1m=`i' in `=`i'+1+(`scoremaxset1'+`scoremaxset2'+2)'
	    qui replace scoreset1p=`i' in `=`i'+1+(`scoremaxset1'+`scoremaxset2'+2)*2'
	}
	forvalues i=`=`scoremaxset1'+2'/`=`scoremaxset1'+`scoremaxset2'+2' {
	    qui replace scoreset2=`i'-`scoremaxset1'-2 in `i'
	    qui replace scoreset2m=`i'-`scoremaxset1'-2 in `=`i'+(`scoremaxset1'+`scoremaxset2'+2)'
	    qui replace scoreset2p=`i'-`scoremaxset1'-2 in `=`i'+(`scoremaxset1'+`scoremaxset2'+2)*2'
	}
	local s=0
	local eqset1b
	foreach i in `eqset1' {
		qui gen s1_`i'=0
		forvalues m=1/`modamax`i'' {
		   qui gen s1_`i'_`m'=0 in 1/`=`scoremaxset1'+1'
		   qui replace s1_`i'_`m'=1 if scoreset1>`s' in 1/`=`scoremaxset1'+1'
		   qui replace s1_`i'=s1_`i'+s1_`i'_`m'
		   local ++s
		}
	    local eqset1b `eqset1b' s1_`i'
	}
	local s=0
	local eqset2b
	foreach i in `eqset2' {
		qui gen s2_`i'=0
		forvalues m=1/`modamax`i'' {
		   qui gen s2_`i'_`m'=0 in `=`scoremaxset1'+2'/`=`scoremaxset1'+`scoremaxset2'+2'
		   qui replace s2_`i'_`m'=1 if scoreset2>`s' in `=`scoremaxset1'+2'/`=`scoremaxset1'+`scoremaxset2'+2'
		   qui replace s2_`i'=s2_`i'+s2_`i'_`m'
		   local ++s
		}
		local eqset2b `eqset2b' s2_`i'
	}
	tokenize `varlist'
	tempname diffset1 diffset2
	*matrix list `diffmat'
	forvalues t=1/2 {
		qui matrix `diffset`t''=J(`nbset`t'',`modamax',.)
		local n=1
		local listset`t'
		foreach j in `eqset`t'' {
			forvalues i=1/`nbitems' {
				if "`j'"=="``i''" {
					local listset`t' `listset`t'' `i'
					forvalues m=1/`modamax' {
						qui matrix `diffset`t''[`n',`m']=`diffmat'[`i',`m']
					}
					local ++n						
				}
			}
		}
	}
	*matrix list `diffset1'
	*matrix list `diffset2'
	local var=`covariates'[1,1]
	qui gen lt=.
	*qui gen selt=.
	forvalues t=1/2 {
		tempname matscorelt`t'
		qui pcm `eqset`t'b', diff(`diffset`t'') var(`var') minsize(1)
		qui matrix `matscorelt`t''=r(matscorelt)
		*di "matscorelt`t':"
		*matrix list `matscorelt`t''
		forvalues i=0/`scoremaxset`t'' {
			qui replace lt=`matscorelt`t''[`=`i'+1',2] if scoreset`t'==`i'
			qui replace lt=`matscorelt`t''[`=`i'+1',2]+1.96*`matscorelt`t''[`=`i'+1',3] if scoreset`t'p==`i'
			qui replace lt=`matscorelt`t''[`=`i'+1',2]-1.96*`matscorelt`t''[`=`i'+1',3] if scoreset`t'm==`i'
			*qui replace selt=`matscorelt`t''[`=`i'+1',3] if scoreset`t'==`i'
		}
		qui ipolate scoreset`t' lt, gen(score`t') epolate
	}
	qui ipolate scoreset1 lt, gen(score1bis) epolate


	*list
	forvalues t=1/2 {
		qui replace score`t'=scoreset`t'm if scoreset`t'm!=.
		qui replace score`t'=scoreset`t'p if scoreset`t'p!=.
		qui replace score1=score1bis if score1==.
		qui replace score`t'=0 if score`t'<0
		qui replace score`t'=`scoremaxset`t'' if score`t'>`scoremaxset`t''
	}
	forvalues t=1/2 {
		tempname matscore`t'
		qui matrix `matscore`t''=J(`=`scoremaxset`t''+1',7,.)
		forvalues s=0/`scoremaxset`t'' {
		    qui matrix `matscore`t''[`=`s'+1',1]=`s'
			qui su lt if scoreset`t'==`s'
		    qui matrix `matscore`t''[`=`s'+1',2]=r(mean)
			qui su lt if scoreset`t'm==`s'
		    qui matrix `matscore`t''[`=`s'+1',3]=r(mean)
			qui su lt if scoreset`t'p==`s'
		    qui matrix `matscore`t''[`=`s'+1',4]=r(mean)
			qui su score`=3-`t'' if scoreset`t'==`s'
		    qui matrix `matscore`t''[`=`s'+1',5]=r(mean)
			qui su score`=3-`t'' if scoreset`t'm==`s'
		    qui matrix `matscore`t''[`=`s'+1',6]=r(mean)
			qui su score`=3-`t'' if scoreset`t'p==`s'
		    qui matrix `matscore`t''[`=`s'+1',7]=r(mean)
		}
		matrix colnames `matscore`t'' =score`t' lt lt- lt+ score`=3-`t'' score`=3-`t''- score`=3-`t''+
		*matrix list `matscore`t''
		di
		di "{hline 78}"
		di "EQUATING SET`t' TO SET`=3-`t''"
		di "{hline 78}"
		di "Set`t' : `eqset`t''"
		di "Set`=3-`t'' : `eqset`=3-`t'''"
		di "{hline 78}"
		di            _col(20) "<----- Latent trait ----->" _col(52) "<------- Score `=3-`t'' --------->"
		di "Score`t'" _col(20) "Estimated" _col(39) "[95%IC]" _col(52) "Estimated" _col(72) "[95%IC]"
		di "{hline 78}"
		forvalues s=0/`scoremaxset`t'' {
			di %4.0f `matscore`t''[`=`s'+1',1] _col(24) %5.2f `matscore`t''[`=`s'+1',2] _col(33) "[" %5.2f `matscore`t''[`=`s'+1',3] ";" %5.2f `matscore`t''[`=`s'+1',4] "]" _col(56) %5.2f `matscore`t''[`=`s'+1',5] _col(66) "[" %5.2f `matscore`t''[`=`s'+1',6] ";" %5.2f `matscore`t''[`=`s'+1',7] "]"
		}
		di "{hline 78}"
		return matrix score`t'_to_`=3-`t''=`matscore`t''
		if "`eqgraph'"!="" {
			*twoway (line lt scoreset1) (line lt scoreset2), name(eq1)
			*twoway (line scoreset1 scoreset2 lt) , name(eq2)
			*twoway (line score1 scoreset1m scoreset1p score2), name(eq3)
			*twoway (line score2 scoreset2m scoreset2p score1), name(eq4)
			twoway (line score`t' score`=3-`t'' if scoreset`=3-`t''!=.) (line score`t' score`=3-`t'' if scoreset`=3-`t''m!=.) (line score`t' score`=3-`t'' if scoreset`=3-`t''p!=.), title("Equating score of the Set `t' from the Set `=3-`t''") ytitle("Score `t'") xtitle("Score `=3-`t''") ylabel(0(1)`scoremaxset`t'') xlabel(0(1)`scoremaxset`=3-`t''') name(eq`t'to`=3-`t'')
			*twoway (line score2 score1 if scoreset1!=.) (line score2 score1 if scoreset1m!=.) (line score2 score1 if scoreset1p!=.), name(eq6)
		}
	}
	*save prout, replace
	*clear
	
	qui use `fileeq',clear
}


/*************************************************************************************************************
RETOUR AU FICHIER INITIAL ET SAUVEGARDE DES NOUVELLES VARIABLES
*************************************************************************************************************/
if "`visit'"!="" {
	tempfile sauv
	qui keep `latent'* `selatent'* `id' `visit'
    qui reshape wide , i(`id') j(`visit')
	qui sort `id'	
	qui save `sauv', replace
	restore,preserve
	*su
	if "`replace'"!=""&("`genlt'"!=""|"`genfit'"!="") {
    	capture drop `genlt' se`genlt' `geninf'
    }
	tempname idorder
	qui gen `idorder'=_n
	qui sort `id'
	qui merge 1:1 `id' using `sauv'
	qui sort `idorder'
    qui drop `idorder'	
}
else {
	restore,not
}

/*************************************************************************************************************
CREATION DU DOCX
*************************************************************************************************************/

if "`docx'"!="" {
    putdocx clear
    putdocx begin
	putdocx paragraph
	putdocx text ("General informations") , bold underline font(,14) smallcaps
	putdocx paragraph
	putdocx text ("Number of individuals: `nbobs'")
	putdocx paragraph
	putdocx text ("Number of complete individuals: `nbobsssmd'")
	putdocx paragraph
    putdocx text ("Number of items: `nbitems'")
	putdocx paragraph
    putdocx text ("List of items: `varlist'")
	putdocx paragraph
    putdocx text ("Date: $S_DATE, $S_TIME")
	putdocx paragraph
	local model Partial Credit Model (PCM)
	if "`rsm'"!="" {
	   local model Rating Scale Model (RSM)
	}
	putdocx text ("Model: `model'")
	putdocx paragraph
	putdocx text ("Marginal log-likelihood: `ll'")

	putdocx paragraph
	putdocx text ("Estimation of the parameters") , bold underline font(,14) smallcaps
	putdocx table tablename = matrix(`diff') , nformat(%9.3f) rownames colnames border(start, nil) border(insideH, nil) border(insideV, nil) border(end, nil)   
	qui putdocx table tablename = matrix(`covariates') , nformat(%9.3f) rownames colnames border(start, nil) border(insideH, nil) border(insideV, nil) border(end, nil)   
	putdocx paragraph
	putdocx text ("Fit indexes for items") , bold underline font(,14) smallcaps
	qui putdocx table tablename = matrix(`fit') , nformat(%9.3f) rownames colnames border(start, nil) border(insideH, nil) border(insideV, nil) border(end, nil)   
	local extension png
}

/*************************************************************************************************************
SAUVEGARDE DES GRAPHIQUES
*************************************************************************************************************/

*set trace on
if "`filesave'"!="" {
    if "`graphs'"!="" {
		if "`docx'"!="" {
		    putdocx pagebreak
		    putdocx paragraph
	        putdocx text ("General graphs") , bold underline font(,14) smallcaps
		}
		foreach i in TCC TCCeo TIC IIC map {
			if "`extension'"!="" {
			    qui graph export "`dirsave'//`i'.`extension'", replace name(`i')     
            }
			*graph display `i' 
			*qui graph save "`dirsave'//`i'", replace    
			if "`docx'"!="" {
			    putdocx paragraph
			    putdocx image "`dirsave'//`i'.png", height(10cm)
			}
		}
	    *discard
		if "`graphitems'"=="" {
			forvalues i=1/`nbitems' {
				if "`docx'"!="" {
				    putdocx paragraph
				    putdocx text ("Graphs for ``i''") , bold underline font(,14) smallcaps
				}
				foreach j in CCC ICC residuals {
					*graph display `j'``i'' 
					*qui graph save "`dirsave'//`j'_``i''", replace    

	 			    if "`extension'"!="" {
					    qui graph export "`dirsave'//`j'_``i''.`extension'", replace name(`j'``i'') 
					}
					if "`docx'"!="" {
						putdocx paragraph
						putdocx image "`dirsave'//`j'_``i''.png" , height(10cm)
					}
				}
			}
		}
	}
}
if "`docx'"!="" {
    putdocx save "`dirsave'//`docx'.docx", replace
}


/*************************************************************************************************************
RETURNS
*************************************************************************************************************/


matrix colnames `diff'=Estimate "s.e." z p lb ul
matrix colnames `covariates'=Estimate "s.e." z p lb ul
matrix rownames `diff'=`diffname'
matrix rownames `covariates'=Variance `continuous' `catname'
return matrix difficulties=`diff'
return matrix covariates=`covariates'
return matrix matscorelt=`matscorelt'
capture restore, not
end