************************************************************************************************************
* Stata program : Raschfit
* The Raschfit and the Raschfit-fast procedures to construct sub-scales of items
* Release 2 : June 8, 2004
*
* Historic
* Version 1 (2004-05-06) [Jean-Benoit Hardouin]
*
* Jean-benoit Hardouin, Regional Health Observatory of Orléans - France
* jean-benoit.hardouin@neuf.fr
*
* Use the Stata programs mmsrm, raschtest and gammasym who can be download on http://anaqol.free.fr
* News about this program :http://anaqol.free.fr
*
* FreeIRT Project website : http://freeirt.free.fr
*
* Copyright 2004 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 raschfit,rclass
version 7.0
syntax varlist(min=2 numeric) [,KERnel(integer 0) NBSCales(integer 1) ITEMSorder(string) nofast]
if "`itemsorder'"=="" {
        local itemsorder mspinv
}
local nbitemstot : word count `varlist'
tokenize `varlist'

tempfile raschfitfile
qui save `raschfitfile',replace
preserve

tempname affect
matrix define `affect'=J(1,`nbitemstot',0)
matrix colnames  `affect'=`varlist'

tempname rep item  matbetadim1 matbetadim2
if `kernel'!=0 {
   local listkernel
   forvalues i=1/`kernel' {
             local listkernel `listkernel' `rep'`i'
             matrix `affect'[1,`i']=1
   }
}

local dim=0
local nbitemsnosel=`nbitemstot'
local nbitemstotdim1=`nbitemstot'
local nbitemsnoselukernel=`nbitemstot'-`kernel'


tempvar id betadim1 betadim2
forvalues i=1/`nbitemstot' {
      qui drop if ``i''==.
      rename ``i'' `rep'`i'
}
qui gen `id'=_n
tempfile filescale
qui save `filescale',replace


