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.
535 lines
20 KiB
Plaintext
535 lines
20 KiB
Plaintext
************************************************************************************************************
|
|
* GEEkel2d: GEE for estimation of unidimensional or 2-dimensional Latent Trait models (Kelderman and Rijkes 1994)
|
|
* Version 4 : June 8, 2004
|
|
*
|
|
* Historic:
|
|
* Version 1 (2003-06-23): Jean-Benoit Hardouin
|
|
* Version 2 (2003-08-13): Jean-Benoit Hardouin
|
|
* version 3 (2003-11-06): Jean-Benoit Hardouin
|
|
*
|
|
* Use the ghquadm program (findit ghquadm)
|
|
*
|
|
* Jean-benoit Hardouin, Regional Health Observatory of Orléans - France
|
|
* jean-benoit.hardouin@neuf.fr
|
|
*
|
|
* News about this program : http://anaqol.free.fr
|
|
* FreeIRT Project : http://freeirt.free.fr
|
|
*
|
|
* Copyright 2003, 2004 Jean-Benoit Hardouin
|
|
*
|
|
* This program is free software; you can redistribute it and/or modify
|
|
* it under the terms of the GNU General Public License as published by
|
|
* the Free Software Foundation; either version 2 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with this program; if not, write to the Free Software
|
|
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
|
*
|
|
*************************************************************************************************************
|
|
|
|
program define GEEkel2d ,rclass
|
|
version 7.0
|
|
syntax varlist(min=2 numeric) [, coef(string) novar ll nbit(integer 30) critconv(real 1e-15) quad(integer 12)]
|
|
|
|
preserve
|
|
local nbitems: word count `varlist'
|
|
tokenize `varlist'
|
|
|
|
qui count
|
|
local N=r(N)
|
|
|
|
forvalues i=1/`nbitems' {
|
|
qui drop if ``i''==.
|
|
}
|
|
|
|
qui count
|
|
local Naf=r(N)
|
|
|
|
di _col(3) in green "Initial step (N=" in yellow `Naf' in green ")"
|
|
di _col(3) in yellow `=`N'-`Naf'' in green " observations are not used for missing values"
|
|
di
|
|
|
|
qui count
|
|
local N=r(N)
|
|
tempname B Q
|
|
if "`coef'"!="" {
|
|
matrix `B'=`coef'
|
|
}
|
|
else {
|
|
matrix `B'=J(`nbitems',1,1)
|
|
}
|
|
|
|
scalar `Q'=colsof(`B')
|
|
|
|
/* CALCUL INITIAUX DES PARAMETRES DELTA ET SIGMA ET DE LA MATRICE BETA*/
|
|
|
|
local sigmath11=.25
|
|
|
|
if `Q'==2 {
|
|
local sigmath22=.25
|
|
local sigmath12=0.125
|
|
}
|
|
|
|
tempname beta
|
|
matrix `beta'=J(`nbitems'+`Q'*(`Q'+1)/2,1,0)
|
|
forvalues j=1/`nbitems' {
|
|
qui count if ``j''==1
|
|
local s``j''=r(N)
|
|
local delta``j''=-log(`s``j'''/(`N'-`s``j'''))
|
|
matrix `beta'[`j',1]=`delta``j'''
|
|
}
|
|
|
|
|
|
forvalues j=1/`nbitems' {
|
|
forvalues l=1/`nbitems' {
|
|
quiet count if ``j''==1&``l''==1
|
|
local s``j''``l''=r(N)
|
|
}
|
|
}
|
|
|
|
if "`var'"==""{
|
|
forvalues i=1/`nbitems' {
|
|
forvalues j=1/`nbitems' {
|
|
forvalues k=1/`nbitems' {
|
|
quiet count if ``i''==1&``j''==1&``k''==1
|
|
local s``i''``j''``k''=r(N)
|
|
}
|
|
}
|
|
}
|
|
|
|
forvalues i=1/`nbitems' {
|
|
forvalues j=1/`nbitems' {
|
|
forvalues k=1/`nbitems' {
|
|
forvalues l=1/`nbitems' {
|
|
quiet count if ``i''==1&``j''==1&``k''==1&``l''==1
|
|
local s``i''``j''``k''``l''=r(N)
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
local l=`nbitems'+1
|
|
matrix `beta'[`l',1]=`sigmath11'
|
|
|
|
if `Q'==2 {
|
|
local l=`nbitems'+2
|
|
matrix `beta'[`l',1]=`sigmath22'
|
|
|
|
local l=`nbitems'+3
|
|
matrix `beta'[`l',1]=`sigmath12'
|
|
}
|
|
|
|
|
|
|
|
tempname variat V11 V12 V21 V22 D11 D12 D21 D22 V D
|
|
matrix `variat'=(1)
|
|
|
|
local compteur=0
|
|
local conv=1
|
|
|
|
|
|
/*********ITERATIONS******************/
|
|
|
|
|
|
|
|
while (`variat'[1,1]>`critconv'&`compteur'<=`nbit'&`conv'==1) {
|
|
|
|
if `compteur'==0{
|
|
di in green _col(3) "First iteration"
|
|
}
|
|
else {
|
|
di in green _col(3) "iteration:" in yellow _col(14) "`compteur'" in green _col(25) "Convergence index:" in yellow _col(44) %10.7e "`macrovariat'"
|
|
}
|
|
|
|
local compteur=`compteur'+1
|
|
|
|
|
|
forvalues j=1/`nbitems' {
|
|
|
|
/* CALCUL DES DERIVEES 1 A 6 POUR CHAQUE ITEM*/
|
|
|
|
local l1``j''=1/(1+exp(`beta'[`j',1]))
|
|
local l2``j''=exp(`beta'[`j',1])/(1+exp(`beta'[`j',1]))^2
|
|
local l3``j''=exp(`beta'[`j',1])*(exp(`beta'[`j',1])-1)/(1+exp(`beta'[`j',1]))^3
|
|
local l4``j''=exp(`beta'[`j',1])*(exp(2*`beta'[`j',1])-4*exp(`beta'[`j',1])+1)/(1+exp(`beta'[`j',1]))^4
|
|
local l5``j''=exp(`beta'[`j',1])*(exp(3*`beta'[`j',1])-11*exp(2*`beta'[`j',1])+11*exp(`beta'[`j',1])-1)/(1+exp(`beta'[`j',1]))^5
|
|
local l6``j''=exp(`beta'[`j',1])*(exp(4*`beta'[`j',1])-26*exp(3*`beta'[`j',1])+66*exp(2*`beta'[`j',1])-26*exp(`beta'[`j',1])+1)/(1+exp(`beta'[`j',1]))^6
|
|
|
|
if `Q'==2 {
|
|
local H2``j''=`B'[`j',1]^2*`sigmath11'+`B'[`j',1]*`B'[`j',2]*`sigmath12'+`B'[`j',2]^2*`sigmath22'
|
|
local H4``j''=3*`B'[`j',1]^4*`sigmath11'^2+12*`B'[`j',1]^3*`B'[`j',2]*`sigmath11'*`sigmath12'+6*`B'[`j',1]^2*`B'[`j',2]^2*(`sigmath11'*`sigmath22'+2*`sigmath12'^2)+12*`B'[`j',1]*`B'[`j',2]^3*`sigmath22'*`sigmath12'+3*`B'[`j',2]^4*`sigmath22'^2
|
|
}
|
|
if `Q'==1 {
|
|
local H2``j''=`B'[`j',1]^2*`sigmath11'
|
|
local H4``j''=3*`B'[`j',1]^4*`sigmath11'^2
|
|
}
|
|
|
|
|
|
/* CALCUL DES MOMENTS D'ORDRE 1 ET 2 ET DE LA MATRICE V11*/
|
|
|
|
local mu``j''=`l1``j'''+`H2``j'''/2*`l3``j'''+`H4``j'''/24*`l5``j'''
|
|
local sigma``j''``j''=`l2``j'''+`H2``j'''/2*(`l3``j'''-2*`l1``j'''*`l3``j''')+`H4``j'''/24*(`l5``j'''-2*(`l3``j''')^2-2*`l1``j'''*`l5``j''')
|
|
}
|
|
|
|
|
|
matrix `V11'=J(`nbitems',`nbitems',0)
|
|
|
|
forvalue j=1/`nbitems' {
|
|
forvalue l=`j'/`nbitems' {
|
|
|
|
if `j'!=`l' {
|
|
if `Q'==2 {
|
|
local H2``j''``l''=`B'[`j',1]*`B'[`l',1]*`sigmath11'+(`B'[`j',1]*`B'[`l',2]+`B'[`j',2]*`B'[`l',1])*`sigmath12'+`B'[`j',2]*`B'[`l',2]*`sigmath22'
|
|
local H4``j''1``l''3=3*`B'[`j',1]*`B'[`l',1]^3*`sigmath11'^2+3*(3*`B'[`j',1]*`B'[`l',1]^2*`B'[`l',2]+`B'[`j',2]*`B'[`l',1]^3)*`sigmath11'*`sigmath12'+(3*`B'[`j',1]*`B'[`l',1]*`B'[`l',2]^2+3*`B'[`j',2]*`B'[`l',1]^2*`B'[`l',2])*(`sigmath11'*`sigmath22'+2*`sigmath12'^2)+3*(`B'[`j',1]*`B'[`l',2]^3+3*`B'[`j',2]*`B'[`l',1]*`B'[`l',2]^2)*`sigmath22'*`sigmath12'+3*`B'[`j',2]*`B'[`l',2]^3*`sigmath22'^2
|
|
local H4``j''2``l''2=3*`B'[`j',1]^2*`B'[`l',1]^2*`sigmath11'^2+6*(`B'[`j',1]^2*`B'[`l',1]*`B'[`l',2]+`B'[`j',1]*`B'[`j',2]*`B'[`l',1]^2)*`sigmath11'*`sigmath12'+(`B'[`j',1]^2*`B'[`l',2]^2+4*`B'[`j',1]*`B'[`j',2]*`B'[`l',1]*`B'[`l',2]+`B'[`j',2]^2*`B'[`l',1]^2)*(`sigmath11'*`sigmath22'+2*`sigmath12'^2)+6*(`B'[`j',1]*`B'[`j',2]*`B'[`l',2]^2+`B'[`j',2]^2*`B'[`l',1]*`B'[`l',2])*`sigmath22'*`sigmath12'+3*`B'[`j',2]^2*`B'[`l',2]^2*`sigmath22'^2
|
|
local H4``j''3``l''1=3*`B'[`l',1]*`B'[`j',1]^3*`sigmath11'^2+3*(3*`B'[`l',1]*`B'[`j',1]^2*`B'[`j',2]+`B'[`l',2]*`B'[`j',1]^3)*`sigmath11'*`sigmath12'+(3*`B'[`l',1]*`B'[`j',1]*`B'[`j',2]^2+3*`B'[`l',2]*`B'[`j',1]^2*`B'[`j',2])*(`sigmath11'*`sigmath22'+2*`sigmath12'^2)+3*(`B'[`l',1]*`B'[`j',2]^3+3*`B'[`l',2]*`B'[`j',1]*`B'[`j',2]^2)*`sigmath22'*`sigmath12'+3*`B'[`l',2]*`B'[`j',2]^3*`sigmath22'^2
|
|
}
|
|
if `Q'==1 {
|
|
local H2``j''``l''=`B'[`j',1]*`B'[`l',1]*`sigmath11'
|
|
local H4``j''1``l''3=3*`B'[`j',1]*`B'[`l',1]^3*`sigmath11'^2
|
|
local H4``j''2``l''2=3*`B'[`j',1]^2*`B'[`l',1]^2*`sigmath11'^2
|
|
local H4``j''3``l''1=3*`B'[`l',1]*`B'[`j',1]^3*`sigmath11'^2
|
|
}
|
|
|
|
|
|
local sigma``j''``l''=`H2``j''``l'''*(`l2``j'''*`l2``l''')+`H4``j''1``l''3'/6*`l2``j'''*`l4``l'''+`H4``j''3``l''1'/6*`l4``j'''*`l2``l'''+(`H4``j''2``l''2'-`H2``j'''*`H2``l''')/4*`l3``j'''*`l3``l'''
|
|
|
|
|
|
}
|
|
|
|
matrix `V11'[`j',`l']=`sigma``j''``l'''
|
|
matrix `V11'[`l',`j']=`sigma``j''``l'''
|
|
}
|
|
}
|
|
|
|
/* DEFINITION DE LA MATRICE COMPOCARRE*/
|
|
|
|
tempname compocarre m
|
|
local carre=`nbitems'*(`nbitems'-1)/2
|
|
matrix `compocarre'=J(2,`carre',0)
|
|
|
|
scalar `m'=0
|
|
forvalue j=1/`nbitems' {
|
|
local lbis=`j'+1
|
|
forvalue l=`lbis'/`nbitems' {
|
|
scalar `m'=`m'+1
|
|
matrix `compocarre'[1,`m']=`j'
|
|
matrix `compocarre'[2,`m']=`l'
|
|
}
|
|
}
|
|
|
|
/* CALCUL DE LA MATRICE V22*/
|
|
|
|
matrix `V22'=J(`carre',`carre',0)
|
|
forvalue k=1/`carre' {
|
|
local j=`compocarre'[1,`k']
|
|
local l=`compocarre'[2,`k']
|
|
|
|
matrix `V22'[`k',`k']=(1-2*`mu``j''')*(1-2*`mu``l''')*`sigma``j''``l'''+`sigma``j''``j'''*`sigma``l''``l'''-`sigma``j''``l'''^2
|
|
}
|
|
|
|
/* CALCUL DES MATRICES V12, V21 ET V*/
|
|
|
|
matrix `V12'=J(`nbitems',`carre',0)
|
|
forvalue k=1/`carre' {
|
|
local j=`compocarre'[1,`k']
|
|
local l=`compocarre'[2,`k']
|
|
|
|
matrix `V12'[`j',`k']=(1-2*`mu``j''')*`sigma``j''``l'''
|
|
matrix `V12'[`l',`k']=(1-2*`mu``l''')*`sigma``j''``l'''
|
|
}
|
|
matrix `V21'=`V12' '
|
|
matrix `V'=(`V11',`V12' \ `V21',`V22')
|
|
|
|
|
|
/*CALCUL DES MATRICES D11*/
|
|
|
|
matrix `D11'=J(`nbitems',`nbitems',0)
|
|
matrix `D12'=J(`nbitems',`Q'*(`Q'+1)/2,0)
|
|
|
|
forvalue j=1/`nbitems' {
|
|
matrix `D11'[`j',`j']=-`l2``j'''-`H2``j'''/2*`l4``j'''-`H4``j'''/24*`l6``j'''
|
|
|
|
if `Q'==2 {
|
|
matrix `D12'[`j',1]=`B'[`j',1]^2*`l3``j'''/2+(`B'[`j',1]^4*`sigmath11'+2*`B'[`j',1]^3*`B'[`j',2]*`sigmath12'+`B'[`j',1]^2*`B'[`j',2]^2*`sigmath22')/4*`l5``j'''
|
|
matrix `D12'[`j',2]=`B'[`j',2]^2*`l3``j'''/2+(`B'[`j',2]^4*`sigmath22'+2*`B'[`j',2]^3*`B'[`j',1]*`sigmath12'+`B'[`j',2]^2*`B'[`j',1]^2*`sigmath11')/4*`l5``j'''
|
|
matrix `D12'[`j',3]=`B'[`j',1]*`B'[`j',2]*`l3``j'''+(`B'[`j',1]^3*`B'[`j',2]*`sigmath11'+2*`B'[`j',1]^2*`B'[`j',2]^2*`sigmath12'+`B'[`j',1]*`B'[`j',2]^3*`sigmath22')/2*`l5``j'''
|
|
}
|
|
if `Q'==1 {
|
|
matrix `D12'[`j',1]=`B'[`j',1]^2*`l3``j'''/2+`B'[`j',1]^4*`sigmath11'/4*`l5``j'''
|
|
}
|
|
}
|
|
|
|
/*CALCUL DES MATRICES D21, D22 et D*/
|
|
|
|
|
|
matrix `D21'=J(`carre',`nbitems',0)
|
|
matrix `D22'=J(`carre',`Q'*(`Q'+1)/2,0)
|
|
|
|
forvalue k=1/`carre' {
|
|
local j=`compocarre'[1,`k']
|
|
local l=`compocarre'[2,`k']
|
|
|
|
matrix `D21'[`k',`j']=-`H2``j''``l'''*`l3``j'''*`l2``l'''-`H4``j''1``l''3'/6*`l3``j'''*`l4``l'''-`H4``j''3``l''1'/6*`l5``j'''*`l2``l'''-(`H4``j''2``l''2'-`H2``j'''*`H2``l''')/4*`l4``j'''*`l3``l'''
|
|
matrix `D21'[`k',`l']=-`H2``j''``l'''*`l2``j'''*`l3``l'''-`H4``j''3``l''1'/6*`l4``j'''*`l3``l'''-`H4``j''1``l''3'/6*`l2``j'''*`l5``l'''-(`H4``j''2``l''2'-`H2``j'''*`H2``l''')/4*`l3``j'''*`l4``l'''
|
|
|
|
if `Q'==2 {
|
|
scalar tmp1=`B'[`j',1]*`B'[`l',1]*`l2``j'''*`l2``l'''+(2*`B'[`j',1]*`B'[`l',1]^3*`sigmath11'+(3*`B'[`j',1]*`B'[`l',1]^2*`B'[`l',2]+`B'[`j',2]*`B'[`l',1]^3)*`sigmath12'+(`B'[`j',1]*`B'[`l',1]*`B'[`l',2]^2+`B'[`j',2]*`B'[`l',1]^2*`B'[`l',2])*`sigmath22')/2*`l2``j'''*`l4``l'''
|
|
scalar tmp2=(2*`B'[`j',1]^3*`B'[`l',1]*`sigmath11'+(3*`B'[`j',1]^2*`B'[`j',2]*`B'[`l',1]+`B'[`j',1]^3*`B'[`l',2])*`sigmath12'+(`B'[`j',1]*`B'[`j',2]^2*`B'[`l',1]+`B'[`j',1]^2*`B'[`j',2]*`B'[`l',2])*`sigmath22')/2*`l4``j'''*`l2``l'''
|
|
scalar tmp3=(`B'[`j',1]^2*`B'[`l',1]^2*`sigmath11'+(`B'[`j',1]^2*`B'[`l',1]*`B'[`l',2]+`B'[`j',1]*`B'[`j',2]*`B'[`l',1]^2)*`sigmath12'+`B'[`j',1]*`B'[`j',2]*`B'[`l',1]*`B'[`l',2]*`sigmath22')*`l3``j'''*`l3``l'''
|
|
|
|
matrix `D22'[`k',1]=tmp1+tmp2+tmp3
|
|
|
|
scalar tmp1=`B'[`j',2]*`B'[`l',2]*`l2``j'''*`l2``l'''+(2*`B'[`j',2]*`B'[`l',2]^3*`sigmath22'+(3*`B'[`j',2]*`B'[`l',2]^2*`B'[`l',1]+`B'[`j',1]*`B'[`l',2]^3)*`sigmath12'+(`B'[`j',2]*`B'[`l',2]*`B'[`l',1]^2+`B'[`j',1]*`B'[`l',2]^2*`B'[`l',1])*`sigmath11')/2*`l2``j'''*`l4``l'''
|
|
scalar tmp2=(2*`B'[`j',2]^3*`B'[`l',2]*`sigmath22'+(3*`B'[`j',2]^2*`B'[`j',1]*`B'[`l',2]+`B'[`j',2]^3*`B'[`l',1])*`sigmath12'+(`B'[`j',2]*`B'[`j',1]^2*`B'[`l',2]+`B'[`j',2]^2*`B'[`j',1]*`B'[`l',1])*`sigmath11')/2*`l4``j'''*`l2``l'''
|
|
scalar tmp3=(`B'[`j',1]^2*`B'[`l',1]^2*`sigmath22'+(`B'[`j',1]^2*`B'[`l',1]*`B'[`l',2]+`B'[`j',1]*`B'[`j',2]*`B'[`l',1]^2)*`sigmath12'+`B'[`j',1]*`B'[`j',2]*`B'[`l',1]*`B'[`l',2]*`sigmath11')*`l3``j'''*`l3``l'''
|
|
|
|
matrix `D22'[`k',2]=tmp1+tmp2+tmp3
|
|
|
|
scalar tmp1=(`B'[`j',1]*`B'[`l',2]+`B'[`j',2]*`B'[`l',1])*`l2``j'''*`l2``l'''+((3*`B'[`j',1]*`B'[`l',1]^2*`B'[`l',2]+`B'[`j',2]*`B'[`l',1]^3)*`sigmath11'+4*(`B'[`j',1]*`B'[`l',1]*`B'[`l',2]^2+`B'[`j',2]*`B'[`l',1]^2*`B'[`l',2])*`sigmath12'+(`B'[`j',1]*`B'[`l',2]^3+3*`B'[`j',2]*`B'[`l',1]*`B'[`l',2]^2)*`sigmath22')/2*`l2``j'''*`l4``l'''
|
|
scalar tmp2=((3*`B'[`j',1]^2*`B'[`j',2]*`B'[`l',1]+`B'[`j',1]^3*`B'[`l',2])*`sigmath11'+4*(`B'[`j',1]*`B'[`j',2]^2*`B'[`l',1]+`B'[`j',1]^2*`B'[`j',2]*`B'[`l',2])*`sigmath12'+(`B'[`j',2]^3*`B'[`l',1]+3*`B'[`j',1]*`B'[`j',2]^2*`B'[`l',2])*`sigmath22')/2*`l4``j'''*`l2``l'''
|
|
scalar tmp3=((`B'[`j',1]^2*`B'[`l',1]*`B'[`l',2]+`B'[`j',1]*`B'[`j',2]*`B'[`l',2]^2)*`sigmath11'+(`B'[`j',1]^2*`B'[`l',2]^2+2*`B'[`j',1]*`B'[`j',2]*`B'[`l',1]*`B'[`l',2]+`B'[`j',2]^2*`B'[`l',1]^2)*`sigmath12'+(`B'[`j',1]*`B'[`j',2]*`B'[`l',2]^2+`B'[`j',2]^2*`B'[`l',1]*`B'[`l',2])*`sigmath22')*`l3``j'''*`l3``l'''
|
|
|
|
matrix `D22'[`k',3]=tmp1+tmp2+tmp3
|
|
}
|
|
if `Q'==1 {
|
|
scalar tmp1=`B'[`j',1]*`B'[`l',1]*`l2``j'''*`l2``l'''+(2*`B'[`j',1]*`B'[`l',1]^3*`sigmath11')/2*`l2``j'''*`l4``l'''
|
|
scalar tmp2=(2*`B'[`j',1]^3*`B'[`l',1]*`sigmath11')/2*`l4``j'''*`l2``l'''
|
|
scalar tmp3=(`B'[`j',1]^2*`B'[`l',1]^2*`sigmath11')*`l3``j'''*`l3``l'''
|
|
|
|
matrix `D22'[`k',1]=tmp1+tmp2+tmp3
|
|
}
|
|
}
|
|
|
|
matrix `D'=(`D11',`D12' \ `D21',`D22')
|
|
|
|
/*CALCUL DE LA MATRICE CHSI*/
|
|
|
|
tempname chsi
|
|
matrix `chsi'=J(`nbitems'+`carre',1,0)
|
|
|
|
forvalue j=1/`nbitems' {
|
|
matrix `chsi'[`j',1]=(`s``j'''-`N'*`mu``j''')/`N'
|
|
}
|
|
|
|
|
|
forvalue k=1/`carre' {
|
|
local j=`compocarre'[1,`k']
|
|
local l=`compocarre'[2,`k']
|
|
local tmp=`nbitems'+`k'
|
|
|
|
matrix `chsi'[`tmp',1]=(`s``j''``l'''-`s``j'''*`mu``l'''-`s``l'''*`mu``j'''+`N'*`mu``j'''*`mu``l'''-`N'*`sigma``j''``l''')/`N'
|
|
}
|
|
|
|
/*CALCUL DE L'ETAPE k*/
|
|
tempname betaold
|
|
matrix `betaold'=`beta'
|
|
matrix `beta'=`betaold'+inv(`D''*inv(`V')*`D')*`D''*inv(`V')*`chsi'
|
|
|
|
local l=`nbitems'+1
|
|
local sigmath11=`beta'[`l',1]
|
|
local l=`nbitems'+2
|
|
local sigmath22=`beta'[`l',1]
|
|
local l=`nbitems'+3
|
|
local sigmath12=`beta'[`l',1]
|
|
|
|
tempname epsilon variatold
|
|
scalar `variatold'=`variat'[1,1]
|
|
matrix `epsilon'=`betaold'-`beta'
|
|
|
|
matrix `variat'=(`epsilon''*`epsilon')
|
|
|
|
if `variat'[1,1]>`variatold' {
|
|
matrix `beta'=`betaold'
|
|
local l=`nbitems'+1
|
|
local sigmath11=`beta'[`l',1]
|
|
local l=`nbitems'+2
|
|
local sigmath22=`beta'[`l',1]
|
|
local l=`nbitems'+3
|
|
local sigmath12=`beta'[`l',1]
|
|
local conv=0
|
|
}
|
|
else {
|
|
local macrovariat=`variat'[1,1]
|
|
}
|
|
|
|
}
|
|
|
|
/*************************CALCUL des STANDARDS ERRORS DES PARAMETRES *********************/
|
|
|
|
if "`var'"==""{
|
|
tempname xicarreA xicarreB xicarreC xicarre
|
|
matrix `xicarreA'=J(`nbitems',`nbitems',0)
|
|
matrix `xicarreB'=J(`nbitems',`carre',0)
|
|
matrix `xicarreC'=J(`carre',`carre',0)
|
|
forvalues i=1/`nbitems' {
|
|
forvalues j=1/`nbitems' {
|
|
matrix `xicarreA'[`i',`j']=`s``i''``j'''-`s``i'''*`mu``j'''-`s``j'''*`mu``i'''+`N'*`mu``i'''*`mu``j'''
|
|
}
|
|
}
|
|
|
|
forvalues i=1/`nbitems' {
|
|
forvalues col=1/`carre' {
|
|
local j=`compocarre'[1,`col']
|
|
local k=`compocarre'[2,`col']
|
|
matrix `xicarreB'[`i',`col']=`s``i''``j''``k'''-`mu``i'''*`s``j''``k'''-`mu``j'''*`s``i''``k'''-`mu``k'''*`s``i''``j'''+`mu``i'''*`mu``j'''*`s``k'''+`mu``i'''*`mu``k'''*`s``j'''+`mu``j'''*`mu``k'''*`s``i'''-`N'*`mu``i'''*`mu``j'''*`mu``k'''-`sigma``j''``k'''*`s``i'''+`N'*`mu``i'''*`sigma``j''``k'''
|
|
}
|
|
}
|
|
|
|
forvalues row=1/`carre' {
|
|
forvalues col=1/`carre' {
|
|
local i=`compocarre'[1,`row']
|
|
local j=`compocarre'[2,`row']
|
|
local k=`compocarre'[1,`col']
|
|
local l=`compocarre'[2,`col']
|
|
matrix `xicarreC'[`row',`col']=`s``i''``j''``k''``l'''-`mu``i'''*`s``j''``k''``l'''-`mu``j'''*`s``i''``k''``l'''-`mu``k'''*`s``i''``j''``l'''-`mu``l'''*`s``i''``j''``k'''+`mu``i'''*`mu``j'''*`s``k''``l'''+`mu``i'''*`mu``k'''*`s``j''``l'''+`mu``i'''*`mu``l'''*`s``j''``k'''+`mu``j'''*`mu``k'''*`s``i''``l'''+`mu``j'''*`mu``l'''*`s``i''``k'''+`mu``k'''*`mu``l'''*`s``i''``j'''-`mu``i'''*`mu``j'''*`mu``k'''*`s``l'''-`mu``i'''*`mu``j'''*`mu``l'''*`s``k'''-`mu``i'''*`mu``k'''*`mu``l'''*`s``j'''-`mu``j'''*`mu``k'''*`mu``l'''*`s``i'''-`sigma``i''``j'''*`s``k''``l'''-`sigma``k''``l'''*`s``i''``j'''+`sigma``i''``j'''*`mu``k'''*`s``l'''+`sigma``i''``j'''*`mu``l'''*`s``k'''+`sigma``k''``l'''*`mu``i'''*`s``j'''+`sigma``k''``l'''*`mu``j'''*`s``i'''+`N'*`mu``i'''*`mu``j'''*`mu``k'''*`mu``l'''-`N'*`sigma``i''``j'''*`mu``k'''*`mu``l'''-`N'*`sigma``k''``l'''*`mu``i'''*`mu``j'''+`N'*`sigma``i''``j'''*`sigma``k''``l'''
|
|
}
|
|
}
|
|
|
|
|
|
matrix `xicarre'=(`xicarreA',`xicarreB' \ `xicarreB' ',`xicarreC')
|
|
|
|
tempname A1 A2 W
|
|
matrix `A1'=`D' '*inv(`V')*`D'
|
|
matrix `A2'=`D' '*inv(`V')*`xicarre'*inv(`V')*`D'
|
|
matrix `W'=1/`N'^2*inv(`A1')*`A2'*inv(`A1')
|
|
|
|
}
|
|
|
|
/*****************************DISPLAY THE RESULTS***************************************/
|
|
|
|
|
|
local compteur=`compteur'-1
|
|
di ""
|
|
di ""
|
|
if `compteur'==0 {
|
|
noi di in red _col(8) "The algorithm can not converge"
|
|
return scalar error=1
|
|
exit
|
|
}
|
|
if `variat'[1,1]<=`critconv'&`compteur'>0 {
|
|
noi di in green _col(8) "The algorithm converge at the `compteur'th iteration"
|
|
}
|
|
if `compteur'==`nbit'&`variat'[1,1]>`critconv' {
|
|
noi di in green _col(8) "The algorithm is stopped at the `compteur'th iteration"
|
|
}
|
|
if `conv'==0&`compteur'>0 {
|
|
noi di in green _col(8) "The algorithm converge no more after the `compteur'th iteration"
|
|
}
|
|
di ""
|
|
if "`var'"=="" {
|
|
noi di in green _col(30) "Parameters" in green _col(43) "Standard errors"
|
|
forvalue j=1/`nbitems' {
|
|
noi di in green _col(20) "``j'': " in yellow _col(30) %10.6f `beta'[`j',1] in yellow _col(50) %8.6f sqrt(`W'[`j',`j'])
|
|
}
|
|
di ""
|
|
noi di in green _col(20) "var1: " in yellow _col(30) %10.6f `beta'[`nbitems'+1,1] in yellow _col(50) %8.6f sqrt(`W'[`nbitems'+1,`nbitems'+1])
|
|
if `Q'==2 {
|
|
noi di in green _col(20) "var2: " in yellow _col(30) %10.6f `beta'[`nbitems'+2,1] in yellow _col(50) %8.6f sqrt(`W'[`nbitems'+2,`nbitems'+2])
|
|
scalar rho=`beta'[`nbitems'+3,1]/sqrt(`beta'[`nbitems'+1,1]*`beta'[`nbitems'+2,1])
|
|
noi di in green _col(20) "covar: " in yellow _col(30) %10.6f `beta'[`nbitems'+3,1] in yellow _col(50) %8.6f sqrt(`W'[`nbitems'+3,`nbitems'+3]) " (rho=" %5.4f rho ")"
|
|
}
|
|
}
|
|
else {
|
|
noi di in green _col(30) "Parameters"
|
|
forvalue j=1/`nbitems' {
|
|
noi di in green _col(20) "``j'': " in yellow _col(30) %10.6f `beta'[`j',1]
|
|
}
|
|
di ""
|
|
noi di in green _col(20) "var1: " in yellow _col(30) %10.6f `beta'[`nbitems'+1,1]
|
|
if `Q'==2 {
|
|
noi di in green _col(20) "var2: " in yellow _col(30) %10.6f `beta'[`nbitems'+2,1]
|
|
scalar rho=`beta'[`nbitems'+3,1]/sqrt(`beta'[`nbitems'+1,1]*`beta'[`nbitems'+2,1])
|
|
noi di in green _col(20) "covar: " in yellow _col(30) %10.6f `beta'[`nbitems'+3,1]
|
|
}
|
|
}
|
|
di ""
|
|
|
|
|
|
|
|
if "`ll'"!="" {
|
|
if `Q'==1 {
|
|
tempname noeuds poids
|
|
ghquadm `quad' `noeuds' `poids'
|
|
tempvar vrais logvrais
|
|
qui gen `vrais'=0
|
|
forvalues u=1/`quad'{
|
|
tempvar vrais`u'
|
|
qui gen `vrais`u''=1/sqrt(_pi)
|
|
forvalues j=1/`nbitems' {
|
|
tempvar Pu`u'j`j'
|
|
qui gen `Pu`u'j`j''=exp(`B'[`j',1]*sqrt(2*`beta'[`nbitems'+1,1])*`noeuds'[1,`u']-`beta'[`j',1])/(1+exp(`B'[`j',1]*sqrt(2*`beta'[`nbitems'+1,1])*`noeuds'[1,`u']-`beta'[`j',1]))
|
|
qui replace `Pu`u'j`j''=1-`Pu`u'j`j'' if ``j''==0
|
|
qui replace `vrais`u''=`vrais`u''*`Pu`u'j`j''
|
|
}
|
|
qui replace `vrais'=`vrais'+`poids'[1,`u']*`vrais`u''
|
|
}
|
|
gen `logvrais'=log(`vrais')
|
|
qui su `logvrais'
|
|
local ll=r(N)*r(mean)
|
|
noi di in green _col(20) "ll: " in yellow _col(30) %12.4f `ll'
|
|
local AIC=-2*`ll'+2*(`nbitems'+1)
|
|
noi di in green _col(20) "AIC: " in yellow _col(30) %12.4f `AIC'
|
|
}
|
|
|
|
if `Q'==2 {
|
|
|
|
tempname noeuds poids
|
|
ghquadm `quad' `noeuds' `poids'
|
|
tempvar vrais logvrais
|
|
qui gen `vrais'=0
|
|
matrix sigma=(`beta'[`nbitems'+1,1],`beta'[`nbitems'+3,1] \ `beta'[`nbitems'+3,1],`beta'[`nbitems'+2,1])
|
|
forvalues u=1/`quad'{
|
|
forvalues v=1/`quad'{
|
|
tempvar vraisu`u'v`v'
|
|
qui gen `vraisu`u'v`v''=1/_pi
|
|
forvalues j=1/`nbitems' {
|
|
tempvar Pu`u'v`v'j`j'
|
|
local A1`u'tilde=sqrt(`beta'[`nbitems'+2,1]/(2*det(sigma)))*`noeuds'[1,`u']
|
|
local A2`v'tilde=(`noeuds'[1,`v']-`beta'[`nbitems'+3,1]/sqrt(det(sigma))*`noeuds'[1,`u'])/(2*`beta'[`nbitems'+2,1])
|
|
qui gen `Pu`u'v`v'j`j''=exp(`B'[`j',1]*`A1`u'tilde'+`B'[`j',2]*`A2`v'tilde'-`beta'[`j',1])/(1+exp(`B'[`j',1]*`A1`u'tilde'+`B'[`j',2]*`A2`v'tilde'-`beta'[`j',1]))
|
|
qui replace `Pu`u'v`v'j`j''=1-`Pu`u'v`v'j`j'' if ``j''==0
|
|
qui replace `vraisu`u'v`v''=`vraisu`u'v`v''*`Pu`u'v`v'j`j''
|
|
}
|
|
qui replace `vrais'=`vrais'+`poids'[1,`u']*`poids'[1,`v']*`vraisu`u'v`v''
|
|
}
|
|
}
|
|
qui gen `logvrais'=log(`vrais')
|
|
qui su `logvrais'
|
|
local ll=r(N)*r(mean)
|
|
noi di in green _col(20) "ll: " in yellow _col(27) %12.4f `ll'
|
|
local AIC=-2*`ll'+2*(`nbitems'+3)
|
|
noi di in green _col(20) "AIC: " in yellow _col(27) %12.4f `AIC'
|
|
}
|
|
}
|
|
|
|
if "`var'"=="" {
|
|
return matrix W= `W'
|
|
}
|
|
return matrix beta= `beta'
|
|
|
|
if "`ll'"!="" {
|
|
return scalar ll= `ll'
|
|
return scalar AIC= `AIC'
|
|
}
|
|
|
|
return scalar J= `nbitems'
|
|
return scalar N= `N'
|
|
|
|
return scalar error=0
|
|
|
|
restore
|
|
|
|
|
|
|
|
end
|
|
|