corentinchoisy
/
ROSALI-SIM
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
main
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'main'
${ noResults }
ROSALI-SIM/Modules/ado/personal/r/old/raschtestv7761.ado

1810 lines
62 KiB
Plaintext
Raw Permalink Normal View History Unescape Escape

Setup initial file structure
9 months ago
*! version 7.6.1 26june2008
*! Jean-Benoit Hardouin
*************************************************************************************************************
* Raschtestv7: Rasch model, fit tests and graphical representations
* Corresponds to the version 7.6 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*/
*
* Needed modules :
* gammasym version 2.1 (http://freeirt.free.fr)
* gausshermite version 1 (http://freeirt.free.fr)
* geekel2d version 4.2 (http://freeirt.free.fr)
* genscore version 1.4 (http://freeirt.free.fr)
* ghquadm (findit ghquadm)
* gllamm version 2.3.11 (ssc describe gllamm)
* gllapred version 2.3.2 (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://anaqol.free.fr
* FreeIRT Project : http://freeirt.free.fr
*
* Copyright 2003-2008 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] , [ 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) GRAph v8 COMp(varname) dif(varlist) time TRace Details nold]
/***********************************************************************************************************
INTRODUCTION
***********************************************************************************************************/
if "`trace'"!="" {
di in green "*** Tests of conformity"
}
if "`v8'"=="" {
version 7.0
}
else {
version 8.0
}
local st = "$S_TIME"
local nbitems : word count `varlist'
marksample touse ,novarlist
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'"!="" {
confirm new variable `genfit' `genlt' `genscore'
if "`genfit'"!="" {
local o:word 1 of `genfit'
forvalues i=1/`nbitems' {
confirm new variable `o'``i''
}
}
}
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
}
/***********************************************************************************************************
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)
tempvar id
gen `id'=_n
qui keep `varlist' `comp' `id' `touse'
tempname rep item
tempvar score realscore
qui genscore `varlist',score(`score')
* di "LD : `ld'"
qui count if `score'==.&`touse'==1
local nbindmiss=r(N)
if "`ld'"=="" {
qui drop if `score'==.
}
forvalues i=1/`nbitems' {
rename ``i'' `rep'`i'
*qui drop if `rep'`i'==.|`rep'`i'>1
*qui replace `score'=`score'+`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"
}
tempname ll coef var beta Vbeta est
if "`method'"=="gee" {
qui geekel2d `rep'1-`rep'`nbitems',ll
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" {
qui xtlogit `rep' `rep'1-`rep'`nbitems', i(`id') nocons
matrix `est'=e(b)
matrix `est'[1,`=`nbitems'+1']=exp(`est'[1,`=`nbitems'+1']/2)
}
else if "`method'"=="cml" {
qui clogit `rep' `rep'1-`rep'`=`nbitems'-1', group(`id')
}
if "`method'"!="gee" {
matrix `coef'=e(b)
scalar `ll'=e(ll)
local nbind=e(N)/`nbitems'
matrix `var'=e(V)
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'+1',`=`nbitems'+1'])*`sig'/2
matrix `beta'=`coef'[1,1..`nbitems']
matrix `Vbeta'=`var'[1..`nbitems',1..`nbitems']
}
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'"
}
local id2=substr("`id'",1,4)
local colnames "`colnames' `id2'1:_cons"
matrix colnames `est'=`colnames'
capture gllamm `rep' `rep'1-`rep'`nbitems', i(`id') nocons family(binom) link(logit) trace from(`est') eval
tempname u
qui gllapred `u',u
tempname theta sdtheta
matrix `theta'=J(1,`=`nbitems'+1',0)
matrix `sdtheta'=J(`=`nbitems'+1',`=`nbitems'+1',0)
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''
}
}
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)
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']=V[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'=sqrt(`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
tempname 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'
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 "`test'"!="NONE" {
local conttest= "_c"
}
di
tempname itemfit globalfit
matrix `globalfit'=J(1,3,0)
if "`method'"=="cml" {
di in green "Estimation method: " in yellow "Conditional maximum likelihood (CML)"
local nbtest=`nbgroups'-1
local line=77
}
else if "`method'"=="mml"{
di in green "Estimation method: " in yellow "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
}
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)"
}
if "`method'"=="cml" {
local nbind=`nbind'+`effscore0'+`effscore`nbgroups''
local cont "_c"
matrix `itemfit'=J(`nbitems',6,0)
}
else {
local cont
matrix `itemfit'=J(`nbitems',5,0)
}
local missing=`N'-`nbind'
di in green "Number of individuals: " in yellow `N'
di in green "Number of individuals with missing values: " in yellow `nbindmiss' _c
if "`ld'"=="" {
di in green " (removed)"
}
else {
di
}
di in green "Number of individuals with nul or perfect score: " in yellow `=`effscore0'+`effscore`nbgroups'''
if "`method'"=="cml" {
di in green "Conditional log-likelihood: " in yellow %-13.4f `ll'
di in green "Log-likelihood: " in yellow %-13.4f `globalll'
}
else {
di in green "Marginal log-likelihood: " in yellow %-13.4f `ll'
di in green "Log-likelihood: " in yellow %-13.4f `globalll'
}
di
noi di in green _col(16) "Difficulty" `conttest'
if "`test'"!="NONE" {
di in green _col(58) "Standardized"
}
noi di in green _col(9) "Items" _col(16) "parameters" _col(28) "std Err." `conttest'
if "`test'"!="NONE" {
local varin=int(2/sqrt(`nbind')*100)/100
local varout=int(6/sqrt(`nbind')*100)/100
di _col(41) "`namewp'" _col(47) "df" _col(50) "p-value" _col(58) "Outfit" _col(66) "Infit" `cont'
if "`method'"=="cml" {
di in green _col(77) "U"
}
}
di in green "{hline `line'}"
forvalues i=1/`nbitems' {
noi di in yellow _col(3) %12s abbrev("``i''`fixed`i''" ,12) in yellow _col(18) %8.5f `beta`i'' _col(29) %6.5f `sd`i'' `conttest'
if "`test'"!="NONE" {
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'
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" {
di in ye _col(72) %6.3f `U`i''
matrix `itemfit'[`i',6]=`U`i''
}
}
}
di in green "{hline `line'}"
if "`test'"!="NONE" {
if "`method'"=="cml" {
local df=(`nbgroups'-2)*(`nbitems'-1)
}
else {
local df=(`nbgroups'-1)*(`nbitems'-1)-1
}
matrix `globalfit'=(`stattest',`df',1-chi2(`df',`stattest'))
noi di in green _col(6) "`descwp'" _col(32) in gr "`namewp'=" _col(36) in ye %8.3f `globalfit'[1,1] _col(46) %3.0f `globalfit'[1,2] _col(51) %6.4f `globalfit'[1,3]
if "`method'"=="cml" {
noi di in green _col(6) "Andersen LR test" _col(34) "Z=" in yellow _col(36) %8.3f `AndersenZ'[1,1] _col(46) %3.0f `AndersenZ'[1,2] _col(51) %6.4f `AndersenZ'[1,3]
}
di in green "{hline `line'}"
}
if "`method'"!="cml" {
di in green _col(10) in yellow "Sigma" _col(18) %8.5f `sig' _col(29) %6.5f `sesig'
di in green "{hline `line'}"
}
if "`method'"=="cml"&"`meandiff'"==""{
di in green "*: The difficulty parameter of this item had been fixed to 0"
}
if "`method'"=="cml"&"`meandiff'"!=""{
di in green "The mean of the difficulty parameters is fixed to 0"
}
if `smallgroup'==1&"`test'"!="NONE" {
di in green "You have groups of scores with less than 30 individuals. The tests can be invalid."
}
if "`method'"!="cml"&"`test'"=="Q" {
di in green "The Q statistics is approximated in the `method' approach. It is preferable to use the R1m test."
}
/*Tabular of the estimated values of the latent trait*/
di
di
noi di in green _col(33) "Ability" _col(60) "Expected"
noi di in green _col(17) "Group" _col(23) "Score" _col(30) "parameters" _col(42) "std Err." _col(54) "Freq." _col(63) "Score"_c
if "`method'"=="cml"&"`test'"!="NONE" {
noi di in green _col(75) "ll"
}
else {
noi di ""
}
if "`method'"=="cml"&"`test'"!="NONE" {
local line=62
}
else {
local line=51
}
di in green _col(17) "{hline `line'}"
local nonul=0
forvalues g=0/`nbgroups' {
if `g'!=0 {
di in green _col(17) "{dup `line':-}"
}
forvalues s=`liminf`g''/`limsup`g'' {
if `s'==`liminf`g'' {
local tmp `ll`g''
local gr `g'
}
else {
local tmp
local gr
}
local expscore`nonul'=0
forvalues i=1/`nbitems' {
local expscore`nonul'=`expscore`nonul''+1/(1+exp(-`theta`nonul''+`beta`i''))
}
if "`method'"=="cml" {
local format1 %8.3f
}
else {
local format1 %8.5f
}
noi di in yellow _col(17) %5s "`gr'"_col(23) %5s "`s'" _col(32) `format1' `theta`nonul'' _col(42) `format1' `sdtheta`nonul'' _col(55) %4.0f `nbrealscore`s'' _col(64) % 4.2f `expscore`nonul'' _col(68) %11.4f `tmp'
local nonul=`nonul'+1
}
}
di in green _col(17) "{hline `line'}"
/***********************************************************************************************************
DETAILS OPTION
************************************************************************************************************/
if "`details'"!="" {
forvalues g=0/`nbgroups' {
if (`g'!=0&`g'!=`nbgroups')|"`method'"!="cml" {
di
di in green "{hline 44}"
di in green "Group: " in ye "`g'" in green " from " in ye "`liminf`g''" in green " to " in ye "`limsup`g''" in green " (n=" in ye "`nbscore`g''" in green")"
di
di _col(3) "Item" _col(10) "Observed" _col(23) "Expected" _col(37) "Scaled"
di in green "{dup 44:-}"
}
if `g'!=0&`g'!=`nbgroups' {
forvalues j=1/`nbitems' {
local tmp=`d2`g''[`j',1]/sqrt(`W2`g''[`j',`j'])
di _col(3) in ye "``j''" _col(14) %4.0f `Obs2'[`j',`g'] _col(25) %6.2f `Th2'[`j',`g'] _col(36) %7.4f `tmp'
}
di in green "{dup 44:-}"
di in green "Contribution to the `namewp' statistics: " %8.4f in ye `test`g''[1,1]
}
else if "`method'"!="cml" {
if `g'==0 {
local h=0
}
else {
local h=`nbitems'
}
local tmp=abs(`nbrealscore`h''-`ci`h'')/sqrt(`ci`h'')
di in ye _col(14) %4.0f `nbrealscore`h'' _col(25) %6.2f `ci`h'' _col(36) %7.4f `tmp'
di in green "{dup 44:-}"
local tmp=`tmp'^2
di in green "Contribution to the `namewp' statistics: " %8.4f in ye `tmp'
}
}
}
/***********************************************************************************************************
OPTION ICC
************************************************************************************************************/
if "`icc'"!="" {
if "`trace'"!="" {
di in green "*** Items Characteristic Curves"
}
tempvar proba propemp propth propthb
qui replace `tl'=. if `realscore'==0|`realscore'==`nbitems'|`realscore'==.
qui gen `propemp'=.
qui gen `propth'=.
qui gen `propthb'=.
label variable `propth' "Expected ICC"
label variable `propemp' "Observed ICC"
label variable `propthb' "Expected ICC"
label variable `tl' "Latent trait"
global iccs
forvalues i=1/`nbitems' {
forvalues s=1/`=`nbitems'-1' {
qui replace `propthb'=exp(`theta`s''-`beta`i'')/(1+exp(`theta`s''-`beta`i'')) if `item'==`i'&`realscore'==`s'
}
qui replace `propemp'=.
qui replace `propth'=.
tempvar propemp`i' propth`i'
qui bysort `realscore' `item' : egen `propth`i''=mean(`propthb') if `item'==`i'
qui bysort `realscore' `item' : egen `propemp`i''=mean(`rep') if `item'==`i'
qui replace `propth'=`propth`i'' if `item'==`i'
qui replace `propemp'=`propemp`i'' if `item'==`i'
local mintl=floor(`theta1')
local maxtl=floor(`theta`=`nbitems'-1'')+1
if "`filessave'"!=""{
local saving "`dirsave'\\icc``i''"
}
if "`v8'"!="" {
graph twoway (line `propemp' `propth' `tl') if `item'==`i' , name(icc``i'',replace) ytitle("") ylabel(0(.25)1) title("Observed and Expected ICC for the item ``i''") xlabel(`mintl'(1)`maxtl') `draw'
pause
if "`filessave'"!="" {
graph save icc``i'' "`saving'" , `replace'
}
}
else {
if "`filessave'"!="" {
local saving "saving(`saving'"
if "`replace'"=="" {
local saving "`saving')"
}
else {
local saving "`saving',replace)"
}
}
qui graph `propemp' `propth' `tl' if `item'==`i' , twoway `saving' c(ll) ylabel(0(.25)1) xlabel(`mintl'(1)`maxtl') title("Observed and Expected ICC for the item ``i''")
pause
}
}
}
/***********************************************************************************************************
OPTION INFORMATION
************************************************************************************************************/
if "`information'"!="" {
if "`trace'"!="" {
di in green "*** Information graph"
}
tempfile saveinfo
qui save "`saveinfo'",replace
tempvar info latent
drop _all
qui set obs 2001
gen `latent'=((_n-1)/1001-1)*3
label variable `latent' "latent trait"
gen `info'=0
forvalues i=1/`nbitems' {
qui replace `info'=`info'+exp(`latent'-`beta`i'')/(1+exp(`latent'-`beta`i''))^2
}
local saving "`dirsave'\\information"
if "`v8'"!="" {
graph twoway (line `info' `latent') , name(information,replace) ytitle("Information") xlabel(-3(1)3) title("Information graph of the scale") `draw'
if "`filessave'"!="" {
graph save information "`saving'" , `replace'
}
pause
}
else {
if "`filessave'"!="" {
local saving "saving(`saving'"
if "`replace'"=="" {
local saving "`saving')"
}
else {
local saving "`saving',replace)"
}
}
if "`filessave'"=="" {
local saving
}
qui graph `info' `latent' , twoway `saving' c(l) xlabel(-3(1)3) title("Information graph of the scale")
pause
}
qui use "`saveinfo'",replace
}
/***********************************************************************************************************
OPTION FITGRAPH
************************************************************************************************************/
if "`fitgraph'"!="" {
if "`trace'"!="" {
di in green "*** Graphical validation of the fit"
}
*set trace on
tempname outfit meanz2 meanv2 sumd infit Voutfit outfitstd Vinfit infitstd
qui egen `outfit'=mean(`y2'),by(`id')
qui egen `Voutfit'=sum(`r'),by(`id')
qui replace `Voutfit'=`Voutfit'/(`nbitems')^2-1/`nbitems'
qui gen `outfitstd'=(`outfit'^(1/3)-1)*(3/sqrt(`Voutfit'))-sqrt(`Voutfit')/3
qui egen `meanz2'=sum(`z2'),by(`id')
qui egen `meanv2'=sum(`v2'),by(`id')
qui gen `infit'=`meanz2'/`meanv2'
qui egen `Vinfit'=sum(`e'),by(`id')
qui replace `Vinfit'=`Vinfit'/(`meanv2')^2
qui gen `infitstd'=(`infit'^(1/3)-1)*(3/sqrt(`Vinfit'))-sqrt(`Vinfit')/3
qui su `outfitstd'
local mino=floor(2*min(`r(min)',-2))/2
local maxo=ceil(2*max(`r(max)',2))/2
qui su `infitstd'
local mini=floor(2*min(`r(min)',-2))/2
local maxi=ceil(2*max(`r(max)',2))/2
if "`filessave'"!=""{
local savingo "`dirsave'\\outfitind"
local savingi "`dirsave'\\infitind"
}
if "`v8'"!="" {
*graph twoway (scatter `outfit' `id'), name(outfit,replace) yline(.6 1.4) title("Outfit indexes") xtitle("Individual indexes") ytitle("Outfit") ylabel(0(.2)2 2(.5)`maxo') `draw'
graph twoway (scatter `outfitstd' `id'), name(outfit,replace) yline(-2 2) title("Outfit indexes") xtitle("Individual indexes") ytitle("Outfit") ylabel(`mino'(.5)`maxo') `draw'
pause
*graph twoway (scatter `infit' `id'), name(infit,replace) yline(.6 1.4) title("Infit indexes") xtitle("Individual indexes") ytitle("Infit") ylabel(0(.2)2 2(.5)`maxi') `draw'
graph twoway (scatter `infitstd' `id'), name(infit,replace) yline(-2 2) title("Infit indexes") xtitle("Individual indexes") ytitle("Infit") ylabel(`mini'(.5)`maxi') `draw'
pause
if "`filessave'"!="" {
graph save outfit "`savingo'" , `replace'
graph save infit "`savingi'" , `replace'
}
}
else {
if "`filessave'"!="" {
local savingo "saving(`savingo'"
local savingi "saving(`savingi'"
if "`replace'"=="" {
local savingo "`savingo')"
local savingi "`savingi')"
}
else {
local savingo "`savingo',replace)"
local savingi "`savingi',replace)"
}
}
graph `outfitstd' `id', `savingo' twoway sy(.) c(.) yline(-2 2) title("Outfit indexes") b2title("Individual indexes") l1("") l2title("Outfit") ylabel(`mino'(.5)`maxo')
pause
graph `infitstd' `id', `savingi' twoway sy(.) c(.) yline(-2 2) title("Infit indexes") b2title("Individual indexes") l1("") l2title("Infit") ylabel(`mini'(.5)`maxi')
pause
}
drop _all
qui set obs `nbitems'
tempvar name betap outfitstdj infitstdj
qui gen str9 `name'=""
qui gen `betap'=.
qui gen `outfitstdj'=.
qui gen `infitstdj'=.
local mino=-2
local maxo=2
local mini=-2
local maxi=2
forvalues j=1/`nbitems' {
qui replace `name'="``j''" in `j'
qui replace `betap'=`beta`j'' in `j'
qui replace `outfitstdj'=`outfitstd`j'' in `j'
qui replace `infitstdj'=`infitstd`j'' in `j'
local mino=floor(min(`mino',`outfitstd`j'')*2)/2
local mini=floor(min(`mini',`infitstd`j'')*2)/2
local maxo=ceil(max(`maxo',`outfitstd`j'')*2)/2
local maxi=ceil(max(`maxi',`infitstd`j'')*2)/2
}
if "`filessave'"!=""{
local savingo "`dirsave'\\outfititems"
local savingi "`dirsave'\\infititems"
}
if "`v8'"!="" {
graph twoway (scatter `outfitstdj' `betap',name(outfititem,replace) mlabel(`name')), yline(-2 2) title("Outfit indexes") xtitle("Difficulty") ytitle("Outfit") ylabel(`mino'(.5)`maxo') `draw'
pause
graph twoway (scatter `infitstdj' `betap',name(infititem,replace) mlabel(`name')), yline(-2 2) title("Infit indexes") xtitle("Difficulty") ytitle("Infit") ylabel(`mini'(.5)`maxi') `draw'
pause
if "`filessave'"!="" {
graph save outfititem "`savingo'" , `replace'
graph save infititem "`savingi'" , `replace'
}
}
else {
if "`filessave'"!="" {
local savingo "saving(`savingo'"
local savingi "saving(`savingi'"
if "`replace'"=="" {
local savingo "`savingo')"
local savingi "`savingi')"
}
else {
local savingo "`savingo',replace)"
local savingi "`savingi',replace)"
}
}
graph `outfitstdj' `betap', `savingo' twoway c(.) yline(-2 2) title("Outfit indexes") b2title("Difficulty") l1title("") l2title("Outfit") ylabel(`mino'(.5)`maxo') sy([`name'])
pause
graph `infitstdj' `betap', `savingi' twoway c(.) yline(-2 2) title("Infit indexes") b2title("Difficulty") l1title("") l2title("Infit") ylabel(`mini'(.5)`maxi') sy([`name'])
pause
}
}
/***********************************************************************************************************
OPTION SPLITTESTS
************************************************************************************************************/
if "`splittests'"!="" {
if "`trace'"!="" {
di in green "*** Splitting tests"
}
forvalues j=1/`nbitems' {
tempname estneg`j' estpos`j'
local listitems
forvalues k=1/`nbitems' {
if `j'!=`k' {
local listitems `listitems' ``k''
}
}
qui use "`saveraschtest'",replace
qui keep if `touse'==1
qui raschtestv7 `listitems' if ``j''==0,test(NONE) meth(cml) `ld'
matrix `estneg`j''=r(beta)
local llneg=r(cll)
qui raschtestv7 `listitems' if ``j''==1,test(NONE) meth(cml) `ld'
matrix `estpos`j''=r(beta)
local llpos=r(cll)
qui raschtestv7 `listitems',test(NONE) meth(cml) `ld'
local llnegpos=r(cll)
local nbcol=colsof(`estneg`j'')
local meanneg=0
local meanpos=0
forvalues k=1/`nbcol' {
local meanneg=`meanneg'+`estneg`j''[1,`k']
local meanpos=`meanpos'+`estpos`j''[1,`k']
}
local meanneg=`meanneg'/`nbitems'
local meanpos=`meanpos'/`nbitems'
forvalues k=1/`nbcol' {
matrix `estneg`j''[1,`k']=`estneg`j''[1,`k']-`meanneg'
matrix `estpos`j''[1,`k']=`estpos`j''[1,`k']-`meanpos'
if "`method'"=="cml" {
matrix `estneg`j''=`estneg`j'',-`meanneg'
matrix `estpos`j''=`estpos`j'',-`meanpos'
}
}
drop _all
qui set obs `=`nbitems'+1'
tempvar neg pos name diag
qui gen `neg'=.
qui gen `pos'=.
qui gen str9 `name'=""
local min=`estneg`j''[1,1]
local max=`estneg`j''[1,1]
forvalues k=1/`=`nbitems'-1' {
qui replace `neg'=`estneg`j''[1,`k'] in `k'
qui replace `pos'=`estpos`j''[1,`k'] in `k'
local min=min(`min',`estneg`j''[1,`k'],`estpos`j''[1,`k'])
local max=max(`max',`estneg`j''[1,`k'],`estpos`j''[1,`k'])
local tmp:word `k' of `listitems'
qui replace `name'="`tmp'" in `k'
}
local min=floor(`min')
local max=floor(`max')+1
qui gen `diag'=.
qui replace `diag'=`min' in `nbitems'
qui replace `diag'=`max' in `=`nbitems'+1'
local Zgr=round(2*(`llneg'+`llpos'-`llnegpos'),0.001)
local Zgr=substr("`Zgr'",1,6)
local pgr=round(1-chi2(`=`nbitems'-1',`Zgr'),0.0001)
local pgr=substr("`pgr'",1,5)
local note="Z=`Zgr', df=`=`nbitems'-1' , p=`pgr'"
if "`filessave'"!=""{
local saving "`dirsave'\\split``j''"
}
if "`v8'"!="" {
graph twoway (scatter `pos' `neg' ,mlabel(`name')) (line `diag' `diag'), ytitle("Positive answer") xtitle("Negative answer") xlabel(`min'(1)`max') ylabel(`min'(1)`max') title(Splitting on ``j'') `draw' legend(off) note("Andersen Z test: `note'") name(split,replace) `draw'
pause
if "`filessave'"!="" {
graph save split "`saving'" , `replace'
}
}
else {
if "`filessave'"!="" {
local saving "saving(`saving'"
if "`replace'"=="" {
local saving "`saving')"
}
else {
local saving "`saving',replace)"
}
}
qui replace `diag'=`neg' if `diag'==.
graph `pos' `diag' `diag', `saving' c(.l) sy([`name']i) l1title("") l2title("Positive answer") b2title("Negative answer") xlabel(`min'(1)`max') ylabel(`min'(1)`max') title("Splitting on ``j'' (`note')") `draw'
pause
}
}
}
/***********************************************************************************************************
OPTION GRAPH
************************************************************************************************************/
if "`graph'"!=""&"`v8'"!="" {
if "`trace'"!="" {
di in green "*** Graph option"
}
tempvar latent2 tl2 delta2 tlth2 labdelta2 labtl2
drop _all
qui set obs 1001
gen `latent2'=(_n-501)/100
label variable `latent2' "latent trait"
qui gen `tl2'=.
qui gen `labtl2'=""
forvalues i=0/`nbitems' {
*if (`i'!=0&`i'!=`nbitems')|"`method'"!="cml" {
qui replace `tl2'=`nbrealscore`i''/`N' if round(`latent2',0.01)==round(`theta`i'',0.01)
qui replace `labtl2'="`i'" if round(`latent2',0.01)==round(`theta`i'',0.01)
*}
}
qui gen `delta2'=.
qui gen `labdelta2'=""
forvalues i=1/`nbitems' {
qui replace `delta2'=-.05 if round(`latent2',0.01)==round(`beta`i'',0.01)
qui replace `labdelta2'="``i''" if round(`latent2',0.01)==round(`beta`i'',0.01)
}
if "`method'"=="mml"|"`method'"=="gee" {
qui gen `tlth2'=1/(2*_pi*`sig')*exp(-.5*(`latent2'/`sig')^2)
label variable `tlth2' "Theorical distribution"
local graphmml line `tlth2' `latent2'
}
label variable `tl2' "Score"
label variable `delta2' "Items"
local saving "`dirsave'\\graph"
local min=-2
local max=2
forvalues j=1/`nbitems' {
if `beta`j''<`min' {
local min=`beta`j''-.5
}
if `beta`j''>`max'&`beta`j''!=. {
local max=`beta`j''+.5
}
}
if "`method'"=="cml" {
if `theta0'<`min' {
local min=`theta0'-.5
}
if `theta`=`nbitems'-1''>`max' {
local max=`theta`=`nbitems'-1''+.5
}
}
else if "`method'"!="cml" {
if `theta0'<`min'&`theta0'!=. {
local min=`theta0'-.5
}
if `theta`nbitems''>`max'&`theta`nbitems''!=. {
local max=`theta`nbitems''+.5
}
}
local min=floor(`min')
local max=floor(`max')+1
qui su `tl2'
local max2=r(max)
if "`method'"!="cml" {
qui su `tlth2'
local max2=max(`max2',`r(max)')
}
qui replace `latent2'=. if `latent2'<`min'|`latent2'>`max'
if "`v8'"!="" {
graph twoway (scatter `tl2' `latent2',mlabel(`labtl2')) (dropline `delta2' `latent2',mlabel(`labdelta2') mlabposition(6)) (`graphmml'), name(graph,replace) ytitle("") xlabel(`min'(1)`max') ylabel(-.1 0(0.05)`max2') title("Individuals/items representations") `draw'
if "`filessave'"!="" {
graph save graph "`saving'" , `replace'
}
pause
}
}
else if "`graph'"!=""&"`v8'"=="" {
di in ye "The graph option is not available with Stata 7"
}
/***********************************************************************************************************
COMPARISON OF TWO POPULATION
************************************************************************************************************/
if "`comp'"!=""&"`method'"!="cml" {
di
di _col(4) "{hline 30}"
di _col(4) in green "Comparison of two populations"
di _col(4) "{hline 30}"
di
di _col(4) in green "`comp'==`mincomp':" _col(20) "N=" _col(22) in yellow %7.0f `Nmin' _col(30) in green "mean= " _col(35) in yellow %8.4f `meanmin' _col(45) in green "variance=" _col(55) in yellow %8.4f `varmin'
di _col(4) in green "`comp'==`maxcomp':" _col(20) "N=" _col(22) in yellow %7.0f `Nmax' _col(30) in green "mean= " _col(35) in yellow %8.4f `meanmax' _col(45) in green "variance=" _col(55) in yellow %8.4f `varmax'
di _col(33) in green "Z= " _col(35) in yellow %8.4f `Zcomp' _col(47) in green "pvalue= " _col(57) in yellow %6.4f `pvalue'
}
/***********************************************************************************************************
TEST DIF
************************************************************************************************************/
if "`method'"=="cml"&"`dif'"!=""&"`test'"!="NONE" {
if "`trace'"!="" {
di in green "*** Tests of DIF"
}
local ssll=0
tempname DIFfile
qui save "`DIFfile'", replace
qui use "`saveraschtest'",replace
qui keep if `touse'==1
unab list:`dif'
local nbdif:word count `list'
tempname DIF
matrix define `DIF'=J(`nbdif',4,0)
local count=1
di
di _col(4) in green "{hline 45}"
di _col(4) in green "Test of Differential Item Functioning (DIF)"
di _col(4) in green "{hline 45}"
foreach j in `list' {
qui inspect `j'
local nbdif=r(N_unique)
qui su `j'
local maxdif=r(max)
if `nbdif'>10&`maxdif'<=10 {
di in ye "The {hi:dif} option is available with variables containing 10 or less modalities (coded from 0 or 1 to k<10)."
di in ye "The variable `j' (`nbdif' modalities) is omitted."
}
else {
local ssll=0
forvalues i=0/10 {
qui count if `j'==`i'
local effdif=r(N)
if `effdif'>0 {
qui raschtestv7 `varlist' if `j'==`i', test(NONE)
local ll`i'=r(cll)
local ssll=`ssll'+(`ll`i'')
}
}
local Z=2*`ssll'-2*`ll'
local ddl=(`nbdif'-1)*(`nbitems'-1)
matrix define `DIF'[`count',1]=(`count',`Z',`ddl',1-chi2(`ddl',`Z'))
di
di _col(4) in green "Variable:" in ye " `j' " in green " Number of groups: " in ye "`nbdif'"
di _col(4) in green "LR Test Z=" in ye %8.3f `Z' in gr " ddl=" in ye %4.0f `ddl' in gr " p-value=" in ye %6.4f `=1-chi2(`ddl',`Z')'
di
}
local ++count
}
use "`DIFfile'",replace
}
/***********************************************************************************************************
RETURN
************************************************************************************************************/
if "`trace'"!="" {
di in green "*** Returns"
}
local colnametheta
forvalues i=0/`nbgroups' {
local colnametheta `colnametheta' score_`i'
}
return clear
if `nbitems'>=3&"`test'"!="NONE" {
matrix colnames `globalfit'=`namewp' df p
matrix rownames `globalfit'=`method'
matrix rownames `itemfit'=`varlist'
matrix roweq `itemfit'=`method'
if "`method'"=="cml" {
matrix colnames `itemfit'=`namewp' df p outfit infit U
matrix rownames `AndersenZ'=`method'
matrix colnames `AndersenZ'=Z df p
return matrix AndersenZ=`AndersenZ'
if "`dif'"!="" {
matrix rownames `DIF'=`list'
matrix colnames `DIF'=num Z df p
return matrix DIF=`DIF'
}
}
else {
matrix colnames `itemfit'=`namewp' df p outfit infit
}
return matrix itemFit=`itemfit'
return matrix globalFit=`globalfit'
}
matrix colnames `theta'=`colnametheta'
matrix rownames `theta'=theta
return matrix theta `theta'
matrix colnames `sdtheta'=`colnametheta'
matrix rownames `sdtheta'=`colnametheta'
return matrix Vartheta `sdtheta'
local varlist2
matrix coleq `beta'=`method'
matrix coleq `Vbeta'=`method'
matrix roweq `Vbeta'=`method'
if "`method'"=="cml" {
forvalues i=1/`=`nbitems'-1' {
local varlist2 `varlist2' ``i''
}
return scalar cll=`ll'
return scalar ll=`globalll'
local AIC=-2*`globalll'+2*(`nbitems'-1)
}
else {
return scalar ll=`ll'
local varlist2 `varlist'
return scalar sigma=`sig'
return scalar sesigma=`sesig'
local AIC=-2*`ll'+2*(`nbitems'+1)
}
if "`comp'"!="" {
return scalar Zcomp=`Zcomp'
return scalar pZcomp=`pvalue'
}
matrix colnames `beta'=`varlist2'
matrix rownames `beta'=beta
return matrix beta `beta'
matrix colnames `Vbeta'=`varlist2'
matrix rownames `Vbeta'=`varlist2'
return matrix Varbeta `Vbeta'
return scalar AIC=`AIC'
return scalar N=`nbind'
if "`pause'"!="" {
pause off
}
drop _all
restore,not
use "`saveraschtest'"
/***********************************************************************************************************
CREATE EVENTUAL NEW VARIABLES
************************************************************************************************************/
if "`genfit'"!=""|"`genlt'"!=""|"`genscore'"!="" {
tempname genlt2 genscore2 outfit2 infit2
qui gen `score'=0 `if'
forvalues i=1/`nbitems' {
qui replace `score'=`score'+``i'' `if'
}
qui gen `genlt2'=. `if'
forvalues i=0/`nbitems' {
qui replace `genlt2'=`theta`i'' `if'&`score'==`i'
}
if "`genscore'"!="" {
qui gen `genscore'=`score' `if'
}
if "`genlt'"!="" {
qui gen `genlt'=`genlt2' `if'
}
if "`genfit'"!="" {
local outfit:word 1 of `genfit'
local infit:word 2 of `genfit'
qui gen `outfit'=0 `if'
qui gen `infit'=0 `if'
tempname infit1 infit2
qui gen `infit1'=0 `if'
qui gen `infit2'=0 `if'
forvalues j=1/`nbitems' {
tempname pi`j'
qui gen `pi`j''=exp(`genlt2'-`beta`j'')/(1+exp(`genlt2'-`beta`j'')) `if'
qui replace `pi`j''=0 `if'
forvalues s=1/`nbitems' {
qui replace `pi`j''=`Pi'[`j',`s'] `if'&`score'==`s'
}
qui gen `outfit'``j''=(``j''-`pi`j'')^2/(`pi`j''*(1-`pi`j'')) `if'
qui replace `infit1'=`infit1'+(``j''-`pi`j'')^2 `if'
qui replace `infit2'=`infit2'+(`pi`j''*(1-`pi`j'')) `if'
qui replace `outfit'=`outfit'+`outfit'``j''/`nbitems' `if'
}
qui replace `infit'=`infit1'/`infit2' `if'
}
}
if "`trace'"!=""|"`time'"!="" {
capture qui elapse `st'
di in green "** Time : " in yellow "$S_elap " in green "seconds"
}
end