while `nbitemsnosel'>2&`dim'<`nbscales' {
      use `filescale',replace
      local iteration=0
      local dim=`dim'+1
      if `dim'>1 {
         local kernel=0
         local listkernel
      }

      di in green "SCALE: " in yellow `dim'
      di in green "{hline 9}"
      di
      tempname result`dim'
      local listitemsnosel
      local varlist`dim'
      tokenize `varlist'
      forvalues i=1/`nbitemstot' {
            if `affect'[1,`i']==0 {
                   local varlist`dim' `varlist`dim'' ``i''
                   local listitemsnosel `listitemsnosel' `rep'`i'
            }
      }
      local nbitemsnosel:word count `listitemsnosel'
      if `dim'>1 {
            local nbitemstotdim`dim':word count `listitemsnosel'
      }

      if `kernel'>=2 {
            local fixed=`kernel'
      }
      else {
            local fixed=2
      }
      matrix define `result`dim''=J(`=`nbitemstotdim`dim''-`fixed'',`=`nbitemstotdim`dim''+7',0)

      tempname order`dim' affect`dim'
      matrix `order`dim''=J(1,`nbitemstotdim`dim'',0)
      matrix `affect`dim''=J(1,`nbitemstotdim`dim'',0)

      if "`itemsorder'"=="msp"|"`itemsorder'"=="mspinv" {
            di in green _col(0) "The program is ordering the items"
            qui msp `listkernel' `listitemsnosel',c(-99)  notest kernel(`kernel')
            local scale1 "`r(scale1)'"
            local scalenum1 "`r(scalenum1)'"
            tokenize `scalenum1'
            local listitemsselnum
            forvalues j=`=`nbitemstotdim`dim''+1-`fixed''/`nbitemstotdim`dim'' {
                matrix `order`dim''[1,``j'']=1
                local k:word `j' of `scalenum1'
                matrix `affect`dim''[1,`k']=1
                local listitemsselnum `listitemsselnum' `k'
            }
            forvalues j=1/`nbitemsnosel' {
                matrix `order`dim''[1,`j']=`=`nbitemsnosel'+1-`j''
            }
            tokenize `scale1'
            local listitemssel ``=`nbitemstotdim`dim''-1'' ``nbitemstotdim`dim'''

            local listitemsnosel
            local listitemsnoselnum
            
            if "`itemsorder'"=="mspinv" {
                forvalues j=1/`=`nbitemstotdim`dim''-`fixed'' {
                      local listitemsnosel `listitemsnosel' ``j''
                      local k:word `j' of `scalenum1'
                      local listitemsnoselnum `listitemsnoselnum' `k'
                }
            }
            else if "`itemsorder'"=="msp"{
                 forvalues j=`=`nbitemstotdim`dim''-`fixed''(-1)1 {
                      local listitemsnosel  `listitemsnosel' ``j''
                      local k:word `j' of `scalenum1'
                      local listitemsnoselnum  `listitemsnoselnum' `k'
                 }
            }
      }
      else if "`itemsorder'"=="order" {
           tokenize `listkernel' `varlist`dim''
           local listitemssel
           local listitemsselnum
           local listitemsnosel
           local listitemsnoselnum
           forvalues j=1/`fixed'{
                     local listitemssel `listitemssel' `rep'`j'
                     local listitemsselnum `listitemsselnum' `j'
                     matrix `affect`dim''[1,`j']=1
           }
           forvalues j=`=`fixed'+1'/`nbitemstotdim`dim'' {
                     local listitemsnosel `listitemsnosel' `rep'`j'
                     local listitemsnoselnum `listitemsnoselnum' `j'
           }
      }

      if `dim'>1 {
             tokenize `varlist`dim''
      }
      else {
             tokenize `varlist'
      }

      local nbitemsnosel:word count `listitemsnosel'
      local list
      *tokenize `varlist'
      forvalues i=1/`=`nbitemsnosel'+`fixed'' {
            local tmp:word `i' of "`listitemsnosel' `listitemssel'"
            local list `list' ``i''
      }

      matrix colnames `result`dim''=`list'  iteration nbitems ll1 AIC1 ll2 AIC2 choose

      di _col(0) in green "The kernel of the scale is " in yellow _continue
      forvalues i=1/`fixed' {
            local inum:word `i' of `listitemsselnum'
            di in yellow "``inum'' " _continue
      }
      di
      tokenize `listitemsnosel'
      forvalues i=1/`=`nbitemsnosel-2'' {
            local iteration=`iteration'+1
            di
            di in green _col(10) "iteration:" in yellow" `iteration'"
            di in green _col(10) "{hline 13}"
            di
            qui use `filescale' , clear
            local i2:word `i' of `listitemsnosel'
            local i2num:word `i' of `listitemsnoselnum'
            qui keep `id' `listitemssel' `i2'
            tempname score1 score2
            qui gen `score2'=0
            tokenize `listitemssel'
            local nbitemssel: word count `listitemssel'
            forvalues j=1/`i' {
                      local j2num:word `j' of `listitemsnoselnum'
                      if `affect`dim''[1,`j2num']==1 {
                            matrix `result`dim''[`iteration',`j']=1
                      }
            }
            forvalues j=1/`nbitemssel' {
                      local j2:word `j' of `listitemssel'
                      local j2num:word `j' of `listitemsselnum'
                      qui replace `score2'=`score2'+`j2'
            }
            tokenize `listitemsnosel'
            qui gen `score1'=`score2'+`i2'
            forvalues j=`=`nbitemsnosel'+1'/`=`nbitemsnosel'+`nbitemssel'' {
                      matrix `result`dim''[`iteration',`j']=1
            }
            matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+1']=`iteration'
            matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+2']=`nbitemssel'
            matrix `result`dim''[`iteration',`i']=2

            if "`fast'"=="" {
                        qui count if `score1'!=`=`nbitemssel'+1'&`score1'!=0&`score2'!=`nbitemssel'
                        local N=r(N)
                        di in green _col(10)"The program is estimating the parameters of the models (N=" in yellow "`N'" in green ")"
                        qui raschtest `listitemssel' `i2' if `score1'!=3&`score1'!=0&`score2'!=2, mml notest

                        matrix `matbetadim1'=e(beta),0
                        qui raschtest `listitemssel' if `score1'!=3&`score1'!=0&`score2'!=2, mml notest
                        matrix `matbetadim2'=e(beta),0
                        qui reshape long `rep' , i(`id') j(`item')
                        qui gen `betadim1'=0
                        qui gen `betadim2'=0

                        forvalues j=1/`nbitemssel' {
                                  local j2num:word `j' of `listitemsselnum'
                                  qui replace `betadim1'=-`matbetadim1'[1,`j'] if `item'==`j2num'
                                  qui replace `betadim2'=-`matbetadim2'[1,`j'] if `item'==`j2num'
                        }
                        qui replace `betadim1'= -`matbetadim1'[1,`=`nbitemssel'+1'] if `item'==`i2num'
                        qui count if `rep'==1&`item'==`i2num'&`score1'!=`=`nbitemssel'+1'&`score1'!=0&`score2'!=`nbitemssel'
                        local nb1=r(N)
                        qui count if `item'==`i2num'&`score1'!=`=`nbitemssel'+1'&`score1'!=0&`score2'!=`nbitemssel'
                        local nb2=r(N)

                        qui replace `betadim2'= log(`nb1'/(`nb2'-`nb1')) if `item'==`i2num'
                        qui xi:logit `rep' i.`score1' if `score1'!=`=`nbitemssel'+1'&`score1'!=0&`score2'!=`nbitemssel', nocons offset(`betadim1') iter(100)
                        local N1=e(N)
                        matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+3']=e(ll)
                        matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+4']=2*(-e(ll)+2*(`nbitemssel'+1)-1)

                        tempvar  newscore2
                        qui gen `newscore2'=`score2'
                        qui replace `newscore2'=99 if `item'==`i2num'
                        qui xi:logit `rep' i.`newscore2' if `score1'!=`=`nbitemssel'+1'&`score1'!=0&`score2'!=`nbitemssel',  offset(`betadim2') iter(100)
                        local N2=e(N)
                        local ll=e(ll)
                        matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+5']=e(ll)
                        matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+6']=2*(-e(ll)+2*(`nbitemssel'+1)-2)
            }
            else {
                        qui count
                        local N=r(N)
                        di in green _col(10)"The program is estimating the parameters of the models (N=" in yellow "`N'" in green ")"
                        qui raschtest `listitemssel' `i2' , mml notest
                        matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+3']=e(ll)
                        matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+4']=2*(-e(ll)+`nbitemssel'+2)

                        local nb1:word count `listitemssel'
                        qui mmsrm `listitemssel' `i2' , part(`nb1' 1)
                        matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+5']=e(ll)
                        matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+6']=2*(-e(ll)+`nbitemssel'+4)
            }
            di in green _col(10)"Unidimensional model: " _col(32) "ll: "  in yellow _col(37) %12.4f `result`dim''[`iteration',`=`nbitemstotdim`dim''+3'] _col(52) in green "AIC: " in yellow _col(57) %12.4f `result`dim''[`iteration',`=`nbitemstotdim`dim''+4']
            di in green _col(10)"Bidimensional model: " _col(32) "ll: "   in yellow _col(37) %12.4f `result`dim''[`iteration',`=`nbitemstotdim`dim''+5'] _col(52) in green "AIC: " in yellow _col(57) %12.4f `result`dim''[`iteration',`=`nbitemstotdim`dim''+6']
            tokenize `listkernel' `varlist`dim''
            if `result`dim''[`iteration',`=`nbitemstot'+4']<=`result`dim''[`iteration',`=`nbitemstot'+6'] {
                 matrix `result`dim''[`iteration',`i']=1
                 matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+7']=1
                 local nbitemssel=`nbitemssel'+1
                 local nbitemsnosel=`nbitemsnosel'-1
                 local listitemssel `listitemssel' `rep'`i2num'
                 local listitemsselnum `listitemsselnum' `i2num'
                 matrix `affect`dim''[1,`i2num']=1
                 di _col(10) in green "The item " in yellow "``i2num''" in green " is selected in the scale " in yellow `dim'
            }
            else {
                 matrix `result`dim''[`iteration',`=`nbitemstotdim`dim''+7']=2
                 di _col(10) in green "The item "in yellow "``i2num''" in green " is not selected in the scale " in yellow `dim'

            }
     }
     return matrix result`dim' `result`dim''
     local j=`kernel'+1
     forvalues i=`=`kernel'+1'/`nbitemstot' {
            if `affect'[1,`i']==0 {
                 if `affect`dim''[1,`j']==1 {
                        matrix `affect'[1,`i']=`dim'
                 }
            local j=`j'+1
            }
     }
}
use `raschfitfile',clear

return matrix affect `affect'

end