*! Version 1.2 June 3, 2013
************************************************************************************************************
* Stata program : raschres
* Analysis of the residuals of the Rasch model
* Version 1 : June 16, 2010
* Version 1.1 : July 15, 2011
* Version 1.2 : June 3, 2013 /*formatting outputs*/
*
*
* Jean-benoit Hardouin, phD, Assistant Professor
* EA4275 - SPHERE
* Team of Biostatistics, Pharmacoepidemiology and Subjective Measures in Health Sciences
* University of Nantes - Faculty of Pharmaceutical Sciences
* France
* jean-benoit.hardouin@univ-nantes.fr
*
* News about this program : http://www.anaqol.org
*
* Copyright 2010-2013 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 raschres,rclass
syntax varlist [if] [in] , diff(string) theta(varname) [genres(string) genfit(string)  vargroup(varname) PCAres DETails]
tempfile saveraschres
if "`if'"!=""|"`in'"!="" {
   *keep `if' `in'
}
*qui save `saveraschres',replace

tempname res

qui count
local ntotal=r(N)
tokenize `varlist'
local nbitems: word count `varlist'
qui tempname score  lt
qui genscore `varlist',score(`score')
qui gen `lt'=`theta'
forvalues i=1/`nbitems' {
   local z2`i'=0
}

forvalues i=1/`nbitems'{
   qui tempname p`i' res`i'
   qui gen `p`i''=exp(`lt'-`diff'[1,`i'])/(1+exp(`lt'-`diff'[1,`i']))
   qui gen `res'`i'=(1-`p`i'')/sqrt(`p`i''*(1-`p`i'')) if ``i''==1
   qui replace `res'`i'=(0-`p`i'')/sqrt(`p`i''*(1-`p`i'')) if ``i''==0
}
su
local chi2=0
forvalues i=1/`nbitems' {
    qui tempname res2
    qui gen `res2'`i'=`res'`i'^2
    qui su `res2'`i'
    local chi2=`chi2'+`r(sum)'
    label variable `res'`i' "``i''"
}
su
tempname cor  Ev
matrix `cor'=J(`nbitems',`nbitems',.)
local list
forvalues i=1/`nbitems' {
   forvalues j=1/`nbitems' {
      qui corr `res'`i' `res'`j'
      matrix `cor'[`i',`j']=r(rho)
   }
   local list "`list' ``i''"
}
matrix colname `cor'=`list'
matrix rowname `cor'=`list'
qui pca `res'*
matrix `Ev'=e(Ev)

tempname lt2 ltcl group  df nbmiss rep total
qui egen `nbmiss'=rowmiss(`varlist')
qui gen `rep'=`nbitems'-`nbmiss'
qui gen `lt2'=`lt' if `score'!=0&`score'!=`rep'
if "`vargroup'"=="" {
   qui gengroup `lt2', newvar(`group') cont      det  minsize(33)
}
else {
   qui `group'=`vargroup'
}
qui su `group'
local maxgroup=r(max)

if "`details'"!="" {
    di "{hline 42}"
    di _col(12) "Min" _col(29) "Max"
    di "{dup 42:-}"
    di "" _col(8) "" _col(14) "Latent" _col(28) "" _col(31) "Latent"  _col(53)  ""
    di "Group" _col(8) "Score" _col(15) "trait" _col(25) "Score" _col(32) "trait"  _col(42)  "n"
    di "{hline 42}"
    forvalues i=1/`maxgroup' {
       qui count if `group'==`i'
       local n`i'=r(N)
       qui su `score' if `group'==`i'
       local scoremin`i'=r(min)
       local scoremax`i'=r(max)
       qui su `lt' if `group'==`i'
       local ltmin`i'=r(min)
       local ltmax`i'=r(max)
       di "`i'" _col(9) %4.0f `scoremin`i'' _col(14)  %6.3f `ltmin`i'' _col(26) %4.0f `scoremax`i'' _col(31) %6.3f `ltmax`i''  _col(39) %4.0f  `n`i''
    }
    di "{hline 42}"
}

local nj=0
forvalues i=1/`nbitems' {
   qui count if ``i''!=.
   local nj=`nj'+r(N)
}
qui gen `df'=(`nj'-(`ntotal' + `nbitems' - 1 ))/`nj'*`rep'
qui gen `total' =`score'
local mult=(`nj'-(`ntotal' + `nbitems' - 1 ))/`nj'


forvalues i=1/`nbitems' {
   qui count if (`total'!=0&`total'!=`rep')&``i''!=.
   local mult`i'=r(N)
   local df`i'=`mult'*`mult`i''
   tempname res2`i'
   qui gen `res2`i''=(`res'`i')^2
   qui su `res2`i''   if `total'!=0&`total'!=`rep'
   local sumsq`i'=r(sum)
   local fit`i'=(`sumsq`i''-`df`i'')/sqrt(2*`df`i'')
}

