*! version 5 : October 22th, 2013
*! Jean-Benoit Hardouin, Myriam Blanchin
************************************************************************************************************
* raschpower: Estimation of the power of the Wald test in order to compare the means of the latent trait in two groups of individuals
*
* Version 1 : January 25, 2010 (Jean-Benoit Hardouin)
* Version 1.1 : January 26, 2010 (Jean-Benoit Hardouin)
* Version 1.2 : November 1st, 2010 (Jean-Benoit Hardouin)
* version 1.3 : May 2th, 2011 (Jean-Benoit Hardouin)
* version 1.4 : July 7th, 2011 (Jean-Benoit Hardouin) : minor corrections
* version 1.5 : July 11th, 2011 (Jean-Benoit Hardouin) : minor corrections
* version 2 : August 30th, 2011 (Jean-Benoit Hardouin, Myriam Blanchin) : corrections
* version 3 : October 18th, 2011 (Jean-Benoit Hardouin, Myriam Blanchin) : Extension to the PCM, -method- option, -nbpatterns- options, changes in the presentation of the results
* version 3.1 : October 25th, 2011 (Jean-Benoit Hardouin, Myriam Blanchin) : POPULATION+GH method
* version 3.2 : February 6th, 2012 (Jean-Benoit Hardouin, Myriam Blanchin) : minor corrections
* version 4 : January 30th, 2013 (Jean-Benoit Hardouin, Myriam Blanchin) : Extension to longitudinal design
* version 4.1 : June 3rd, 2013 (Jean-Benoit Hardouin, Myriam Blanchin) : detail option
* version 5 : October 22th, 2013 (Jean-Benoit Hardouin, Myriam Blanchin) : Extension to longitudinal design with polytomous items
*
* Jean-benoit Hardouin, jean-benoit.hardouin@univ-nantes.fr
* Myriam Blanchin, myriam.blanchin@univ-nantes.fr
* EA 4275 "Biostatistics, Pharmacoepidemiology and Subjectives Measures in Health"
* Faculty of Pharmaceutical Sciences - University of Nantes - France
* http://www.sphere-nantes.org
*
* News about this program : http://www.anaqol.org
* FreeIRT Project : http://www.freeirt.org
*
* Copyright 2010-2013 Jean-Benoit Hardouin, Myriam Blanchin
*
* 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 raschpower5,rclass
syntax [varlist] [, n0(int 100) n1(int 100) Gamma(real .5) Difficulties(string) Var(string) Method(string) NBPatterns(int 2) nodata EXPectedpower(real -1) LONGitudinal freevar DETail]
version 11

tempfile raschpowerfile
capture qui save "`raschpowerfile'",replace
tempname db d dlong matvar
if "`difficulties'"=="" {
   matrix `d'=[-1\-0.5\0\0.5\1]
}
else {
   matrix `d'=`difficulties'
}
local nbitems=rowsof(`d')
local nbmodat=colsof(`d')+1

if "`longitudinal'"==""{
	local nbt=1
	local nbtotitems=`nbitems'
	local nbpatmax=2*`nbmodat'^`nbitems'
	if "`var'"==""{
		local var=1
	}
	else{	
		capture confirm number `var'
		if !_rc {
			local var=`var'
		}
		else{ 
			matrix `matvar'=`var'
			capture confirm matrix `matvar'
			if !_rc & colsof(`matvar')==1 & rowsof(`matvar')==1{
				local var=`matvar'[1,1]
			}
			else{
				di in red "{p}The {hi:var} option should contain a number or a 1x1 matrix for transversal studies{p_end}"
				error 198
				exit
			}
		}
	}
	matrix `dlong'=`d'


	local sumd=0
	forvalues numit=1/`nbitems'{
		local sumd=`sumd'+`d'[`numit',1]
	}
	local sumd=`sumd'/`nbitems'
	if `=abs(`sumd')'>=`=2*sqrt(`var')'{
		di in red "{p}WARNING: It is not recommended to use Raschpower when the gap between the global mean of the latent trait (fixed to 0) and the mean of the item parameters is greater than or equal to 2 standard deviation of the latent trait. {p_end}"
	}
	
}
else{
	local nbt=2
	local nbtotitems=2*`nbitems'
	local nbpatmax=`nbmodat'^(`nbitems'*2)
	local n1=0
	local mean1=0 
	local mean2=`mean1'+`gamma' 
	if "`var'"==""{
		matrix `matvar'=(1,0\0,1)
	}
	else{
		matrix `matvar'=`var'
		capture confirm matrix `matvar'
		if _rc | colsof(`matvar')!=2 | rowsof(`matvar')!=2{
			di in red "{p}The {hi:var} option should contain a 2x2 matrix for longitudinal studies{p_end}"
			error 198
			exit
		}
	}
	matrix `dlong'=J(`nbtotitems',`=`nbmodat'-1',.)
	matrix `dlong'[1,1]=`d'
	matrix `dlong'[`=`nbitems'+1',1]=`d'
}

if substr("`method'",1,3)=="pop" | substr("`method'",1,3)=="POP"{
	local method POPULATION+GH
}

if "`method'"=="MEAN+GH"&`nbpatterns'*(`n1'+`n0')>=`nbpatmax' {
   di in gr "The MEAN+GH will be inefficient compared to GH since the maximal number of pattern's responses"
   di in gr "is lesser than the number of pattern retained by the MEAN+GH method."
   di in gr "The -method- option is replaced by GH."
   local method GH
}
else if ("`method'"=="MEAN+GH" | "`method'"=="MEAN") & `nbpatmax'>1000000{
	di in red "The number of patterns is too high for the chosen method"
	exit
}
else if "`method'"=="" {
    if `nbpatmax'<2000 {
        local method GH
    }
    else  {
        local method POPULATION+GH
    }
}

if "`method'"!="MEAN+GH" & "`method'"!="MEAN" & "`method'"!="GH" & "`method'"!="POPULATION+GH"{
	di in red "Invalid method name"
	exit
}


