*! version 8.7 28 May 2013 *! Jean-Benoit Hardouin ************************************************************************************************************* * Raschtestv7: Rasch model, fit tests and graphical representations * Corresponds to the version 8 of Raschtest (http://freeirt.org) for Stata7 * * Historic: * Version 2.1 (2003-07-10): Jean-Benoit Hardouin * Version 3.1 (2004-01-02) : Jean-Benoit Hardouin * Version 7.1 (2005-03-30) : Jean-Benoit Hardouin * Version 7.2.1 (2005-05-21) : Jean-Benoit Hardouin * Version 7.3.1 (2005-07-02) : Jean-Benoit Hardouin * Version 7.3.2 (2005-10-05) : Jean-Benoit Hardouin /*correction for the *genlt* option with *method(cml)**/ * Version 7.3.3 (2005-12-21) : Jean-Benoit Hardouin /*correction if scores are not represented with method(mml) or method(gee)*/ * Version 7.4.1 (2006-01-17) : Jean-Benoit Hardouin /*Standardized OUTFIT and INFIT, DIF option*/ * Version 7.4.2 (2006-03-21) : Jean-Benoit Hardouin /*Avoids a bug when the temporary files are saved in a directory containing special characters*/ * Version 7.4.3 (2006-11-15) : Jean-Benoit Hardouin /*Avoids a bug when the directory defined in the dirsave option contains special characters*/ * Version 7.5.1 (2007-04-20) : Jean-Benoit Hardouin /*Return results of the comp option*/ * Version 7.6.1 (2008-06-26) : Jean-Benoit Hardouin /*Nold option, corrects a bug with the genlt option*/ * Version 7.6.2 (2009-03-01) : Jean-Benoit Hardouin /*Correction of the SE of the latent trait, graph option, correction of a bug with the if option, PSI in MML*/ * Version 8.1 (2009-06-24) : Jean-Benoit Hardouin /*DIFFICULTIES option, COVARIATES option*/ * Version 8.2 (2009-07-15) : Jean-Benoit Hardouin /*Correction of a bug with CML*/ * Version 8.3 (2009-12-19) : Jean-Benoit Hardouin /*correction of a bug with DIFFICULTIES and COVARIATES options together*/ * Version 8.4 (2010-06-15) : Jean-Benoit Hardouin /*GENRES option */ * Version 8.5 (2011-12-20) : Jean-Benoit Hardouin /*Correction for the ss1 and ss3 suboptions of the COVARIATES option */ * Version 8.6 (2012-11-20) : Jean-Benoit Hardouin /*HTML option*/ * Version 8.7 (2013-05-28) : Jean-Benoit Hardouin /*Correction of a bug in the covariates option with ss1 and ss3*/ * * Needed modules : * gammasym version 2.2 (http://www.freeirt.org) * gausshermite version 1 (http://www.freeirt.org) * geekel2d version 4.3 (http://www.freeirt.org) * genscore version 1.4 (http://www.freeirt.org) * ghquadm (findit ghquadm) * gllamm version 2.3.20 (ssc describe gllamm) * gllapred version 2.3.8 (ssc describe gllapred) * elapse (ssc describe elapse) * * Jean-benoit Hardouin - Department of Biomathematics and Biostatistics - University of Nantes - France * EA 4275 "Biostatistics, Clinical Research and Subjective Measures in Health Sciences" * jean-benoit.hardouin@univ-nantes.fr * * News about this program : http://www.anaqol.org * FreeIRT Project : http://www.freeirt.org * * Copyright 2003-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 ************************************************************************************************************/ /*********************************************************************************************************** DEFINITION / SYNTAX ***********************************************************************************************************/ program define raschtestv7,rclass syntax varlist(min=1 numeric) [if] [in] , ID(varname) [HTML(string) MEANdiff DIRsave(string) FILESsave nodraw PAUse REPlace ICC INFormation SPLITtests FITgraph Method(string) group(numlist >0 ascen) AUTOGroup Test(string) q2 GENLT(string) GENSCOre(string) GENFIT(string) GENRES(string) GRAph v8 COMp(varname) dif(varlist) time TRace Details nold iterate(int 200) DIFFiculties(string) COVariates(string)] local covariables `covariates' if "`if'"=="" { local if "if 1" } /*********************************************************************************************************** INTRODUCTION ***********************************************************************************************************/ if "`trace'"!="" { di in green "*** Tests of conformity" } if "`v8'"=="" { version 7.0 } else { version 8.0 } if "`html'"!="" { set scheme sj local htmlregion "graphregion(fcolor(white) ifcolor(white))" local draw } local st = "$S_TIME" local nbitems : word count `varlist' marksample touse ,novarlist if "`covariables'"!="" { tokenize "`covariables'",parse(" ,") local i=1 local tcovariables while "``i''"!=","&"``i''"!="" { local covariable`i' ``i'' local tcovariables `tcovariables' ``i'' local ++i } local nbcovariables=`i'-1 local ss1 local ss3 if "``i''"=="," { forvalues j=`=`i'+1'/`=`i'+2' { if "``j''"=="ss1" { local ss1 ss1 } else if "``j''"=="ss3" { local ss3 ss3 } else if "``j''"=="" { * di in green "option `j' vide" } else { di in red "Invalid option in the {hi:covariables} option" error 198 } } local i=`i'+3 if "``i''"!="" { di in red "There is too much options in the {hi:covariates} option." error 198 } } local covariables `tcovariables' } else { local nbcovariables=0 } isid `id' tokenize `varlist' local bad forvalues i=1/`nbitems' { qui count if ``i''!=0&``i''!=1&``i''!=. local N=r(N) if `N'>0 { local bad `bad' ``i'' } } if "`bad'"!="" { di in red "The item(s) {hi:`bad'} is(are) not binary item(s) (with responses 0 or 1)." exit } if "`method'"=="" { local method cml } if "`test'"=="" { local test R } local method=lower("`method'") local test=upper("`test'") if "`method'"!="mml"&"`method'"!="cml"&"`method'"!="gee" { di in red "Uncorrect method option." error 198 exit } if "`test'"!="R"&"`test'"!="Q"&"`test'"!="WP"&"`test'"!="NONE" { di in red "Uncorrect test option." error 198 exit } if "`genfit'"!="" { local nbwordgenfit:word count `genfit' if `nbwordgenfit'!=2 { di in red "Uncorrect genfit option." di in red "This option must contain exactly two words" error 198 exit } } if "`genfit'"!=""|"`genlt'"!=""|"`genscore'"!=""|"`genres'"!="" { capture confirm new variable `genfit' `genscore' `genres' if "`genfit'"!="" { local o:word 1 of `genfit' forvalues i=1/`nbitems' { confirm new variable `o'``i'' } } if "`genlt'"!="" { local 0 `genlt' gettoken left 0: 0,parse(",") gettoken right 0: 0,parse(",") local genlt `left' local replacegenlt `0' local length=length("`replacegenlt'") if `length'>=3&substr("`replacegenlt'",1,`length')==substr("replace",1,`length') { local replacegenlt replace } if "`replacegenlt'"!="replace"&"`replacegenlt'"!="" { di "Not allowed option in the {hi:genlt} option" error 198 } else if "`replacegenlt'"=="" { confirm new variable `genlt' } } } preserve tempfile saveraschtest qui save `saveraschtest' qui keep if `touse'==1 if "`autogroup'"!=""&"`group'"!="" { di in green "The autogroup and the group options cannot be defined in the same time" di in green "Only the group option is retained." local autogroup } if "`autogroup'"!="" { tempvar autoscore qui genscore `varlist',score(`autoscore') tempname matscore tmp matrix `matscore'=J(`=`nbitems'-1',3,0) forvalues i=1/`=`nbitems'-1' { matrix `matscore'[`i',1]=`i' matrix `matscore'[`i',2]=`i' qui count if `autoscore'==`i' matrix `matscore'[`i',3]=r(N) } local stop=0 local j=0 while `j'<=`=`nbitems'-3'&`stop'!=1 { local j=`j'+1 local scoretogroup=99999999 local rowtogroup1=0 local rowtogroup2=0 local stop=1 forvalues i=1/`=`nbitems'-`j'' { if `matscore'[`i',1]>=0 & `matscore'[`i',3]<30 & `matscore'[`i',3]<`scoretogroup' { local scoretogroup=`matscore'[`i',3] local rowtogroup1=`i' local stop=0 } } if `stop'!=1 { if `rowtogroup1'>1&`rowtogroup1'<`=`nbitems'-`j'' { if `matscore'[`=`rowtogroup1'-1',1]<`matscore'[`=`rowtogroup1'+1',1] { local rowtogroup2=`rowtogroup1' local rowtogroup1=`rowtogroup1'-1 } else { local rowtogroup2=`rowtogroup1'+1 } } else if `rowtogroup1'==1 { local rowtogroup2=2 } else if `rowtogroup1'==`=`nbitems'-`j'' { local rowtogroup2=`nbitems'-`j' local rowtogroup1=`nbitems'-`j'-1 } matrix `tmp'=`matscore' matrix `matscore'=J(`=`nbitems'-`j'',3,0) if `rowtogroup1'!=1 { matrix `matscore'[1,1]=`tmp'[1..`=`rowtogroup1'-1',1..3] } matrix `matscore'[`rowtogroup1',1]=`tmp'[`rowtogroup1',1] matrix `matscore'[`rowtogroup1',2]=`tmp'[`rowtogroup2',2] matrix `matscore'[`rowtogroup1',3]=`tmp'[`rowtogroup1',3]+`tmp'[`rowtogroup2',3] if `rowtogroup2'!=`=`nbitems'-`j'' { matrix `matscore'[`rowtogroup2',1]=`tmp'[`=`rowtogroup2'+1'..`=`nbitems'-`j'',1..3] } } } local nbrows=rowsof(`matscore')-1 local thresholds forvalues i=1/`nbrows' { local tmp=`matscore'[`i',2] local thresholds `thresholds' `tmp' } local group `thresholds' } if "`group'"=="" { forvalues i=1/`=`nbitems'-1' { local group "`group' `i'" } } local nbgroups:word count `group' local groupmax:word `nbgroups' of `group' if `groupmax'>=`nbitems' { di in red "You cannot form a group with the higher possible score." di in red "The higher possible value of the group option is `=`nbitems'-1'." di in red "Please correct your group option." error 198 exit } else { if `groupmax'!=`=`nbitems'-1' { local group "`group' `=`nbitems'-1'" local nbgroups=`nbgroups'+1 } } local nbgroups=`nbgroups'+1 if "`dirsave'"!=""&"`filessave'"=="" { di in ye "If you want to save yours graphs, use the filessave option" } if "`filessave'"!="" { if "`dirsave'"=="" { local dirsave "`c(pwd)'" } di in ye "The graphs files will be saved in `dirsave'" } if "`dirsave'"!="" { if "`filesave'"=="" { local filesave "filesave" } } if "`pause'"!=""&"`draw'"=="" { pause on } if "`difficulties'"!="" { capture confirm matrix `difficulties' if _rc!=0 { di in red "The vector `difficulties' defined in the {hi:difficulties} option does not exist." error 198 exit } } if "`difficulties'"!=""&"`method'"!="mml" { di in red "The {hi:difficulties} option can be defined only for MML method." error 198 exit } if "`covariables'"!=""&"`method'"!="mml" { di in red "The {hi:covariables} option can be defined only for MML method." error 198 exit } if "`covariables'"!=""&"`comp'"!="" { di in red "The {hi:covariables} and {hi:comp} options can be defined jointly." error 198 exit } /*********************************************************************************************************** POSSIBLE TEST ***********************************************************************************************************/ if ("`icc'"!=""|"`fitgraph'"!=""|"`splittest'"!=""|"`details'"!="")&"`test'"=="NONE" { di in green "You cannot use the {hi:details}, {hi:icc}, {hi:fitgraph} or {hi:splittest} options if you use {hi:test(none)}." di in green "These options are ignored." local details local icc local fitgraph local splittest } if "`comp'"!=""&"`method'"=="cml" { di in green "You cannot compare two populations ({hi:comp} option) with the CML estimation method" di in green "This option is ignored." local comp } if "`method'"!="cml"&"`test'"=="WP" { di in green "The Wright-Panchapakesan test is not authorized with MML or GEE." di in green "The WP tests are replaced by Van den Wollenberg Q tests." local test Q } if "`method'"=="gee"&"`ld'"!="" { di in green "You cannot use the {hi:nold} option with the GEE method of estimation" di in green "This option is ignored." local ld } if "`test'"==""|"`test'"=="R" { local test "R" if "`method'"=="cml" { local namewp "R1c" local descwp "R1c test" } else { local namewp="R1m" local descwp "R1m test" } } if `nbitems'>999999|"`test'"=="WP" { local namewp " Y" local descwp "Wright-Panchapakesan Y test" local q2 } else if "`test'"=="Q" { local namewp " Q1" local descwp "Van den Wollenberg Q1 test" } if "`method'"!="cml"&"`meandiff'"!="" { di in green "The {hi:meandiff} option is not available with MML or GEE." di in green "This option is ignored." local meandiff } if "`method'"!="cml"&"`splittests'"!="" { di in green "The {hi:splittests} option is not available with MML or GEE." di in green "This option is ignored." local splittests } if "`method'"!="cml"&"`dif'"!="" { di in green "The {hi:dif} option is not available with MML or GEE." di in green "This option is ignored." local dif } /*********************************************************************************************************** SCORES AND GROUPS ************************************************************************************************************/ qui count if `touse'==1 local N=r(N) qui keep `varlist' `comp' `covariables' `id' `touse' tempname rep item tempvar score realscore qui genscore `varlist',score(`score') qui count if `score'==.&`touse'==1 local nbindmiss=r(N) if "`ld'"=="" { qui drop if `score'==. } forvalues i=1/`nbitems' { rename ``i'' `rep'`i' } local liminf0=0 local limsup0=0 local liminf`nbgroups'=`nbitems' local limsup`nbgroups'=`nbitems' local recode forvalues i=1/`=`nbgroups'-1' { if `i'!= 1{ local liminf`i' : word `=`i'-1' of `group' } else { local liminf1=0 } local liminf`i'=`liminf`i''+1 local limsup`i':word `i' of `group' local recode "`recode' `liminf`i''/`limsup`i''=`i'" } qui gen `realscore'=`score' qui recode `score' `recode' `nbitems'=`nbgroups' local smallgroup=0 forvalues i=0/`nbgroups' { qui count if `score'==`i' local effscore`i'=r(N) if `i'!=0&`i'!=`nbgroups'&`effscore`i''<30 { local smallgroup=1 } } /*********************************************************************************************************** ESTIMATION OF THE DIFFICULTY PARAMETERS ************************************************************************************************************/ if "`trace'"!="" { di in green "*** Estimation of the difficulty parameters" } if "`covariables'"!=""&"`difficulties'"=="" { forvalues i=1/`nbcovariables' { qui su `covariable`i'' local mean`covariable`i''=r(mean) qui replace `covariable`i''=`covariable`i''-`mean`covariable`i''' } } tempname ll coef var beta Vbeta est if "`method'"=="gee" { qui geekel2d `rep'1-`rep'`nbitems',ll local nbobserv=r(N) if `r(error)'==1 { di "The variance of the latent trait probably is too high to made the GEE a efficient method to estimate the parameters." error 499 exit } scalar `ll'=r(ll) local nbind=r(N) matrix `coef'=r(b) matrix `est'=`coef' matrix `est'[1,`=`nbitems'+1']=sqrt(`est'[1,`=`nbitems'+1']) matrix `var'=r(V) matrix `beta'=`coef'[1,1..`nbitems'] matrix `Vbeta'=`var'[1..`nbitems',1..`nbitems'] local sig=sqrt(`coef'[1,`=`nbitems'+1']) local sesig=sqrt(`var'[`=`nbitems'+1',`=`nbitems'+1'])/sqrt(`nbind')/(`sig'*2) } qui reshape long `rep', i(`id') j(`item') tempvar diff tl gen `diff'=0 forvalues i=1/`nbitems' { qui gen `rep'`i'=`item'==`i' qui replace `rep'`i'=-`rep'`i' } if "`method'"=="mml" { if "`difficulties'"==""{ qui xtlogit `rep' `rep'1-`rep'`nbitems' `covariables', i(`id') nocons iterate(`iterate') matrix `est'=e(b) matrix `est'[1,`=`nbitems'+`nbcovariables'+1']=exp(`est'[1,`=`nbitems'+`nbcovariables'+1']/2) } else { tempname offset b qui gen `offset'=0 forvalues i=1/`nbitems' { qui replace `offset'=-(`difficulties'[1,`i']) if `item'==`i' } qui gllamm `rep' `covariables' , offset(`offset') i(`id') nocons iterate(`iterate') link(logit) fam(bin) matrix `b'=e(b) matrix `est'=`difficulties',`b' matrix `est'[1,`=`nbitems'+`nbcovariables'+1']=exp(`b'[1,`=`nbcovariables'+1']/2) } } else if "`method'"=="cml" { qui clogit `rep' `rep'1-`rep'`=`nbitems'-1' , group(`id') } if "`method'"!="gee" { if "`difficulties'"=="" { matrix `coef'=e(b) matrix `var'=e(V) } else { matrix `coef'=`difficulties' matrix `var'=J(`=`nbitems'+`nbcovariables'+1',`=`nbitems'+`nbcovariables'+1',.) matrix `var'[`=`nbitems'+1',`=`nbitems'+1']=e(V) } scalar `ll'=e(ll) local nbind=e(N)/`nbitems' local nbobserv=e(N) *di "nombre d'obs:`nbobserv'" if "`meandiff'"!="" { matrix `var'=J(`nbitems',`nbitems',.) matrix `coef'=J(1,`nbitems',.) local param local lin `rep'1 forvalues j=2/`=`nbitems'-1' { local lin `lin'+`rep'`j' } local lin (`lin')/`nbitems' forvalues j=1/`=`nbitems'-1' { qui lincom `rep'`j'-`lin' matrix `coef'[1,`j']=`r(estimate)' matrix `var'[`j',`j']=`r(se)'^2 } qui lincom -`lin' matrix `coef'[1,`nbitems']=`r(estimate)' matrix `var'[`nbitems',`nbitems']=`r(se)'^2 } if "`method'"=="mml" { local sig=e(sigma_u) local sesig=sqrt(`var'[`=`nbitems'+`nbcovariables'+1',`=`nbitems'+`nbcovariables'+1'])*`sig'/2 tempname betacov Vbetacov if "`difficulties'"=="" { matrix `beta'=`coef'[1,1..`nbitems'] matrix `Vbeta'=`var'[1..`nbitems',1..`nbitems'] if "`covariables'"!="" { matrix `betacov'=`coef'[1,`=`nbitems'+1'..`=`nbitems'+`nbcovariables''] matrix `Vbetacov'=`var'[`=`nbitems'+1'..`=`nbitems'+`nbcovariables'',`=`nbitems'+1'..`=`nbitems'+`nbcovariables''] } local sig=e(sigma_u) local sesig=sqrt(`var'[`=`nbitems'+`nbcovariables'+1',`=`nbitems'+`nbcovariables'+1'])*`sig'/2 } else { tempname tmp matrix `tmp'=e(b) matrix `beta'=`difficulties' if "`covariables'"!="" { matrix `betacov'=`tmp'[1,1..`nbcovariables'] } local sig=`tmp'[1,`=`nbcovariables'+1'] matrix `tmp'=e(V) local sesig=sqrt(`tmp'[`=`nbcovariables'+1',`=`nbcovariables'+1']) matrix `Vbeta'=J(`nbitems',`nbitems',.) if "`covariables'"!="" { matrix `Vbetacov'=`tmp'[1..`nbcovariables',1..`nbcovariables'] } } } else if "`method'"=="cml"&"`meandiff'"==""{ matrix `beta'=`coef'[1,1..`=`nbitems'-1'] matrix `Vbeta'=`var'[1..`=`nbitems'-1',1..`=`nbitems'-1'] } else if "`method'"=="cml"&"`meandiff'"!=""{ matrix `beta'=`coef' matrix `Vbeta'=`var' } } if ("`method'"=="mml"|"`method'"=="gee") { local colnames forvalues i=1/`nbitems' { local colnames "`colnames' `rep':`rep'`i'" } forvalues i=1/`nbcovariables'{ local tmp:word `i' of `covariables' local colnames "`colnames' `rep':`tmp'" } local id2=substr("`id'",1,4) local colnames "`colnames' `id2'1:_cons" matrix colnames `est'=`colnames' qui gllamm `rep' `rep'1-`rep'`nbitems' `covariables', i(`id') allc nocons family(binom) link(logit) trace from(`est') eval tempname u qui gllapred `u',u fsample forvalues i=1/`nbcovariables' { local tmp:word `i' of `covariables' qui replace `u'm1=`u'm1+`betacov'[1,`i']*(`tmp')/*+`mean`covariable`i''')*/ } tempname theta sdtheta matrix `theta'=J(1,`=`nbitems'+1',0) matrix `sdtheta'=J(`=`nbitems'+1',`=`nbitems'+1',0) qui su `u'm1 if `rep'1==-1 local vartheta=r(Var) qui gen `u's2=`u's1^2 qui su `u's2 if `rep'1==-1 local meanse2=r(mean) local psi=1-`meanse2'/(`meanse2'+`vartheta') forvalues s=0/`nbitems' { qui su `u'm1 if `realscore'==`s' local theta`s'=r(mean) qui su `u's1 if `realscore'==`s' local sdtheta`s'=r(mean) matrix `theta'[1,`=`s'+1']=`theta`s'' matrix `sdtheta'[`=`s'+1',`=`s'+1']=`sdtheta`s'' } if "`genlt'"!="" { tempfile ltsave ltsavetmp qui save `ltsavetmp', replace qui keep if `rep'1==-1 qui keep `u'm1 `u's1 `id' qui sort `id' qui save `ltsave' qui use `ltsavetmp' } } forvalues i=1/`nbitems' { if `i'==`nbitems'&"`method'"=="cml"&"`meandiff'"=="" { local beta`i'=0 local sd`i' . local fixed`i' "*" } else { local beta`i'=`coef'[1,`i'] local sd`i'=sqrt(`var'[`i',`i']) } qui replace `diff'=-`beta`i'' if `item'==`i' } /*********************************************************************************************************** COMPARISON OF TWO POPULATIONS ***********************************************************************************************************/ if "`comp'"!=""&"`method'"!="cml" { if "`trace'"!="" { di in green "*** Test of comparison of two populations" } qui inspect `comp' local unique=r(N_unique) if `unique'== 2 { qui su `comp' local mincomp=r(min) local maxcomp=r(max) tempname bmin bmax qui xtlogit `rep' if `comp'==`mincomp',offset(`diff') i(`id') matrix `bmin'=e(b) local meanmin=`bmin'[1,1] local varmin=e(sigma_u)^2 local llmin=e(ll) local Nmin=e(N_g) qui xtlogit `rep' if `comp'==`maxcomp',offset(`diff') i(`id') matrix `bmax'=e(b) local meanmax=`bmax'[1,1] local varmax=e(sigma_u)^2 local llmax=e(ll) local Nmax=e(N_g) local Zcomp=(`meanmin'-`meanmax')/sqrt(`varmin'/`Nmin'+`varmax'/`Nmax') local pvalue=1-norm(abs(`Zcomp')) } else { di "It is impossible to compare more than two populations" di "The comparison process is not run" local comp } } /*********************************************************************************************************** ESTIMATION OF THE ABILITY PARAMETERS / CML ************************************************************************************************************/ if `nbitems'>=2 { if "`trace'"!="" { di in green "*** Estimation of the ability parameters" } if "`method'"=="cml" { tempfile verytmp qui save `verytmp',replace drop _all qui set obs 20001 qui gen theta=(_n-10001)/1000 qui gen A=1 forvalues j=1/`nbitems' { qui gen u`j'=exp(theta-`beta`j'') qui gen p`j'=u`j'/(1+u`j') qui gen a`j'=1/u`j' qui replace A=A*a`j' qui gen i`j'=u`j'/(1+u`j')^2 qui gen j`j'=(u`j'*(1-u`j'))/(1+u`j')^3 } qui egen P=rsum(p*) qui egen I=rsum(i*) qui egen J=rsum(j*) qui gen V=1/I^2*(I+J^2)/(4*I^2) qui gen V2=1/I tempname theta sdtheta matrix `theta'=J(1,`=`nbitems'+1',0) matrix `sdtheta'=J(`=`nbitems'+1',`=`nbitems'+1',.) forvalues s=0/`nbitems' { qui gen f`s'=abs(`s'-P+.5*J/I) qui sort f`s' matrix `theta'[1,`=`s'+1']=theta[1] matrix `sdtheta'[`=`s'+1',`=`s'+1']=sqrt(V2[1]) } use "`verytmp'",replace } qui gen `tl'=0 forvalues s=0/`nbitems' { local theta`s'=`theta'[1,`=`s'+1'] qui replace `tl'=`theta`s'' if `realscore'==`s' local sdtheta`s'=`sdtheta'[`=`s'+1',`= `s'+1'] } tempname pred qui gen `pred'=log(exp(`rep'*(`tl'+`diff'))/(1+exp(`tl'+`diff'))) qui su `pred' local globalll=r(sum) local nulscore=0 forvalues i=0/`nbgroups' { qui count if `score'==`i'&`item'==1 local nbscore`i'=r(N) if `nbscore`i''==0 { local nulscore=1 } } if `nulscore' { di in green "{p}At least one group of scores concerns none individuals. Tests will be not computed.{p_end}" di in green "{p}Use the {hi:group} or the {hi:autogroup} options.{p_end}" local test "NONE" local details local icc local fitgraph local splittest } forvalues i=0/`nbitems' { qui count if `realscore'==`i'&`item'==1 local nbrealscore`i'=r(N) } } /*********************************************************************************************************** TESTS OF THE FIRST ORDER ************************************************************************************************************/ tempname Pi matrix define `Pi'=J(`nbitems',`=`nbitems'-1',0) if "`test'"!="NONE" { qui drop if `score'==. if "`trace'"!="" { di in green "*** Tests of the first order" } tempname Obs Obs2 Th Th2 matrix define `Obs'=J(`nbitems',`=`nbitems'-1',0) matrix define `Obs2'=J(`nbitems',`=`nbgroups'-1',0) matrix define `Th'=J(`nbitems',`=`nbitems'-1',0) matrix define `Th2'=J(`nbitems',`=`nbgroups'-1',0) local listofitemsc /* Estimation of the gamma symetrical functions*/ local c0 "1" local c`nbitems' "`nbitems'*x" forvalues j=1/`nbitems' { local listini`j' local listofitemsc "`listofitemsc' `beta`j''" local c0 `c0'*(1+exp(x-`beta`j'')) local c`nbitems' `c`nbitems''-`beta`j'' forvalues k=1/`nbitems' { local listini`j'k`k' if `k'!=`j' { local listini`j' "`listini`j'' `beta`k''" } forvalues l=1/`nbitems' { if `l'!=`j'&`l'!=`k' { local listini`j'k`k' "`listini`j'k`k'' `beta`l''" } } } } gammasym `listofitemsc' /*Estimation, for each value of the score s of the probability to respond to each item (Pi) and of the theorical number of positive responses (Ths) and of theorical number of positive respond per item pair (Ws)*/ forvalues s=1/`nbitems' { local denom`s'=r(gamma`s') tempname W`s' matrix define `W`s''=J(`nbitems',`nbitems',0) } tempvar prob prob2 z y v z2 y2 v2 c r q e qui gen `prob'=. qui gen `prob2'=. forvalues j=1/`nbitems' { forvalues s=1/`=`nbitems'-1' { qui count if `rep'==1&`item'==`j'&`realscore'==`s' matrix `Obs'[`j',`s']=r(N) if "`test'"!="WP" { gammasym `listini`j'' local num`j'=r(gamma`=`s'-1') if "`method'"=="cml"|"`test'"=="Q" { matrix `Pi'[`j',`s']=exp(-`beta`j'')*`num`j''/`denom`s'' } else { gausshermite exp(`s'*x)/(`c0'), sigma(`sig') local int`s'=r(int) if "`test'"=="R"&`nbrealscore`s''!=0{ local tmp=exp(-`beta`j'')*`num`j''*`int`s''*`nbind'/`nbrealscore`s'' matrix `Pi'[`j',`s']=`tmp' } else if `nbrealscore`s''==0 { matrix `Pi'[`j',`s']=0 } } if "`test'"=="R" { forvalues k=`=`j'+1'/`nbitems' { if `s'>=2 { gammasym `listini`j'k`k'' local num`j'k`k'=r(gamma`=`s'-2') if "`method'"=="cml" { matrix `W`s''[`j',`k']=`nbrealscore`s''*exp(-`beta`j'')*exp(-`beta`k'')*`num`j'k`k''/`denom`s'' } else { matrix `W`s''[`j',`k']=exp(-`beta`j'')*exp(-`beta`k'')*`num`j'k`k''*`int`s''*`nbind' } matrix `W`s''[`k',`j']=`W`s''[`j',`k'] } } } } else if "`test'"=="WP" { matrix `Pi'[`j',`s']=exp(`theta`s''-`beta`j'')/(1+exp(`theta`s''-`beta`j'')) } matrix `Th'[`j',`s']=`Pi'[`j',`s']*`nbrealscore`s'' qui replace `prob'=`Pi'[`j',`s'] if `realscore'==`s'&`item'==`j' qui replace `prob2'=exp(`theta`s''-`beta`j'')/(1+exp(`theta`s''-`beta`j'')) if `realscore'==`s'&`item'==`j' matrix `W`s''[`j',`j']=`nbrealscore`s''*`Pi'[`j',`s'] if "`test'"!="R" { matrix `W`s''[`j',`j']=`W`s''[`j',`j']*abs(1-`Pi'[`j',`s']) } } } qui gen `v2'=abs(`prob'*(1-`prob')) qui gen `z2'=(`rep'-`prob')^2 qui gen `y2'=`z2'/`v2' qui gen `c'=abs(`prob'*(1-`prob')*(`prob'^3+(1-`prob')^3)) qui gen `r'=`c'/(`v2')^2 qui gen `e'=`c'-(`v2')^2 forvalues j=1/`nbitems' { qui su `y2' if `item'==`j' local outfit`j'=r(mean) qui su `r' if `item'==`j' local Voutfit`j'=r(sum) local Voutfit`j'=`Voutfit`j''/(`nbind')^2-1/`nbind' local outfitstd`j'=(`outfit`j''^(1/3)-1)*(3/sqrt(`Voutfit`j''))-sqrt(`Voutfit`j'')/3 qui su `z2' if `item'==`j' local n=r(sum) qui su `v2' if `item'==`j' local d=r(sum) local infit`j'=`n'/`d' qui su `e' if `item'==`j' local sume=r(sum) local Vinfit`j'=`sume'/(`d')^2 local infitstd`j'=(`infit`j''^(1/3)-1)*(3/sqrt(`Vinfit`j''))-sqrt(`Vinfit`j'')/3 } tempname tmp stattest testitems /*Estimation, by scores group g, of the theorical number of positive response (Th2g) and of the positive respond to each items pair (W2g)*/ scalar `stattest'=0 forvalues g=1/`=`nbgroups'-1' { tempname W2`g' d2`g' matrix define `W2`g''=J(`nbitems',`nbitems',0) forvalues s=`liminf`g''/`limsup`g'' { forvalues j=1/`nbitems' { matrix `Obs2'[`j',`g']=`Obs2'[`j',`g']+`Obs'[`j',`s'] matrix `Th2'[`j',`g']=`Th2'[`j',`g']+`Th'[`j',`s'] if "`test'"=="R" { forvalues k=`=`j'+1'/`nbitems' { matrix `W2`g''[`j',`k']=`W2`g''[`j',`k']+`W`s''[`j',`k'] matrix `W2`g''[`k',`j']=`W2`g''[`j',`k'] } } matrix `W2`g''[`j',`j']=`W2`g''[`j',`j']+`W`s''[`j',`j'] } } /*Estimation of the d2g vectors*/ matrix `d2`g''=`Obs2'[1..`nbitems',`g']-`Th2'[1..`nbitems',`g'] /*Estimation of the influences on the test of each item (testitemsg matrix) and of the test statistic (stattest)*/ tempname test`g' testitems`g' matrix `test`g''=`d2`g'''*syminv(`W2`g'')*`d2`g'' capture matrix `testitems`g''=cholesky(syminv(`W2`g''))*`d2`g'' if _rc!=0 { matrix `tmp'=J(`nbitems',`nbitems',0) forvalues j=1/`nbitems' { matrix `tmp'[`j',`j']=`W2`g''[`j',`j'] } di in green "Group `g' (`liminf`g''/`limsup`g''): The weight matrix is not positive-definite, the diagonal matrix is used to estimate the contribution of the items on the `wpname' statistic" matrix list `tmp' matrix `testitems`g''=cholesky(syminv(`tmp'))*`d2`g'' } else { matrix `testitems`g''=cholesky(syminv(`W2`g''))*`d2`g'' } scalar `stattest'=`stattest'+`test`g''[1,1] } matrix `testitems'=J(`nbitems',1,0) forvalues j=1/`nbitems' { forvalues g=1/`=`nbgroups'-1' { matrix `testitems'[`j',1]=`testitems'[`j',1]+`testitems`g''[`j',1]^2 } } /*Adaptation for the Q1 statistic*/ if "`test'"=="Q" { scalar `stattest'=`stattest'*`=`nbitems'-1'/`nbitems' } /*Correction for R1m and Q1m*/ if ("`test'"=="R"|"`test'"=="Q")&("`method'"=="mml"|"`method'"=="gee") { local c`nbitems' exp(`c`nbitems'')/(`c0') local c0 1/(`c0') gausshermite `c0', sigma(`sig') local ci0=r(int)*`nbind' gausshermite `c`nbitems'',sigma(`sig') local ci`nbitems'=`nbind'*r(int) scalar `stattest'=`stattest'+(`nbrealscore0'-`ci0')^2/`ci0'+(`nbrealscore`nbitems''-`ci`nbitems'')^2/`ci`nbitems'' } /*********************************************************************************************************** TESTS U ************************************************************************************************************/ if "`method'"=="cml" { if "`trace'"!="" { di in green "*** Tests U" } local quartile=`nbind'/4 local c1=0 local n1=0 while `n1'<`quartile' { local c1=`c1'+1 local n1=`n1'+`nbrealscore`c1'' } local c2=`nbitems' local n2=0 while `n2'<`quartile' { local c2=`c2'-1 local n2=`n2'+`nbrealscore`c2'' } forvalues j=1/`nbitems' { local zu1=0 local zu2=0 forvalues s=1/`c1' { local zu1=`zu1'+`nbrealscore`s''*(`Obs'[`j',`s']/`nbrealscore`s''-`Pi'[`j',`s'])/sqrt(`nbrealscore`s''*`Pi'[`j',`s']*(1-`Pi'[`j',`s'])) } forvalues s=`c2'/`=`nbitems'-1' { local zu2=`zu2'+`nbrealscore`s''*(`Obs'[`j',`s']/`nbrealscore`s''-`Pi'[`j',`s'])/sqrt(`nbrealscore`s''*`Pi'[`j',`s']*(1-`Pi'[`j',`s'])) } local U`j'=(`zu1'-`zu2')/sqrt(`c1'+`nbitems'-`c2') } } /*********************************************************************************************************** TESTS OF THE SECOND ORDER /*undocumented in beta test*/ ************************************************************************************************************/ if "`q2'"!="" { if "`trace'"!="" { di in green "*** Tests of the second order" } tempfile Q2file qui save "`Q2file'",replace qui use "`saveraschtest'",replace qui keep if `touse'==1 gen `score'=0 forvalues i=1/`nbitems' { local Q2i`i' rename ``i'' `rep'`i' qui replace `score'=`score'+`rep'`i' } qui recode `score' `recode' forvalues i=1/`nbitems' { local Q2i`i'=0 forvalues j=`=`i'+1'/`nbitems' { local listinci`i'j`j' forvalues k=1/`nbitems' { if `k'!=`i'&`k'!=`j' { local listinci`i'j`j' "`listinci`i'j`j'' `beta`k''" } } } } local Q2tot=0 forvalues k=2/`=`nbitems'-1' { forvalues i=1/`=`nbitems'-1' { forvalues j=`=`i'+1'/`nbitems' { if `k'==1 { local num=0 } else { gammasym `listinci`i'j`j'' local num=r(gamma`=`k'-2') } local athi`i'j`j'k`k'=exp(-`beta`i'')*exp(-`beta`j'')*`num'/`denom`k'' local bthi`i'j`j'k`k'=`nbth`i's`k''-`athi`i'j`j'k`k'' local cthi`i'j`j'k`k'=`nbth`j's`k''-`athi`i'j`j'k`k'' local dthi`i'j`j'k`k'=1-`athi`i'j`j'k`k''-`bthi`i'j`j'k`k''-`cthi`i'j`j'k`k'' } } } forvalues k=1/`=`nbgroups'-1' { local Q2`k'=0 forvalues i=1/`=`nbitems'-1' { forvalues j=`=`j'+1'/`nbitems' { qui count if `rep'`i'==1&`rep'`j'==1&`score'==`k' local aempi`i'j`j'k`k'=r(N) local ath2i`i'j`j'k`k'=0 local bth2i`i'j`j'k`k'=0 local cth2i`i'j`j'k`k'=0 local dth2i`i'j`j'k`k'=0 } } if `limsup`k''!=1 { forvalues l=`liminf`k''/`limsup`k'' { forvalues i=1/`=`nbitems'-1' { forvalues j=`=`i'+1'/`nbitems' { if `l'!=1 { local ath2i`i'j`j'k`k'=`ath2i`i'j`j'k`k''+`athi`i'j`j'k`l''*`nbrealscore`l'' local bth2i`i'j`j'k`k'=`bth2i`i'j`j'k`k''+`bthi`i'j`j'k`l''*`nbrealscore`l'' local cth2i`i'j`j'k`k'=`cth2i`i'j`j'k`k''+`cthi`i'j`j'k`l''*`nbrealscore`l'' local dth2i`i'j`j'k`k'=`dth2i`i'j`j'k`k''+`dthi`i'j`j'k`l''*`nbrealscore`l'' } } } } forvalues i=1/`=`nbitems'-1' { forvalues j=`=`i'+1'/`nbitems' { local d2i`i'j`j'k`k'=(`aempi`i'j`j'k`k''-`ath2i`i'j`j'k`k'')^2 local Q2i`i'j`j'k`k'=`d2i`i'j`j'k`k''/`ath2i`i'j`j'k`k''+`d2i`i'j`j'k`k''/`bth2i`i'j`j'k`k''+`d2i`i'j`j'k`k''/`cth2i`i'j`j'k`k''+`d2i`i'j`j'k`k''/`dth2i`i'j`j'k`k'' local Q2i`i'=`Q2i`i''+`Q2i`i'j`j'k`k'' local Q2i`j'=`Q2i`j''+`Q2i`i'j`j'k`k'' local Q2`k'=`Q2`k''+`Q2i`i'j`j'k`k'' } } } local Q2tot=`Q2tot'+`Q2`k'' } forvalues i=1/`nbitems' { di in green "Item ``i'' : Q2 = `Q2i`i''" } local Q2=`Q2tot'*(`nbitems'-3)/(`nbitems'-1) di in green "Q2 = `Q2tot'" qui use "`Q2file'",replace } /*********************************************************************************************************** TEST LR Z ************************************************************************************************************/ if "`method'"=="cml" { if "`trace'"!="" { di in green "*** Tests LR of Andersen" } local ssll=0 tempfile Zfile qui save "`Zfile'", replace qui use "`saveraschtest'",replace qui keep if `touse'==1 gen `score'=0 forvalues j=1/`nbitems' { qui replace `score'=`score'+``j'' } qui recode `score' `recode' `nbitems'=`nbgroups' forvalues i=1/`=`nbgroups'-1' { if `effscore`i''>0 { qui raschtestv7 `varlist' if `score'==`i', test(NONE) `ld' id(`id') local ll`i'=r(cll) local ssll=`ssll'+(`ll`i'') } } local Z=2*`ssll'-2*`ll' use "`Zfile'",replace tempname AndersenZ matrix define `AndersenZ'=(`Z',`=(`nbgroups'-2)*(`nbitems'-1)',1-chi2(`=(`nbitems'-1)*(`nbgroups'-2)',`Z')) } } /*********************************************************************************************************** DISPLAYING RESULTS WITH TESTS ************************************************************************************************************/ if "`html'" != "" { di "" } if "`test'"!="NONE" { local conttest= "_c" } di tempname itemfit globalfit matrix `globalfit'=J(1,3,0) if "`method'"=="cml" { if "`html'" == "" { di in green "Estimation method: " in yellow "Conditional maximum likelihood (CML)" } else { di "
Estimation method: Conditional maximum likelihood (CML)
" } local nbtest=`nbgroups'-1 local line=77 } else if "`method'"=="mml"{ if "`html'" == "" { di in green "Estimation method: " in yellow "Marginal maximum likelihood (MML)" } else { di "Estimation method: Marginal maximum likelihood (MML)
" } local nbtest=`nbgroups'+1 local line=70 } else if "`method'"=="gee" { di in green "Estimation method: " in yellow "Generalized Estimating Equations (GEE)" local nbtest=`nbgroups'+1 local line=70 } if "`test'"=="NONE" { local line=35 } if "`html'" == "" { di in green "Number of items: " in yellow `nbitems' di in green "Number of groups: " in yellow `=`nbgroups'+1' `conttest' if "`test'"!="NONE" { di in green " (" in yellow "`nbtest'" in green " of them are used to compute the statistics of test)" } } else { di "Number of items: " `nbitems' "Difficulty | " if "`test'"!="NONE" { di " | Standardized | " } di "||||||
---|---|---|---|---|---|---|---|---|
Items | parameters | std Err. | " if "`test'"!="NONE" { local varin=int(2/sqrt(`nbind')*100)/100 local varout=int(6/sqrt(`nbind')*100)/100 di "`namewp' | df | p-value | Outfit | Infit | " if "`method'"=="cml" { di "U | " } } di "
" %12s abbrev("``i''`fixed`i''" ,12) " | " %8.5f `beta`i'' " | " %6.5f `sd`i'' " | " } if "`test'"!="NONE" { if "`html'" == "" { di _col(36) %8.3f `testitems'[`i',1] _col(46) %3.0f `=`nbgroups'-2' _col(51) %6.4f 1-chi2(`=`nbgroups'-2',`testitems'[`i',1]) _col(58) %6.3f `outfitstd`i'' _col(65) %6.3f `infitstd`i'' /*_col(72) %6.3f `outfit`i'' _col(79) %6.3f `infit`i''*/ `cont' } else { di "" %8.3f `testitems'[`i',1] " | " %3.0f `=`nbgroups'-2' " | " %6.4f 1-chi2(`=`nbgroups'-2',`testitems'[`i',1]) " | " %6.3f `outfitstd`i'' " | " %6.3f `infitstd`i'' " | " } matrix `itemfit'[`i',1]=`testitems'[`i',1] matrix `itemfit'[`i',2]=`=`nbgroups'-2' matrix `itemfit'[`i',3]=1-chi2(`=`nbgroups'-2',`testitems'[`i',1]) matrix `itemfit'[`i',4]=`outfitstd`i'' matrix `itemfit'[`i',5]=`infitstd`i'' if "`method'"=="cml" { if "`html'" == "" { di in ye _col(72) %6.3f `U`i'' } else { di "" %6.3f `U`i'' " | " } matrix `itemfit'[`i',6]=`U`i'' } if "`html'" != "" { di "
`descwp' | `namewp'= | " %8.3f `globalfit'[1,1] " | " %3.0f `globalfit'[1,2] " | " %6.4f `globalfit'[1,3] " | ||||
Andersen LR test | Z= | " %8.3f `AndersenZ'[1,1] " | " %3.0f `AndersenZ'[1,2] " | " %6.4f `AndersenZ'[1,3] " |
Parameters | Coef. | std Err. | z | P>|z| |
---|---|---|---|---|
Sigma | " %8.5f `sig' " | " %6.5f `sesig' " | " %6.3f `zsig' " | " %7.4f `pzsig' " |
Group | Score | Ability parameters | std Err. | Freq. | Expected Score | "
if "`method'"=="cml"&"`test'"!="NONE" {
di "ll | " } di "
---|---|---|---|---|---|---|
" %5s "`gr'" " | " %5s "`s'" " | " `format1' `theta`nonul'' " | " `format1' `sdtheta`nonul'' " | " %4.0f `nbrealscore`s'' " | " % 4.2f `expscore`nonul'' " | " %11.4f `tmp' " | " } local nonul=`nonul'+1 } } if "`html'" == "" { di in green _col(17) "{hline `line'}" } else { di "
Mean variance of the error | " %10.4f `meanse2' " |
---|---|
Estimated variance of the latent trait | " %10.4f `sig2' " |
Personal Separation Index (PSI) | " %10.4f `psi2' " |
Adjusted PSI on covariates (PSIadj) | " %10.4f `psi' " |
Item | Observed | Expected | Scaled |
---|---|---|---|
``j'' | " %4.0f `Obs2'[`j',`g'] " | " %6.2f `Th2'[`j',`g'] " | " %7.4f `tmp' " |
Contribution to the `namewp' statistics: | " %8.4f `test`g''[1,1] " | ||
" %4.0f `nbrealscore`h'' " | " %6.2f `ci`h'' " | " %7.4f `tmp' " | |
Contribution to the `namewp' statistics: | " %8.4f `tmp' " |