qui gen `ltcl'=.
forvalues g=1/`maxgroup' {
   qui su `lt' if `group'==`g'
   local ltgroup`g'=r(mean)
   local N`g'=r(N)
   qui replace `ltcl'=`ltgroup`g'' if `group'==`g'
}
forvalues i=1/`nbitems' {
   local chi`i'=0
   forvalues g=1/`maxgroup' {
      qui su ``i'' if `group'==`g'
      local po`i'_`g'=r(mean)
      local pe`i'_`g'=exp(`ltgroup`g''-`diff'[1,`i'])/(1+exp(`ltgroup`g''-`diff'[1,`i']))
      local std`i'_`g'=2*`N`g''*(`po`i'_`g''-`pe`i'_`g'')^2/sqrt(.5*(`po`i'_`g''+`pe`i'_`g'')*(1-.5*(`po`i'_`g''+`pe`i'_`g'')))
      local chi`i'=`chi`i''+`std`i'_`g''
   }
   di "chi item ``i'' = `chi`i''"
}

local chisq=0
tempname y2 vy2
qui gen `y2'=0
qui gen `vy2'=0
forvalues i=1/`nbitems' {
   qui replace `y2'=`y2'+`res'`i'^2  if ``i''!=.
   qui replace `vy2'=`vy2'+2*tanh((`lt2'-`diff'[1,`i'])/2)*sinh(`lt2'-`diff'[1,`i']) if ``i''!=.
   tempname y2`i' vy2`i'
   qui gen `y2`i''=`res'`i'^2
   qui su `y2`i'' if `total'!=0&`total'!=`rep'&``i''!=.
   local sumy2`i'=r(sum)
   qui tab `ltcl'
   qui gen `vy2`i''=2*tanh((`lt2'-`diff'[1,`i'])/2)*sinh(`lt2'-`diff'[1,`i'])
   di "item ``i''"
    su `vy2`i'' if `total'!=0&`total'!=`rep'&``i''!=.
   local sumvy2`i'=r(sum)
   local fit`i'=(`sumy2`i''-`df`i'')/sqrt(`sumvy2`i'')
   local fit2`i'=`df`i''*(ln(`sumy2`i'')-ln(`df`i''))/sqrt(`sumvy2`i'')
   local chisq`i'=0
   forvalues g=1/`maxgroup' {
      qui su `y2`i'' if `group'==`g'&``i''!=.
      local sumy2`i'_g`g'=r(sum)
      local n`i'g`g'=r(N)
      qui su `vy2`i'' if `group'==`g'&``i''!=.
      local sumvy2`i'_g`g'=r(sum)
      local fit`i'_g`g'=(`sumy2`i'_g`g''-`n`i'g`g'')/sqrt(`sumvy2`i'_g`g'')
      local chisq`i'=`chisq`i''+(`fit`i'_g`g'')^2
   }
   local chisq=`chisq'+`chisq`i''
   /*fit2 est l'équivalent de FitResid et chisq est correct*/
}
local RMSEA=sqrt(`chisq'/((`nbitems'*(`maxgroup'-1)-1)*(`ntotal'-1)))

di in gr "{hline 55}"
di in gr "Items" _col(18) "FitResid" _col(31) "ChiSq" _col(42) "df" _col(55) "p"
di in gr "{hline 55}"
forvalues i=1/`nbitems' {
   di in ye abbrev("``i''",18) _col(19) %7.3f `fit2`i'' _col(30) %6.2f `chisq`i'' _col(40) %4.0f  `=`maxgroup'-1' _col(50) %6.4f 1-chi2(`=`maxgroup'-1',`chisq`i'')
}
di in gr "{dup 55:-}"
di in ye "Total"  _col(30) %6.2f `chisq' _col(40) %4.0f  `=(`maxgroup'-1)*`nbitems'' _col(50) %6.4f 1-chi2(`=(`maxgroup'-1)*`nbitems'',`chisq`i'')
di in ye "RMSEA"  _col(30) %6.4f `RMSEA'
di in gr "{hline 55}"
di

if "`pcares'"!="" {
    di "{hline 36}"
    di "Correlation matrix between residuals"
    di "{hline 36}"
    matrix list `cor',noheader format(%5.3f)
    di
    di "{hline 28}"
    di "Eigenvalues of the residuals"
    di "{hline 28}"
    matrix rowname `Ev'=Eigenvalues
    tempname Evp E
    matrix `Evp'=`Ev'/`nbitems'
    matrix `E'=`Ev' \ `Evp'
    matrix rowname `E'=Eigenvalues Proportion
matrix list `E',noheader format(%5.3f)
}

*use `saveraschres', clear
local dfind=`df'/`ntotal'*`nbitems'
if "`genfit'"!="" {
   qui gen `genfit'=`df'*(ln(`y2')-log(`df'))/sqrt(`vy2') `if' `in'
}
if "`genres'"!="" {
   forvalues i=1/`nbitems' {
      qui gen `genres'`i'=`res'`i'  `if' `in'
   }
}

end