*! version 1 : February 15th, 2012
*! Myriam Blanchin
************************************************************************************************************
* raschlong: Estimation of the parameters of a model of the Rasch family in a longitudinal setting
*
* Version 1 : February 15th, 2012: pcm et rasch, contraintes sur la matrice de variance-covariance et difficultés d'items fixées 
*
************************************************************************************************************/

program define raschlong,eclass
syntax varlist [, NBT(int 2) DIFFiculties(string) VAR(string) ]

*the dataset should be in wide format
*the matrix of difficulties should contain `nbit' rows and `nbmodpos' columns
*the item parameters are assumed to be constant with time
*same number of modalities for each item
*varlist should contain for each timepoint the list of variables containing the answers to the items - example: T1item1 T1item2 T1item3 T2item1 T2item2 T2item3
preserve

tempname diff varcov item id temps one obs estbeta estvar 
/*di "item" "`item'"
di "temps" "`temps'"
di "obs" "`obs'"
di "diff" "`diff'"
di "varcov" "`varcov'"
di "id" "`id'"
di "one" "`one'"*/
*verif varlist
tokenize `varlist'
local nbittot:word count `varlist'
local nbit=`nbittot'/`nbt'
if `=int(`nbit')'!=`nbit'{
	di in red "The number of variables should be equal to (the number of items * the number of timepoints). Please correct it."
	error 198
	exit
}


/*verif nb diff items*/
if "`difficulties'"!=""{
	matrix `diff'=`difficulties'
	if `nbit'!=`=rowsof(`diff')'{
		di in red "The number of rows of the matrix of item parameters should be equal to the number of items. Please correct it."
		error 198
		exit
	 }
	local nbmodat=colsof(`diff')+1
}
else{
	local modcount=word("`varlist'",1)
	qui tab `modcount'
    local nbmodat=r(r)
}

*verif matcov
if "`var'"!=""{
	matrix `varcov'=`var'
	if `=rowsof(`var')'!=`nbt'{
		di in red "The covariance matrix is incorrectly specified. Please correct it."
		error 198
		exit
	}
}

qui{
	*reshape of the dataset
	gen `one'=1
	collapse(sum) wt2=`one', by (`varlist')
	*list in 1/3
	forvalues t=1/`nbt'{
		forvalues i=1/`nbit'{
			local vartemp=word("`varlist'",`=(`t'-1)*`nbit'+`i'')
			gen tps`t'item`i'=`vartemp'
		}
	}

	gen `id'=_n
	local list=""
	forvalues i=1/`nbit'{
		local list="`list' tps@item`i'"
	}
	reshape long `list', i(`id') j(`temps')
	*list in 1/20
	reshape long tpsitem, i(`id' `temps') j(`item')
	*list in 1/20
	drop if tpsitem==.
	gen `obs'=_n
	expand `nbmodat'
	sort `id' `temps' `item' `obs'
	*list in 1/20
	tab `temps', gen(t)
	forvalues t=1/`nbt'{
		by `obs', sort: gen x`t'=(_n-1)*t`t'
	}
	gen chosen=.
	forvalues t=1/`nbt'{
		replace chosen=tpsitem==x`t' if `temps'==`t'
	}
	tab `item', gen(it)

	if "`difficulties'"==""{
		forvalues i=1/`nbit'{
			forvalues g=1/`=`nbmodat'-1'{
				gen d`i'_`g'=.
			}
		}
		forvalues t=1/`nbt'{
			forvalues i=1/`nbit'{
				forvalues g=1/`=`nbmodat'-1'{
					replace d`i'_`g'=-1*it`i'*(x`t'>=`g') if `temps'==`t'
				}
			}
		}
	}
	else{
		gen offset=0
		forvalues t=1/`nbt'{
			forvalues i=1/`nbit'{
				local sumdiff=0 
				forvalues g=1/`=`nbmodat'-1'{
					local sumdiff=`sumdiff'-`diff'[`i',`g']
					replace offset=`sumdiff' if it`i'==1 & x`t'==`g'
				}
			}			
		}
	}
	*list in 1/25
	local listeq ""
	forvalues t=1/`nbt'{
		eq slope`t':x`t'
		local listeq "`listeq' slope`t'"
	}
	if "`var'"!=""{
		matrix C=cholesky(`varcov')
		constraint define 1 [__01_1]x1=`=C[1,1]'
		constraint define 2 [__01_2]x2=`=C[2,2]'
		constraint define 3 [__01_2_1]_cons=`=C[2,1]'
	}
}
tempname b V
if "`difficulties'"==""{
	if "`var'"==""{
		gllamm x1-x`nbt' d1_1-d`nbit'_`=`nbmodat'-1',i(`id') eqs(`listeq') link(mlogit) expand(`obs' chosen o) weight(wt) nocons nrf(`nbt') adapt
	}
	else{
		gllamm x1-x`nbt' d1_1-d`nbit'_`=`nbmodat'-1',i(`id') eqs(`listeq') link(mlogit) constraints(1 2 3)  expand(`obs' chosen o) weight(wt) nocons nrf(`nbt') adapt
	}
	matrix `b'=e(b)
	matrix EstDiff=J(`nbit',`=`nbmodat'-1',.)
	forvalues j=1/`nbit'{
		matrix EstDiff[`j',1]=`b'[1,`=`nbt'+(`j'-1)*(`nbmodat'-1)'..`=`nbt'-1+`j'*(`nbmodat'-1)']
	}

}
else{
	if "`var'"==""{
		gllamm x1-x`nbt',i(`id') eqs(`listeq') link(mlogit) expand(`obs' chosen o) weight(wt) nocons nrf(`nbt') adapt offset(offset)
	}
	else{
		gllamm x1-x`nbt',i(`id') eqs(`listeq') link(mlogit) constraints(1 2 3) expand(`obs' chosen o) weight(wt) nocons nrf(`nbt') adapt offset(offset)
	}
	matrix `b'=e(b)
}
matrix `V'=e(V)
matrix EstMu=`b'[1,1..`=`nbt'-1']
matrix Var=e(chol)*e(chol)'
matrix se2=vecdiag(`V'[1..`=`nbt'-1',1..`=`nbt'-1'])

/*if test
test x2=x3=0
 r(chi2)
  r(df)
  r(p)
*/
di
di
di in gr "{hline 91}"
di in gr _col(39) "RASCH FAMILY MODEL"
di in gr "{hline 91}"
di in gr "Number of timepoints:   " in ye `nbt'
di in gr "Number of items:   " in ye `nbit'
di in gr "Number of modalities per item:   " in ye `nbmodat'
di in gr "{hline 91}"

local collist ""
forvalues i=1/`=`nbmodat'-1'{
	local collist "`collist' dj_`i'"
}
local rowlist ""
forvalues i=1/`nbit'{
	local rowlist "`rowlist' item`i'"
}
	
if "`difficulties'"==""{
	di in gr "Item parameters estimations:"
	matrix colnames EstDiff = `collist'
	matrix rownames EstDiff = `rowlist'
	matrix list EstDiff,noheader format(%7.3f)
	ereturn matrix Diff=EstDiff
}
else{
	di in gr "Item parameters fixed to:"
	matrix colnames `diff' = `collist'
	matrix rownames `diff' = `rowlist'
	matrix list `diff',noheader format(%7.3f)
	ereturn matrix Diff=`diff'
}
di
di
if "`var'"==""{
	di in gr "Covariance matrix estimations:"
}
else{
	di in gr "Covariance matrix fixed to:"
}
matrix list Var,noheader nonames format(%7.3f)
di
di
di in gr "Time effect estimations:"
di in gr _col(19) "Coef" _col(29) "S.E."
forvalues t=1/`=`nbt'-1'{
	di in gr "Time" `=`t'+1' in ye _col(16) %7.2f EstMu[1,`t'] _col(26) %7.2f `=sqrt(se2[1,`t'])'
}

ereturn matrix Var=Var



end