if "`longitudinal'"==""{
	di in gr "Number of individuals in the first group:  " in ye `n0'
	di in gr "Number of individuals in the second group: " in ye `n1'
	di in green "Group effect: " in ye `gamma'
	di in  gr "Variance of the latent trait: " in ye `var'
}
else{
	di in gr "Number of individuals at each time:  " in ye `n0'
	di in green "Time effect: " in ye `gamma'
	di in  gr "Variance matrix of the latent trait: " 
	matrix list `matvar',noheader
}
di in gr "Number of items: " in ye `nbitems'
di in green "Difficulties parameters of the items: "
tempname dd
matrix `dd'=`d''
local items
forvalues i=1/`nbitems' {
   local items "`items' item`i'"
}
local modalities
forvalues i=1/`=`nbmodat'-1' {
   local modalities "`modalities' delta_`i'"
}

matrix colnames `dd'=`items'
matrix rownames `dd'=`modalities'
matrix list `dd',noblank nohalf noheader
if "`detail'"!=""{
	di in gr "Method:  " in ye "`method'"
	di in gr "Number of studied response's patterns: " in ye `nbpatmax'
}
matrix `dd'=`d'
local gamma=`gamma'

local tmp=1
qui matrix `db'=J(`=`nbitems'*(`nbmodat'-1)',1,.)
forvalues j=1/`nbitems' {
	forvalues m=1/`=`nbmodat'-1' {
		qui matrix `db'[`tmp',1]=`d'[`j',`m']
		local ++tmp
	}
}

if "`data'"=="" {
	clear
	if "`method'"!="POPULATION+GH"{
		local temp=`nbmodat'^(`nbtotitems')
		qui range x 0 `=`temp'-1' `temp'
		qui g t=x
		loc i=`nbtotitems'
		qui count if t>0
		loc z=r(N)
		qui while `z'>0 {
			qui gen item`=`nbtotitems'-`i'+1'=floor(t/`nbmodat'^`=`i'-1')
			qui replace t=mod(t,`nbmodat'^`=`i'-1')
			qui count if t>0
			loc z=r(N)
			loc i=`i'-1
		}
		drop t
		if "`longitudinal'"==""{
			qui expand 2
			qui gen group=0 in 1/`temp'
			qui replace group=1 in `=`temp'+1'/`=2*`temp''
		}	
	}
	else {
		if "`longitudinal'"==""{
			qui simirt, clear pcm(`d') cov(`var') group(`=`n1'/(`n1'+`n0')') deltagroup(`gamma') nbobs(1000000)
			qui drop lt1
			qui contract item* group, freq(freq)
			qui count
			local patobs=r(N)
			if `patobs'>=`=(`n0'+`n1')*`nbpatterns''{
				qui gen keep=0
				qui gsort +group -freq
				qui replace keep=1 in 1/`=`nbpatterns'*`n0''
				qui gsort -group -freq
				qui replace keep=1 in 1/`=`nbpatterns'*`n1''
			}
			else{
				qui gen keep=1
			}
		}
		else{
			qui simirt, clear pcm(`dlong') covmat(`matvar') mu(`mean1' `mean2') dim(`nbitems' `nbitems') nbobs(1000000)
			forvalues j=1/`nbitems'{
				rename itemA`j' item`j'
				rename itemB`j' item`=`nbitems'+`j''
			}
			qui drop lt1 lt2
			qui contract item*, freq(freq)
			local patobs=r(N)
			if `patobs'>=`=`n0'*`nbpatterns''{
				qui gen keep=0
				qui gsort -freq
				qui replace keep=1 in 1/`=`nbpatterns'*`n0''
			}
			else{
			  qui gen keep=1
			}
		}
		qui keep if keep==1
		qui count
		local retain=r(N)
		di "Number of kept patterns:`retain'"
		local method GH
	}
	if "`longitudinal'"==""{
		qui gen mean1=-`n1'*`gamma'/(`n0'+`n1') if group==0
		qui replace mean1=`n0'*`gamma'/(`n0'+`n1') if group==1
	}
	else{
		qui gen mean1=`mean1'
		qui gen mean2=`mean2'
	}
	if "`method'"=="GH" {
		local temp=`nbmodat'^(`nbtotitems')
		local diff0=0
		qui gen proba=.
		local dixj=10
		qui count
		local tmp=r(N)

		forvalues i=1/`tmp' {
			local dix=floor(`tmp'/10)
			if mod(`i',`dix')==0 {
				if "`dixj'"!="10" {
					di ".." _c
				}
				di "`dixj'%" _c
				local dixj=`dixj'+10
			}
			forvalues t=1/`nbt'{
				local int`t'=1			 
				forvalues j=`=(`t'-1)*`nbitems'+1'/`=`t'*`nbitems'' {
					qui su item`j' in `i'
					local rep=r(mean)
					local diff0=0
					local diff1=`dlong'[`j',1]
					local sum "1+exp(x`t'-`diff1')"
					forvalues m=2/`=`nbmodat'-1' {
					   local diff`m'=`diff`=`m'-1''+`dlong'[`j',`m']
					   local sum "`sum'+exp(`m'*x`t'-`diff`m'')"
					}
					local int`t' "(`int`t''*exp(`rep'*x`t'-`diff`rep''))/(`sum')"
				}
			}	
			if "`longitudinal'"==""{
				qui su mean1 in `i'
				local mean=r(mean)
				qui gausshermite `int1',mu(`mean') sigma(`=sqrt(`var')') display name(x1)
				qui replace proba=r(int) in `i'
			}
			else{
				local int "`int1'*`int2'"
				matrix mean=(`mean1',`mean2')
				qui gausshermite `int',mu(mean) var(`matvar') display name(x1 x2)
				qui replace proba=r(int) in `i'
			}
		}
		di
	}
	else {
		qui gen proba=1
		forvalues t=1/`nbt'{
			forvalues i=`=(`t'-1)*`nbitems'+1'/`=`t'*`nbitems'' {
				local diff0=0
				local diff1=`dlong'[`i',1]
				qui gen eps0=1
				qui gen eps1=exp(mean`t'-`diff1')
				qui gen d=eps0+eps1
				forvalues m=2/`=`nbmodat'-1' {
					local diff`m'=`diff`=`m'-1''+`dlong'[`i',`m']
					qui gen eps`m'=exp(`m'*mean`t'-`diff`m'')
					qui replace d=d+eps`m'
				}
				local listeps
				forvalues m=0/`=`nbmodat'-1' {
					qui replace proba=proba*eps`m'/d if item`i'==`m'
					local listeps `listeps' eps`m'
				}
				qui drop `listeps' d
			}
		}
		if "`method'"=="MEAN+GH" {
			set tracedepth 1
			if "`longitudinal'"==""{
				qui gen keep=0
				qui gsort -group -proba
				local min=min(`=`nbmodat'^`nbitems'',`=`n1'*`nbpatterns'')
				qui replace keep=1 in 1/`min'
				qui gsort +group -proba
				local min=min(`=`nbmodat'^`nbitems'',`=`n0'*`nbpatterns'')
				qui replace keep=1 in 1/`min'
				qui keep if keep==1
				qui su proba if group==0
				local sumproba0=r(sum)*100
				qui su proba if group==1
				local sumproba1=r(sum)*100
			}
			else{
				qui gen keep=0
				qui gsort -proba
				local min=min(`nbpatmax',`=`n0'*`nbpatterns'')
				qui replace keep=1 in 1/`min'
				qui keep if keep==1
			}

			qui drop keep proba
			local diff0=0
			qui gen proba=.
			qui count
			local nnew=r(N)
			di in gr "Number of studied response's patterns for the GH step: " in ye `nnew'
			if "`longitudinal'"==""{
				di in gr "(" in ye %6.2f `sumproba0' in gr "% of the group 0 and " in ye %6.2f `sumproba1' in gr "% of the group 1)"
			}
			local dixj=10
			forvalues i=1/`nnew' {
				local dix=floor(`nnew'/10)
				if mod(`i',`dix')==0 {
					if "`dixj'"!="10" {
						di ".." _c
					}
					di "`dixj'%" _c
					local dixj=`dixj'+10
				}
				forvalues t=1/`nbt'{
					local int`t'=1
					forvalues j=`=(`t'-1)*`nbitems'+1'/`=`t'*`nbitems'' {
						qui su item`j' in `i'
						local rep=r(mean)
						local diff0=0
						local diff1=`dlong'[`j',1]
						local sum "1+exp(x`t'-`diff1')"
						forvalues m=2/`=`nbmodat'-1' {
							local diff`m'=`diff`=`m'-1''+`dlong'[`j',`m']
							local sum "`sum'+exp(`m'*x`t'-`diff`m'')"
						}
						local int`t' "(`int`t''*exp(`rep'*x`t'-`diff`rep''))/(`sum')"
					}
				}
				if "`longitudinal'"==""{
					qui su mean1 in `i'
					local mean=r(mean)
					qui gausshermite `int1',mu(`mean') sigma(`=sqrt(`var')') display name(x1)
					qui replace proba=r(int) in `i'
				}
				else{
					local int "`int1'*`int2'"
					matrix mean=(`mean1',`mean2')
					qui gausshermite `int',mu(mean) var(`matvar') display
					qui replace proba=r(int) in `i'
				}
			}
		}
	}
	qui gen eff=proba
	qui gen eff2=.  
   
	if "`longitudinal'"==""{
		qui keep item* eff* group proba
		local p1=1/`n1'
		local p0=1/`n0'
		qui replace eff2=floor(eff/`p1') if group==1
		qui replace eff2=floor(eff/`p0') if group==0
		qui replace eff=eff-eff2*(`p1'*group+`p0'*(1-group))
		qui su eff2 if group==1
		local aff1=r(sum)
		qui su eff2 if group==0
		local aff0=r(sum)

		local unaff1=`n1'-`aff1'
		local unaff0=`n0'-`aff0'
		qui gen efftmp=eff2
		qui gsort + group - eff
		qui replace eff2=eff2+1 in 1/`unaff0'
		qui gsort - group - eff
		qui replace eff2=eff2+1 in 1/`unaff1'

		qui drop if eff2==0
		gsort group item*
		qui expand eff2
	}
	else{
		qui keep item* eff* proba
		local p0=1/`n0'
		qui replace eff2=floor(eff/`p0')
		qui replace eff=eff-eff2*`p0'
		qui su eff2
		local aff0=r(sum)

		local unaff0=`n0'-`aff0'
		qui gen efftmp=eff2
		qui gsort - eff
		qui replace eff2=eff2+1 in 1/`unaff0'

		qui drop if eff2==0
		gsort item*
		qui expand eff2
	}
   qui drop proba eff eff2
}

qui alpha item*
local alpha=r(alpha)

if "`longitudinal'"==""{
	qui gen groupc=-`n1'/(`n0'+`n1') if group==0
	qui replace groupc=`n0'/(`n0'+`n1') if group==1
	if `nbmodat'==2 {
		qui gen i=_n

		tempname diff
		matrix `diff'=`dd''

		qui reshape long item, i(i)
		qui rename item rep
		qui rename _j item

		qui gen  offset=0
		forvalues i=1/`nbitems' {
			qui replace offset=-`diff'[1,`i'] if item==`i'
		}

		constraint 1 _cons=`=ln(`var')'
		qui xtlogit rep groupc ,nocons i(i) offset(offset) constraint(1)
		tempname b V
	}
	else {
		matrix `db'=`db''
		qui pcm item*, fixed(`db') covariates(groupc) fixedmu fixedvar(`var')

	}
	tempname b V
	matrix `b'=e(b)
	matrix `V'=e(V)
	local gammaest=`b'[1,1]
	local se=`V'[1,1]^.5
}
else{
	qui raschlong item*, nbt(`nbt') diff(`d') var(`matvar')
	tempname b V
	matrix `b'=e(b)
	matrix `V'=e(V)
	local gammaest=`b'[1,1]
	local se=`V'[1,1]^.5
	local n1=`n0'
}
local poweruni=1-normal(1.96-`gamma'/`se')


	di
	di
	if "`detail'"!=""{
		di in gr "{hline 91}"
		di _col(60)  "Estimation with the "
		di _col(50)  "Cramer-Rao bound" _col(75) "classical formula"
		di in gr "{hline 91}"
	}

	if "`longitudinal'"==""{
		local power=1-normal(1.96-`gamma'/`se')+normal(-1.96-`gamma'/`se')
		local clpower=normal(sqrt(`n1'*`gamma'^2/((`n1'/`n0'+1)*`var'))-1.96)
		local clnsn=(`n1'/`n0'+1)/((`n1'/`n0')*(`gamma'/sqrt(`var'))^2)*(1.96+invnorm(`poweruni'))^2
		local ratio=(`n0'+`n1')/(`clnsn'*(1+`n1'/`n0'))
		if "`detail'"==""{
			di in gr "{hline 65}"
			di in green "Estimation of the variance of the group effect" _col(56) in ye %10.4f `=`se'^2'
			di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni'
			di in gr "{hline 65}"
		}
		else{
			if "`gammafixed'"==""  {
			   di in green "Estimated value of the group effect" _col(59) in ye  %7.2f `gammaest'
			}
			di in green "Estimation of the s.e. of the group effect" _col(59) in ye %7.2f `se'
			di in green "Estimation of the variance of the group effect" _col(56) in ye %10.4f `=`se'^2'
			di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni' _col(86) in ye %6.4f `clpower'
			di in green "Number of patients for a power of" %6.2f `=`poweruni'*100' "%" _col(59) in ye `n0' "/" `n1' _col(77) in ye %7.2f `clnsn' "/" %7.2f `=`clnsn'*`n1'/`n0''
			di in green "Ratio of the number of patients" in ye %6.2f _col(68)`ratio'
			di in gr "{hline 91}"
		}
	}
	else{
		local clpower=normal(sqrt(`n0'*`gamma'^2)/(2*(`=`matvar'[1,1]'-`=`matvar'[2,1]'))-1.96)
		local clnsn=2*(`=`matvar'[1,1]'-`=`matvar'[2,1]')*(1.96+invnorm(`poweruni'))^2/(`gamma'^2)
		local ratio=(`n0'+`n1')/(`clnsn'*(1+`n1'/`n0'))
		if "`detail'"==""{
			di in gr "{hline 65}"
			di in green "Estimation of the variance of the group effect" _col(56) in ye %10.4f `=`se'^2'
			di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni'
			di in gr "{hline 65}"
		}
		else{
			if "`gammafixed'"==""  {
			   di in green "Estimated value of the time effect" _col(59) in ye  %7.2f `gammaest'
			}
			di in green "Estimation of the s.e. of the time effect" _col(59) in ye %7.2f `se'
			di in green "Estimation of the variance of the time effect" _col(56) in ye %10.4f `=`se'^2'
			di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni' _col(86) in ye %6.4f `clpower'
			di in green "Number of patients for a power of" %6.2f `=`poweruni'*100' "%" _col(63) in ye `n0' _col(85) in ye %7.2f `clnsn'
			di in green "Ratio of the number of patients" in ye %6.2f _col(68)`ratio'
			di in gr "{hline 91}"
		}
	}
	/* ceci est un commentaire */
	if `expectedpower'!=-1 {
	   qui sampsi `=-`gamma'/2' `=`gamma'/2', sd1(`=sqrt(`var')') sd2(`=sqrt(`var')') alpha(0.05) power(`expectedpower') ratio(`=`n1'/`n0'')
	   local expn_1=r(N_1)
	   local expn_2=r(N_2)
	   local expn2=`expn_1'*`ratio'
	   di in green "Number of patients for a power of" %6.2f `=`expectedpower'*100' "%" _col(51) in ye %7.2f `expn2' "/" %7.2f `=`expn2'*`n1'/`n0'' _col(77) in ye %7.2f `expn_1' "/" %7.2f `expn_2'
	}

	return scalar EstGamma=`gammaest'
	return scalar CRbound=`=`se'^2'
	return scalar CRPower=`poweruni'
	return scalar ClPower=`clpower'
	return scalar ClSS=`clnsn'
	return scalar Ratio=`ratio'
	return scalar CronbachAlpha=`alpha'



capture qui use `raschpowerfile',clear

end