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.

1197 lines
41 KiB
Plaintext

program define pcmodelplus, eclass
syntax varlist [, COVariates(string) DIFficulties(string) ITerate(string) Group(string) adapt robust nodis]
/* syntaxe : liste des variables correspondant aux items, cov(liste des covariables de groupe, qualitatives codées numériquement), dif(liste de vecteurs lignes ayant le même nom que les variables d'items, correspondant aux difficultés des modalités d'items) it(nombre d'itérations maximum pour la procédure d'estimation par convergeance) adapt : utilisation de quadratures "adaptative" robust : utilisation de l'estimateur de Huber/White/sandwich pour la matrice de covariance */
preserve
tokenize `varlist' `covariates'
local nbit: word count `varlist'
local nbcov: word count `covariates'
local nbit2:word count `difficulties'
tempfile bddini
qui q memory
local matsizeini=r(matsize)
/*************************************/
/* Verifications */
/*************************************/
if "`iterate'"==""{
local it=""
}
else{
local iterateII="it(`iterate')"
}
tempvar one id item reponse obs wt x choix it covariable covverytemp
qui save `bddini', replace
if "`difficulties'"!=""{
if `nbit2'!=`nbit'{
noi di in red "Not the same number of difficuly vectors and of items, Grrr!"
use `bddini',replace
error 100
}
forvalues i=1/`nbit'{
if strpos("`varlist'",word("`difficulties'",`i'))==0{
noi di in red "The item difficulty vectors must have the same names than the item variables, Grrr! be carefull please!!!"
use `bddini',replace
error 100
}
}
}
/*************************************/
/* Estimation of the parameters */
/*************************************/
if `nbcov'!=0{
local nbtotmodacov=0
forvalues i=1/`nbcov'{
qui tab ``=`nbit'+`i''', matrow(nom) matcell(val)
local nbModCov`i'=r(r)
forvalues k=1/`nbModCov`i''{
local valmod`k'cov`i'=nom[`k',1]
local nbmod`k'cov`i'=val[`k',1]
}
local nbtotmodacov=`nbtotmodacov'+`nbModCov`i''
if `nbModCov`i''>10{
noi di in red "``=`nbit'+`i''' has too much modalities (>10). Hey, why not a continuous variable, while you are at it?"
use `bddini',replace
error 100
}
local tot`i'=r(N)
forvalues k=1/`nbModCov`i''{
local sum=0
forvalues k2=1/`nbModCov`i''{
if `k2'!=`k'{
local sum=`sum'+val[`k2',1]
}
}
local a`i'_`k'=round((`sum'-`tot`i'')/`tot`i'',0.01)
local b`i'_`k'=round((`sum')/`tot`i'',0.01)
}
}
}
qui{
keep `varlist' `covariates'
local nbdifftot=0
if "`difficulties'"==""{
forvalues i=1/`nbit'{
gen `reponse'`i' = ``i''
drop ``i''
su `reponse'`i'
local moda`i'=`=r(max)+1'
local nbdifftot=`nbdifftot'+`moda`i''-1
}
}
else{
forvalues i=1/`nbit'{
gen `reponse'`i' = ``i''
drop ``i''
local moda`i'=colsof(``i'')+1
}
}
gen `one'=1
su `one'
local Nbid=r(N)
collapse (sum) `wt'2=`one', by(`reponse'1-`reponse'`nbit' `covariates')
gen `id'=_n
reshape long `reponse', i(`id') j(`item')
drop if `reponse'==.
gen `obs'=_n
local ddlssIII=_N-1
forvalues i=1/`nbit'{
expand `moda`i'' if `item'==`i'
}
by `obs', sort: gen `x'=_n-1
gen `choix'=`reponse'==`x'
tab `item', gen(`it')
forvalues i=1/`nbit'{
forvalues g=1/`=`moda`i''-1'{
gen d_``i''_`g'=(-1)*`it'`i'*(`x'>=`g')
}
}
if "`difficulties'"!=""{
gen difficulties=0
forvalues i=1/`nbit'{
forvalues g=1/`=`moda`i''-1'{
replace difficulties=difficulties+``i''[1,`g']*d_``i''_`g'
}
}
}
if `nbcov'!=0{
forvalues i=1/`nbcov'{
gen `covverytemp'=``=`nbit'+`i'''
drop ``=`nbit'+`i'''
rename `covverytemp' ``=`nbit'+`i'''
tab ``=`nbit'+`i''', gen(``=`nbit'+`i'''__) matrow(nom)
local nbModCov`i'=r(r)
forvalues k=1/`nbModCov`i''{
replace ``=`nbit'+`i'''__`k'=``=`nbit'+`i'''__`k'*`x'
local ident`i'_`k'=nom[`k',1]
if `k'==1 & `i'==1{
local ident1=`ident`i'_`k''
}
rename ``=`nbit'+`i'''__`k' ``=`nbit'+`i'''_`ident`i'_`k''
}
}
forvalues i=1/`nbcov'{
drop ``=`nbit'+`i'''
}
}
rename `id' theta
rename `x' estimates
}
eq slope:estimates
gen obs=`obs'
gen choix=`choix'
gen wt=`wt'
if `nbcov'!=0{
if "`difficulties'"==""{
/* Recherche de valeurs initiales */
matrix a0=J(1,`=`nbdifftot'+`nbtotmodacov'-`nbcov'',0)
qui{
tempfile bddencours
save `bddencours', replace
use `bddini', replace
local countdif=1
forvalues i=1/`nbit'{
ologit ``i''
matrix ttt=e(b)
matrix a0[1,`countdif']=ttt/2
local countdif=`countdif'+e(k)
}
use `bddencours', replace
}
local colnom ""
forvalues i=1/`nbit'{
forvalues j=1/`=`moda`i''-1'{
local colnom "`colnom' d_``i''_`j'"
}
}
}
constraint drop _all
local contrainte ""
local subgroup0 "``=`nbit'+1'':`ident1_1'"
forvalues i=1/`nbcov'{
local contrainte "`contrainte' `i'"
constraint `i' ``=`nbit'+`i'''_`ident`i'_1'=0
if "`difficulties'"==""{
forvalues j=2/`nbModCov`i''{
local colnom "`colnom' ``=`nbit'+`i'''_`ident`i'_`j''"
}
}
}
if `nbcov'>1{
forvalues i=2/`nbcov'{
local subgroup0 "`subgroup0', ``=`nbit'+`i''':`ident`i'_1'"
}
}
if "`difficulties'"==""{
if wordcount("`colnom'")!=`=`nbdifftot'+`nbtotmodacov'-`nbcov''{
di in red "ERRUER MAJEUR!!! (trop de covariables, d'items, longueur de nom des variables d'items ou de covariable trop longs... bref, réduire le nombre de caractères)"
use `bddini',replace
error 100
}
matrix coln a0=`colnom'
matrix coleq a0=estimates
matrix a1=J(1,1,0)
matrix coln a1=estimates
matrix coleq a1=theta1_1
matrix a=a0,a1
}
if "`difficulties'"!=""{
gen cons=estimates
}
ereturn clear
if "`difficulties'"==""{
gllamm estimates d_`1'_1-d_``nbit''_`=`moda`nbit''-1' ``=`nbit'+1''_`ident1'-``=`nbit'+`nbcov'''_`ident`nbcov'_`nbModCov`nbcov''', i(theta) eqs(slope) link(mlogit) expand(`obs' `choix' o) weight(`wt') `adapt' `robust' nocons `iterateII' nodis from(a) constraint(`contrainte') copy
matrix estimations = e(b)
local sdmu=estimations[1,`=colsof(estimations)']
}
else{
gllamm estimates cons ``=`nbit'+1''_`ident1'-``=`nbit'+`nbcov'''_`ident`nbcov'_`nbModCov`nbcov''', offset(difficulties) i(theta) eqs(slope) link(mlogit) expand(`obs' `choix' o) weight(`wt') `iterateII' `adapt' `robust' nocons nodis constraint(`contrainte')
}
matrix eB=e(b)
local nbc=colsof(eB)
matrix eV=e(V)
/* SS type III*/
tempfile bddencours
qui save `bddencours', replace
qui use `bddini', replace
local varianceRes=eB[1,`=colsof(eB)']^2
qui{
local contv=1
local ddlssIIIc=`ddlssIII'
forvalues c=`=1+`nbit''/`=`nbcov'+`nbit''{
local ddlssIIIc=`ddlssIIIc'-(`nbModCov`=`c'-`nbit'''-1)
tab ``c'', matcell(cell) matrow(row)
local muG=0
local denommuG=0
local muGcarre=0
if "`difficulties'"==""{
forvalues cm=1/`nbModCov`=`c'-`nbit'''{
local muG=`muG'+cell[`cm',1]*eB[1,`=`nbdifftot'+`contv'']
local denommuG=`denommuG'+cell[`cm',1]
local muGcarre=`muGcarre'+(`varianceRes'+eB[1,`=`nbdifftot'+`contv'']^2)*cell[`cm',1]
local contv=`contv'+1
}
}
else{
forvalues cm=1/`nbModCov`=`c'-`nbit'''{
local muG=`muG'+cell[`cm',1]*eB[1,`=1+`contv'']
local denommuG=`denommuG'+cell[`cm',1]
local muGcarre=`muGcarre'+(`varianceRes'+eB[1,`=1+`contv'']^2)*cell[`cm',1]
local contv=`contv'+1
}
}
local Var_`c'=(`muGcarre'/`denommuG') - (`muG'/`denommuG')^2
}
local contv=1
local Varttot=`varianceRes'
forvalues c=`=1+`nbit''/`=`nbcov'+`nbit''{
local ddlssIIIs`c'=`ddlssIIIc'+(`nbModCov`=`c'-`nbit'''-1)
tab ``c'', matcell(cell) matrow(row)
local muG=0
local denommuG=0
local muGcarre=0
if "`difficulties'"==""{
forvalues cm=1/`nbModCov`=`c'-`nbit'''{
local muG=`muG'+cell[`cm',1]*eB[1,`=`nbdifftot'+`contv'']
local denommuG=`denommuG'+cell[`cm',1]
local muGcarre=`muGcarre'+(`Varttot'+eB[1,`=`nbdifftot'+`contv'']^2)*cell[`cm',1]
local contv=`contv'+1
local Varttot=(`muGcarre'/`denommuG') - (`muG'/`denommuG')^2
}
}
else{
forvalues cm=1/`nbModCov`=`c'-`nbit'''{
local muG=`muG'+cell[`cm',1]*eB[1,`=1+`contv'']
local denommuG=`denommuG'+cell[`cm',1]
local muGcarre=`muGcarre'+(`Varttot'+eB[1,`=1+`contv'']^2)*cell[`cm',1]
local contv=`contv'+1
local Varttot=(`muGcarre'/`denommuG') - (`muG'/`denommuG')^2
}
}
}
forvalues c=`=1+`nbit''/`=`nbcov'+`nbit''{
local sst3_`c'=(`Var_`c''*`ddlssIIIs`c''-`varianceRes'*`ddlssIIIc')
local VarExp_`c'=(`Var_`c''*`ddlssIIIs`c''-`varianceRes'*`ddlssIIIc')/(`Varttot'*`ddlssIII')
local df_`c'=`nbModCov`=`c'-`nbit'''-1
}
local contv=1
forvalues c=`=1+`nbit''/`=`nbcov'+`nbit''{
if "`difficulties'"==""{
forvalues cc=1/`nbModCov`=`c'-`nbit'''{
if `cc'==1{
local contv=`contv'+1
}
else{
local ESit`c'mod`cc' = `=eB[1,`=`nbdifftot'+`contv'']/sqrt(`Varttot')'
local contv=`contv'+1
}
noi di in gr ""
}
}
else{
forvalues cc=1/`nbModCov`=`c'-`nbit'''{
if `cc'==1{
local contv=`contv'+1
}
else{
/*
noi di `contv'
noi matrix list eB
*/
local ESit`c'mod`cc' = `=eB[1,`=1+`contv'']/sqrt(`Varttot')'
local contv=`contv'+1
}
noi di in gr ""
}
}
}
}
qui use `bddencours', replace
}
else{
ereturn clear
if "`difficulties'"==""{
matrix a0=J(1,`nbdifftot',0)
qui{
tempfile bddencours
save `bddencours', replace
use `bddini', replace
local countdif=1
forvalues i=1/`nbit'{
ologit ``i''
matrix ttt=e(b)
matrix a0[1,`countdif']=ttt/2
local countdif=`countdif'+e(k)
}
use `bddencours', replace
}
local colnom ""
forvalues i=1/`nbit'{
forvalues j=1/`=`moda`i''-1'{
local colnom "`colnom' d_``i''_`j'"
}
}
matrix coln a0=`colnom'
matrix coleq a0=estimates
matrix a1=J(1,1,0)
matrix coln a1=estimates
matrix coleq a1=theta1_1
matrix a=a0,a1
gllamm estimates d_`1'_1-d_``nbit''_`=`moda`nbit''-1', i(theta) eqs(slope) link(mlogit) expand(`obs' `choix' o) weight(`wt') `iterateII' `adapt' `robust' nocons nodis from(a) copy
matrix estimations = e(b)
local sdmu=estimations[1,`=colsof(estimations)']
}
else{
gen cons=estimates
gllamm estimates cons, offset(difficulties) i(theta) eqs(slope) link(mlogit) expand(`obs' `choix' o) weight(`wt') `iterateII' `adapt' `robust' nocons nodis
}
matrix eB=e(b)
local nbc=colsof(eB)
matrix eV=e(V)
}
clear
use `bddini'
if "`difficulties'"==""{
/*************************************/
/* test R1m */
/*************************************/
noi di in gr "Performing R1m test"
qui{
set matsize `matsizeini'
tempfile matU
tempname dm cell row
egen `dm'=rowmiss(`varlist')
replace `dm'=`dm'>0
tab `dm',matcell(`cell') matrow(`row')
local nbmissing=`cell'[2,1]
local percentmissing=`cell'[2,1]/(`cell'[1,1]+`cell'[2,1])
if `nbmissing'==.{
local nbmissing=0
local percentmissing=0
}
keep if `dm'==0
drop `dm'
keep `varlist' `covariates'
if `nbcov'!=0{
local mumodifie=0
local compteuri=`nbdifftot'
forvalues c=1/`nbcov'{
forvalues m=1/`nbModCov`c''{
local compteuri=`compteuri'+1
local mumodifie=`mumodifie'+(`nbmod`m'cov`c''/`Nbid')*estimations[1,`compteuri']
}
}
drop `covariates'
}
forvalues i=1/`nbit'{
rename ``i'' pre_``i''
}
forvalues i=1/`nbit'{
gen it`i'=pre_``i''
drop pre_``i''
}
forvalues i=1/`nbit'{
qui tab it`i'
local nbmoda_`i'=r(r)
}
local score "it1"
forvalues i=2/`nbit'{
local score "`score'+it`i'"
}
gen score=`score'
save `bddini'_se, replace
local N=_N
/* autogroup */
if "`group'"==""{
tab score, matcell(freq) matrow(val)
if r(N)<60{
di in red "Pas assez d'individus pour grouper de façon automatique"
}
local nbscore=r(r)
local ini=0
local pb=1
local list ""
forvalues i=1/`nbscore'{
local ini=`ini'+freq[`i',1]
if `ini'<30{
local lim`i'=""
local pb=1
}
else{
local lim`i'=" `=val[`i',1]'"
local ini=0
local pb=0
}
local list "`list' `lim`i''"
}
local list "`=itrim("`list'")'"
if `pb'==1{
local list "`=subinword("`list'",word("`list'",-1),"",.)'"
local list "`=itrim("`list' `=val[`nbscore',1]'")'"
}
if wordcount("`list'")<2{
di in red "Autogroup impossible : pas de possibilité de constituer plusieurs groupes de plus de 30 individus à partir du score"
}
local nbgroup=wordcount("`list'")
}
else{
local list "`group'"
local nbgroup=wordcount("`list'")
tab score, matcell(freq) matrow(val)
local nbscore=r(r)
}
local nbmodagpe=0
forvalues i=1/`nbit'{
local nbmodagpe=`nbmodagpe'+`nbmoda_`i''
}
gen g=.
local liminf1=-1
forvalues i=1/`nbgroup'{
local limsup`i'=real("`=word("`list'",`i')'")
if `i'>1{
local liminf`i'=real("`=word("`list'",`=`i'-1')'")
}
replace g=`i' if score<=`limsup`i'' & score>`liminf`i''
}
forvalues i=1/`nbgroup'{
su score if g==`i'
local N`i'=r(N)
}
save `bddini'_se, replace
save `bddini'_se_g, replace
use `bddini'_se, replace
local count=0
forvalues i=1/`nbit'{
forvalues j=1/`nbmoda_`i''{
if `j'==1{
local diffest`i'_`j'=0
}
else{
local count=`count'+1
local diffest`i'_`j'=round(estimations[1,`count'],0.001)+`diffest`i'_`=`j'-1''
}
local num`i'_`j' "exp(`=`j'-1'*x-`diffest`i'_`j'')"
if `j'==1{
local denom`i' "1"
}
else{
local denom`i' "`denom`i''+`num`i'_`j''"
}
}
forvalues j=1/`nbmoda_`i''{
local PCM`i'_`j' "(`num`i'_`j'')/(`denom`i'')"
}
}
local nbcomb=1
forvalues i=1/`nbit'{
local nbcomb=`nbcomb'*`nbmoda_`i''
}
clear
set obs `nbcomb'
local facteur=1
gen score=0
local order "score"
forvalues i=1/`nbit'{
gen it`i'=floor(mod(_n/`facteur',`nbmoda_`i''))
local facteur=`facteur'*`nbmoda_`i''
replace score=score+it`i'
tab it`i', gen(i`i'c)
local order "`order' it`i'"
}
order `order', first
sort score
/* Probabilité des paterns de réponse */
tab score, gen(score)
local nbscore=r(r)
gen P=1
local gdN=_N
forvalues l=1/`gdN'{
local exp ""
forvalues i=1/`nbit'{
forvalues m=1/`nbmoda_`i''{
if i`i'c`m'[`l']==1{
if "`exp'"==""{
local exp "`PCM`i'_`m''"
}
else{
local exp "`exp'*`PCM`i'_`m''"
}
}
}
}
if `nbcov'==0{
gausshermite `exp' , sigma(`sdmu') mu(0) display
}
else{
local sdgauss=sqrt(`Varttot')
gausshermite `exp' , sigma(`sdgauss') mu(`mumodifie') display
}
replace P=r(int) if _n==`l'
di "`exp'"
}
tab score
local Pscoresl=r(r)
matrix Pscores=J(`Pscoresl',2,.)
forvalues i=1/`Pscoresl'{
su P if score==`=`i'-1'
matrix Pscores[`i',2]=r(sum)
matrix Pscores[`i',1]=`i'-1
}
save `matU'_2, replace
local sort "score"
forvalues i=1/`nbit'{
local sort "`sort' it`i'"
forvalues j=1/`nbmoda_`i''{
replace i`i'c`j'=i`i'c`j'*P
}
}
forvalues i=1/`nbscore'{
replace score`i'=score`i'*P
}
save `matU'_1, replace
/* Matrices U */
local order ""
forvalues i=1/`=wordcount("`list'")'{
local nbscoreg`i'= `nbscore'
forvalues j=0/`=`nbscore'-1'{
if `j'>`limsup`i'' | `j'<=`liminf`i''{
local nbscoreg`i'=`nbscoreg`i''-1
capture drop s`=`j'+1'g`i'
}
}
local order "`order' *g`i'"
use `matU'_1, replace
keep if score<=`limsup`i'' & score>`liminf`i''
gen g=`i'
forvalues a=1/`nbit'{
forvalues b=1/`nbmoda_`a''{
capture rename i`a'c`b' i`a'c`b'g`i'
}
}
forvalues a=1/`nbscore'{
capture rename score`a' s`a'g`i'
forvalues j=1/`=wordcount("`list'")'{
if `=`a'-1'>`limsup`i'' | `=`a'-1' <= `liminf`i''{
capture drop s`a'g`j'
}
}
}
tempfile bloc`i'
sort `sort'
save `bloc`i''_1, replace
}
local order ""
forvalues i=1/`=wordcount("`list'")'{
local order "`order' *g`i'"
use `matU'_2, replace
keep if score<=`limsup`i'' & score>`liminf`i''
gen g=`i'
forvalues a=1/`nbit'{
forvalues b=1/`nbmoda_`a''{
capture rename i`a'c`b' i`a'c`b'g`i'b
}
}
forvalues a=1/`nbscore'{
capture rename score`a' s`a'g`i'b
forvalues j=1/`=wordcount("`list'")'{
if `=`a'-1'>`limsup`i'' | `=`a'-1' <= `liminf`i''{
capture drop s`a'g`j'
}
}
}
sort `sort'
save `bloc`i''_2, replace
describe
local mat_`i'=r(k)-3-`nbit'
}
/* Matrice W */
local nbscoreg0=0
forvalues g=1/`nbgroup'{
qui q memory
if `mat_`g''>r(matsize) & `mat_`g''<=r(m_matsize){
set matsize `mat_`g''
}
if `mat_`g''>r(m_matsize){
di in red "Error, not enought memory to allocate to the W matrix calculation"
}
matrix W_`g'=J(`mat_`g'',`mat_`g'',.)
use `bloc`g''_2, replace
local count=1
forvalues i=1/`nbit'{
forvalues j=1/`nbmoda_`i''{
rename i`i'c`j'g`g'b v`count'b
local count=`count'+1
}
}
forvalues s=1/`=`nbscore'+1'{
capture rename s`s'g`g'b v`count'b
if _rc ==0{
local count=`count'+1
}
}
local Ncol_`g'_=`count'-1
save `bloc`g''_2b, replace
use `bloc`g''_1, replace
local count=1
forvalues i=1/`nbit'{
forvalues j=1/`nbmoda_`i''{
rename i`i'c`j'g`g' v`count'
local count=`count'+1
}
}
forvalues s=1/`=`nbscore'+1'{
capture rename s`s'g`g' v`count'
if _rc ==0{
local count=`count'+1
}
}
local Ncol_`g'__=`count'-1
forvalues i=1/`=`count'-1'{
replace v`i'=P if round(v`i',0.000001)!=round(0,0.000001)
}
local countW_`g'=`count'-1
save `bloc`g''_1b, replace
merge `sort' using `bloc`g''_2b
drop _merge
forvalues i=1/`mat_`g''{
forvalues j=1/`mat_`g''{
gen prod=v`i'*v`j'b
qui su prod
matrix W_`g'[`i',`j']=r(sum)
drop prod
}
}
capture matrix W_`g'i=invsym(W_`g')
if _rc!=0{
capture matrix W_`g'i=inv(W_`g')
}
if _rc!=0{
di in red "Error while computing the Wg matrix"
}
mata: A=J(1,1,.)
mata: A[1,1]=rank(st_matrix("W_`g'"))
mata: st_matrix("rank_",A)
local rank`g'=rank_[1,1]
}
/* Matrice poid */
forvalues g=1/`nbgroup'{
use `bddini'_se_g, replace
keep if g==`g'
gen Pd=1
collapse (sum) nb=Pd, by(`sort')
sort `sort'
save `bddini'_se_`g'obs, replace
use `bloc`g''_2b, replace
sort `sort'
merge `sort' using `bddini'_se_`g'obs
drop _merge
replace nb=0 if nb==.
gen nbth=P*`N'
gen Di=(nb-nbth)/sqrt(`N')
sort `sort'
save `bloc`g''_2bgi, replace
if `countW_`g''>10{
set matsize `countW_`g''
}
else{
set matsize 10
}
matrix deviation_`g'=J(`countW_`g'',1,.)
use `bloc`g''_2bgi, replace
gen Pobs=nb/`N'
gen Diff=Pobs-P
forvalues j=1/`countW_`g''{
replace v`j'b=v`j'b*Diff
su v`j'b
matrix deviation_`g'[`j', 1]=r(sum)*sqrt(`N')
}
}
local sum=0
local ddl=0
local rank ""
qui q memory
if `matsizeini'>r(matsize) & 10000 > r(m_matsize){
set matsize `matsizeini'
}
else{
set matsize 10000
}
forvalues g=1/`nbgroup'{
capture matrix G_`g'=deviation_`g''*W_`g'i*deviation_`g'
if _rc!=0{
noi di in red "Error, not enought memory to allocate to the W matrix calculation
}
local sum=`sum'+G_`g'[1,1]
local ddl=`ddl'+`rank`g''
local rank "`rank'`rank`g'' "
}
local ddl=`ddl'-`nbmodagpe'+`nbit'-2
set matsize `matsizeini'
}
local R1m_group "`list'"
local R1m_Rank `rank'
local R1M_ddl=`ddl'
local R1M_stat=`sum'
local R1M_p=`=1-chi2(`ddl', `sum')'
/*************************************/
/* tests Si */
/*************************************/
qui{
local countbase=0
forvalues ii=1/`nbit'{
forvalues j=1/`nbmoda_`ii''{
local countbase=`countbase'+1
}
}
forvalues s=1/`nbscore'{
local countbase=`countbase'+1
}
matrix DGibase=J(6,`countbase',.)
use `bddini'_se_g, replace
local countbase=0
forvalues ii=1/`nbit'{
tab it`ii', gen(i`ii'c)
forvalues j=1/`nbmoda_`ii''{
local countbase=`countbase'+1
matrix DGibase[2,`countbase']=`ii'
matrix DGibase[3,`countbase']=`j'
su i`ii'c`j'
matrix DGibase[4,`countbase']=r(sum)
}
drop i`ii'c*
}
tab score, gen(s)
forvalues s=1/`nbscore'{
capture su s`s'
if _rc!=0{
gen s`s'=0
}
}
forvalues s=1/`nbscore'{
local countbase=`countbase'+1
matrix DGibase[2,`countbase']=`s'
su s`s'
matrix DGibase[4,`countbase']=r(sum)
}
use `matU'_2, replace
local countbase=0
forvalues ii=1/`nbit'{
forvalues j=1/`nbmoda_`ii''{
local countbase=`countbase'+1
su P if i`ii'c`j'==1
matrix DGibase[5,`countbase']=r(sum)*`N'
matrix DGibase[6,`countbase']=(DGibase[4,`countbase']-DGibase[5,`countbase'])/sqrt(`N')
}
}
forvalues s=1/`nbscore'{
local countbase=`countbase'+1
su P if score`s'==1
matrix DGibase[5,`countbase']=r(sum)*`N'
matrix DGibase[6,`countbase']=(DGibase[4,`countbase']-DGibase[5,`countbase'])/sqrt(`N')
}
local colbase=`countbase'
}
forvalues i=1/`nbit'{
qui{
use `matU'_2, replace
gen g=.
local liminf1=-1
forvalues g=1/`nbgroup'{
replace g=`g' if score<=`limsup`g'' & score>`liminf`g''
}
forvalues g=1/`nbgroup'{
forvalues j=1/`nbmoda_`i''{
gen testi`i'c`j'g`g'=0
replace testi`i'c`j'g`g'=i`i'c`j' if g==`g'
}
}
local count=1
forvalues ii=1/`nbit'{
forvalues j=1/`nbmoda_`ii''{
rename i`ii'c`j' v`count'b
local count=`count'+1
}
}
forvalues s=1/`=`nbscore'+1'{
capture rename score`s' v`count'b
if _rc ==0{
local count=`count'+1
}
}
forvalues g=1/`nbgroup'{
forvalues j=1/`nbmoda_`i''{
rename testi`i'c`j'g`g' v`count'b
local count=`count'+1
}
}
local nbarrive`i'=`count'-1
sort `sort'
save `matU'_2i`i', replace
use `matU'_1, replace
gen g=.
local liminf1=-1
forvalues g=1/`nbgroup'{
replace g=`g' if score<=`limsup`g'' & score>`liminf`g''
}
forvalues g=1/`nbgroup'{
forvalues j=1/`nbmoda_`i''{
gen testi`i'c`j'g`g'=0
replace testi`i'c`j'g`g'=i`i'c`j' if g==`g'
}
}
local count=1
forvalues ii=1/`nbit'{
forvalues j=1/`nbmoda_`ii''{
rename i`ii'c`j' v`count'
local count=`count'+1
}
}
forvalues s=1/`=`nbscore'+1'{
capture rename score`s' v`count'
if _rc ==0{
local count=`count'+1
}
}
forvalues g=1/`nbgroup'{
forvalues j=1/`nbmoda_`i''{
rename testi`i'c`j'g`g' v`count'
local count=`count'+1
}
}
sort `sort'
save `matU'_1i`i', replace
merge `sort' using `matU'_2i`i'
drop _merge
if `nbarrive`i''>r(matsize) & `nbarrive`i''<=r(m_matsize){
set matsize `nbarrive`i''
}
if `nbarrive`i''>r(m_matsize){
di in red "Error, not enought memory to allocate to the W matrix calculation"
}
matrix i`i'W=J(`nbarrive`i'',`nbarrive`i'',.)
forvalues ii=1/`nbarrive`i''{
forvalues j=1/`nbarrive`i''{
gen prod=v`ii'*v`j'b
qui su prod
matrix i`i'W[`ii',`j']=r(sum)
drop prod
}
}
capture matrix i`i'Wi=invsym(i`i'W)
if _rc!=0{
capture matrix i`i'Wi=inv(i`i'W)
}
if _rc!=0{
di in red "Error while computing the Wg matrix"
}
mata: A=J(1,1,.)
mata: A[1,1]=rank(st_matrix("i`i'W"))
mata: st_matrix("rank_",A)
local i`i'rank=rank_[1,1]
matrix plus`i'=J(6,`=`nbarrive`i''-`colbase'',.)
matrix DGi_`i'=DGibase,plus`i'
save `matU'_i`i', replace
local count=`colbase'
use `bddini'_se_g, replace
tab it`i', gen(i`i'c)
forvalues g=1/`nbgroup'{
forvalues j=1/`nbmoda_`i''{
local count=`count'+1
matrix DGi_`i'[2,`count']=`i'
matrix DGi_`i'[3,`count']=`j'
matrix DGi_`i'[1,`count']=`g'
su i`i'c`j' if g==`g'
matrix DGi_`i'[4,`count']=r(sum)
}
}
use `matU'_i`i', replace
forvalues t=`=`colbase'+1'/`nbarrive`i''{
su P if v`t'b==1
matrix DGi_`i'[5,`t']=r(sum)*`N'
matrix DGi_`i'[6,`t']=(DGi_`i'[4,`t']-DGi_`i'[5,`t'])/sqrt(`N')
}
set matsize `matsizeini'
matrix tDGi_`i'=DGi_`i'[6,1..`=colsof(DGi_`i')']'
matrix G=tDGi_`i''*i`i'Wi*tDGi_`i'
local sum=G[1,1]
local ddl=`i`i'rank'
local ddl=`ddl'-`nbmodagpe'+`nbit'-2
local rank "`i`i'rank'"
}
local S`i'_rank=`rank'
local S`i'_ddl=`ddl'
local S`i'_stat=`sum'
local S`i'_p=`=1-chi2(`ddl', `sum')'
}
}
/*************************************/
/* print */
/*************************************/
di in gr ""
di in gr " log likelihood : " in ye "`=e(ll)'"
local lll=e(ll)
local ccnn=e(cn)
di in gr " Condition number : " in ye "`=e(cn)'"
di in gr " Number of individuals : "in ye "`Nbid'"
di in gr " Number of items : "in ye "`nbit'"
if `nbcov'!=0{
di in gr " Number of covariates : "in ye "`nbcov'"
}
di in gr ""
di ""
if "`difficulties'"==""{
di in gr "Global tests of the fit : test R1m"
di in gr " groups : " in ye "`R1m_group'"
di in gr " Number of individuals with missing data : " in ye "`nbmissing' " in gr "(" in ye "`=round(`=`percentmissing'*100',0.01)'%" in gr ")"
di ""
di in gr "{hline 46}"
di _col(22) in gr "test value" _col(35) in gr "df" _col(40) in gr "p-value"
di in gr "{hline 46}"
di _col(3) in gr "R1m" _col(24) in ye %8.5f `R1M_stat' _col(33) in ye %4.0f `R1M_ddl' _col(39) in ye %8.5f `R1M_p'
di in gr "{hline 46}"
di in gr ""
di ""
matrix globalFit=(`R1M_stat',`R1M_ddl',`R1M_p')
matrix coln globalFit=R1m df p
matrix rown globalFit=Global
}
else{
di in gr "Tests of the fit (tests R1m and Si) not performed when difficulties are fixed"
di ""
}
di in gr "Parameters of the Latent trait distribution :"
di ""
if "`difficulties'"==""{
if `nbcov'!=0{
di in gr " Latent trait mean in the subgroup `subgroup0' : fixed to 0"
}
else{
di in gr " Overall latent trait mean : fixed to 0"
}
}
di in gr " Variance of the Latent trait : Sigma²=" in ye %8.5f (`=`=eB[1,`nbc']'')^2 in gr " (SE:"in ye %8.5f `= sqrt((2*`=eB[1,`nbc']')^2 *`=eV[`nbc',`nbc']' )' in gr ")"
local varTheta=(`=`=eB[1,`nbc']'')^2
local Varvartheta=(2*`=eB[1,`nbc']')^2 *`=eV[`nbc',`nbc']'
di ""
if `nbcov'!=0{
di in gr "Latent trait group effect :"
di ""
di in gr "{hline 101}"
if "`difficulties'"==""{
local compteur=`nbdifftot'
di _col(21) in gr "Coef." _col(31) in gr "S.E." _col(41) in gr "z" _col(45) in gr "P>|z|" _col(58) in gr "[95% C.I.]" _col(71) in gr "D.St." _col(83) in gr "SS.III" _col(91) in gr "df" _col(95) in gr "V.expl."
di in gr "{hline 101}"
}
else{
di _col(22) in gr "Coef." _col(31) in gr "S.E." _col(41) in gr "z" _col(45) in gr "P>|z|" _col(58) in gr "[95% C.I.]" _col(71) in gr "D.St." _col(83) in gr "SS.III" _col(91) in gr "df" _col(95) in gr "V.expl."
di in gr "{hline 101}"
di _col(1) in gr "_cons" _col(19) in ye %7.4f `=eB[1,1]' _col(28) in ye %7.4f `=sqrt(eV[1,1])' _col(37) in ye %5.2f `=eB[1,1]/sqrt(eV[1,1])' _col(44) in ye %6.3f `=2*(1-normal(abs(eB[1,1]/sqrt(eV[1,1]))))' _col(52) in ye %7.4f `=eB[1,1]-1.96*sqrt(eV[1,1])' _col(61) in ye %7.4f `=eB[1,1]+1.96*sqrt(eV[1,1])' _col(75) in ye "."
di in gr "{hline 101}"
local compteur=1
}
forvalues i=1/`nbcov'{
local compteur=`compteur'+1
noi{
di _col(1) in gr "``=`nbit'+`i''' :" _col(82) in ye %6.3f `sst3_`=`i'+`nbit''' _col(91) in ye %2.0f `df_`=`i'+`nbit''' _col(95) in ye %6.3f `=`VarExp_`=`i'+`nbit'''*100' "%"
di _col(4) in gr "``=`nbit'+`i''' :" in ye " `ident`i'_1'" _col(25) in ye "0" _col(34) in ye "." _col(41) in ye "." _col(49) in ye "." _col(58) in ye "." _col(67) in ye "." _col(75) in ye "."
}
forvalues k=2/`nbModCov`i''{
local compteur=`compteur'+1
local estimate_e=eB[1,`compteur']
local se_e=sqrt(`=eV[`compteur',`compteur']')
noi{
di _col(4) in gr "``=`nbit'+`i''' :" in ye " `ident`i'_`k''" _col(19) in ye %7.4f `estimate_e' _col(28) in ye %7.4f `se_e' _col(37) in ye %5.2f `=`estimate_e'/`se_e'' _col(44) in ye %6.3f `=2*(1-normal(abs(`estimate_e'/`se_e')))' _col(52) in ye %7.4f `=`estimate_e'-1.96*`se_e'' _col(61) in ye %7.4f `=`estimate_e'+1.96*`se_e'' _col(70) in ye %6.3f `ESit`=`i'+`nbit''mod`k''
}
}
di in gr "{hline 101}"
}
di ""
}
else if "`difficulties'"!=""{
di in gr "Latent trait distribution"
di in gr "{hline 85}"
di _col(23) in gr "Coef." _col(36) in gr "S.E." _col(48) in gr "z" _col(57) in gr "P>|z|" _col(66) in gr "[95% Conf. Interval]"
di in gr "{hline 85}"
di _col(1) in gr "Mu" _col(20) in ye %8.5f `=eB[1,1]' _col(32) in ye %8.5f `=sqrt(eV[1,1]) ' _col(44) in ye %5.2f `=`=eB[1,1]'/`=sqrt(eV[1,1]) '' _col(56) in ye %6.3f `=2*(1-normal(abs(`=eB[1,1]'/`=sqrt(eV[1,1]) ')))' _col(66) in ye %8.5f `=`=eB[1,1]'-1.96*`=sqrt(eV[1,1]) '' _col(78) in ye %8.5f `=`=eB[1,1]'+1.96*`=sqrt(eV[1,1]) ''
di _col(1) in gr "Sigma" _col(20) in ye %8.5f `=eB[1,`nbc']' _col(32) in ye %8.5f `=sqrt(`=eV[`nbc',`nbc']')' _col(44) in ye %5.2f `=eB[1,`nbc']/sqrt(`=eV[`nbc',`nbc']')' _col(56) in ye %6.3f `=2*(1-normal(abs(eB[1,`nbc']/sqrt(`=eV[`nbc',`nbc']'))))' _col(66) in ye %8.5f `=eB[1,`nbc']-1.96*sqrt(`=eV[`nbc',`nbc']')' _col(78) in ye %8.5f `=eB[1,`nbc']+1.96*sqrt(`=eV[`nbc',`nbc']')'
di in gr "{hline 85}"
di ""
}
if "`difficulties'"==""{
matrix itemFit=J(`nbit',3,.)
local itemFitname ""
matrix coln itemFit=R1m df p
matrix rown globalFit=Global
di in gr "Items difficulty parameters :"
di ""
di in gr "{hline 80}"
di _col(1) in gr "Item" _col(21) in gr "Coef." _col(31) in gr "S.E." _col(43) in gr "[95% C.I.]" _col(66) in gr"Si" _col(70) in gr "df" _col(74) in gr "p-value"
di in gr "{hline 80}"
local compteur=1
forvalues i=1/`nbit'{
local itemFitname "`itemFitname' ``i'' "
di _col(1) in gr "``i'' :" _col(60) in ye %8.5f `S`i'_stat' _col(65) in ye %4.0f `S`i'_ddl' _col(73) in ye %8.5f `S`i'_p'
matrix itemFit[`i',1]=`S`i'_stat'
matrix itemFit[`i',2]=`S`i'_ddl'
matrix itemFit[`i',3]=`S`i'_p'
forvalues g=1/`=`moda`i''-1'{
di _col(4) in gr "modality :" in ye "`g'" _col(19) in ye %7.4f `=eB[1,`compteur']' _col(28) in ye %7.4f `=sqrt(`=eV[`compteur',`compteur']')' _col(37) in ye %7.4f `=eB[1,`compteur']-1.96*sqrt(`=eV[`compteur',`compteur']')' _col(46) in ye %7.4f `=eB[1,`compteur']+1.96*sqrt(`=eV[`compteur',`compteur']')'
local compteur=`compteur'+1
}
di in gr "{hline 80}"
}
matrix rown itemFit=`itemFitname'
di ""
}
else{
di in gr "Items difficulty parameters : fixed for the analysis"
}
/*************************************/
/* ereturn */
/*************************************/
if "`difficulties'"==""{
matrix Varbeta=e(V)
local Varsigma=Varbeta[`=colsof(Varbeta)',`=colsof(Varbeta)']
matrix Varbeta=Varbeta[1..`nbit',1..`nbit']
matrix beta=e(b)
local sigma=beta[1,`=colsof(beta)']
matrix beta=beta[1,1..`nbit']
if `nbcov'!=0{
matrix Vartheta=e(b)
matrix Vartheta=beta[`nbit'..`=colsof(Vartheta)-1',`nbit'..`=colsof(Vartheta)-1']
matrix theta=e(b)
matrix theta=beta[1,`nbit'..`=colsof(theta)-1']
}
else{
local theta=0
local Vartheta=0
}
}
else{
if `nbcov'!=0{
matrix Vartheta=e(V)
local Varsigma=Vartheta[`=colsof(Vartheta)',`=colsof(Vartheta)']
matrix Vartheta=Vartheta[1..`=colsof(Vartheta)-1',1..`=colsof(Vartheta)-1']
matrix theta=e(b)
local sigma=theta[1,`=colsof(theta)']
matrix theta=theta[1,1..`=colsof(theta)-1']
}
else{
matrix theta=e(b)
local sigma=theta[1,`=colsof(theta)']
matrix Vtheta=e(b)
local Varsigma=Vtheta[`=colsof(Vtheta)',`=colsof(Vtheta)']
local theta=theta[1,1]
local vartheta=Vtheta[1,1]
}
}
ereturn clear
ereturn scalar lll=`lll'
ereturn scalar cn=`ccnn'
ereturn scalar N=`Nbid'
ereturn scalar Nit=`nbit'
if `nbcov'!=0{
ereturn scalar Ncov=`nbcov'
}
if "`difficulties'"==""{
ereturn matrix Varbeta=Varbeta
ereturn matrix beta=beta
ereturn scalar sigma=`sigma'
ereturn scalar Varsigma=`Varsigma'
if `nbcov'!=0{
ereturn matrix Vartheta=Vartheta
ereturn matrix theta=theta
}
else{
matrix Vartheta=(0)
matrix coleq Vartheta=estimates
matrix coln Vartheta=mu
matrix roweq Vartheta=estimates
matrix rown Vartheta=mu
matrix theta=(0)
matrix coleq theta=estimates
matrix coln theta=mu
matrix rown theta=y1
ereturn matrix Vartheta=Vartheta
ereturn matrix theta=theta
}
ereturn matrix itemFit=itemFit
ereturn matrix globalFit=globalFit
}
else{
if `nbcov'!=0{
ereturn scalar sigma=`sigma'
ereturn scalar Varsigma=`Varsigma'
ereturn matrix Vartheta=Vartheta
ereturn matrix theta=theta
}
else{
ereturn scalar sigma=`sigma'
ereturn scalar Varsigma=`Varvartheta'
matrix Vartheta=(`vartheta')
matrix coleq Vartheta=estimates
matrix coln Vartheta=mu
matrix roweq Vartheta=estimates
matrix rown Vartheta=mu
ereturn matrix Vartheta=Vartheta
matrix theta=(`theta')
matrix coleq theta=estimates
matrix coln theta=mu
matrix rown theta=y1
ereturn matrix theta=theta
}
}
use `bddini', replace
restore
end