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