corentinchoisy/simul_these
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

File restructure #2

Browse Source
This commit is contained in:
Corentin Choisy
2024-04-19 16:48:36 +02:00
parent ecac05b9c4
commit 0e9e01eca8
702 changed files with 272561 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

113
Modules/ado/personal/g/gammasym.ado Normal file
View File

@ -0,0 +1,113 @@
*! version 2.2 09march2007
*! Jean-Benoit Hardouin
************************************************************************************************************
* Gammasym : Symmetric gamma function
*
* Historic:
* Version 1 (2004-01-29): Jean-Benoit Hardouin
* Version 2 (2004-02-01): Jean-Benoit Hardouin
* Version 2.1 (2005-04-04) : Jean-benoit Hardouin
* Version 2.2 (2007-03-09) : Jean-benoit Hardouin /*Weights*/
*
* Jean-benoit Hardouin, University of Nantes - France
* jean-benoit.hardouin@univ-nantes.fr
*
* News about this program : http://www.anaqol.org
* FreeIRT Project : http://www.freeirt.org
*
* Copyright 2004-2007 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 gammasym,rclass
version 7
syntax anything [, Weight(string) DISplay]
local B `weight'
if "`B'"!="" {
numlist "`B'"
local B=r(numlist)
local nbw:word count `B'
local score=0
forvalues i=1/`nbw' {
local B`i':word `i' of `B'
local score=`score'+`B`i''
}
}
local nbgroups:word count `anything'
if "`B'"=="" {
forvalues i=1/`nbgroups' {
local B`i'=1
}
local nbw=`nbgroups'
local score=`nbgroups'
}
if `nbgroups'!=`nbw' {
di in red "There is not the same number of values and of weights. Please correct that."
error
}
tokenize `anything'
local list
forvalues i=1/`nbgroups' {
local tmp=``i''
local list "`list' `tmp'"
}
tokenize `list'
tempname G
matrix `G'=J(`=`score'+4',`nbgroups',0)
local poids=`score'+2
local epsilon=`score'+3
local max=`score'+4
forvalues j=1/`nbgroups' {
matrix `G'[`poids',`j']=`B`j''
matrix `G'[`epsilon',`j']=exp(-``j'')^`B`j''
if `j'!=1 {
matrix `G'[`max',`j']=`G'[`max',`=`j'-1']+`B`j''
}
matrix `G'[1,`j']=1
}
matrix `G'[`=`B1'+1',1]=exp(-`1')^`B1'
forvalues j=2/`nbgroups' {
local w=`B`j''
forvalues s=2/`=`score'+1' {
matrix `G'[`s',`j']=`G'[`s',`j']+`G'[`s',`=`j'-1']
}
forvalues s=`=`w'+1'/`=`score'+1' {
matrix `G'[`s',`j']=`G'[`s',`j']+exp(-``j'')^`w'*`G'[`=`s'-`w'',`=`j'-1']
}
}
forvalues r=0/`score' {
if "`display'"!="" {
display in green "Level: " in ye `r' in green _col(25) " Gamma=" in ye %10.5f `G'[`=`r'+1',`nbgroups']
}
return scalar gamma`r'= `G'[`=`r'+1',`nbgroups']
}
*return matrix G=`G'
end

40
Modules/ado/personal/g/gammasym.hlp Normal file
View File

@ -0,0 +1,40 @@
{smcl}
{* 9march2007}{...}
{hline}
help for {hi:gammasym}{right:Jean-Benoit Hardouin}
{hline}
{title:The symmetric gamma function}
{p 8 14 2}{cmd:gamma} {it:list} [,{cmdab:w:eight}({it:string}) {cmdab:dis:play}]
{p 4 4 2}{it:list} is a list of the values to use.
{title:Description}
{p 4 4 2}{cmd:gammasym} compute the values of the symmetric gamma function for all the possible level.
{title:Options}
{p 4 4 2}{cmd:weight} defines weights for each value (by default, all the values have a weight of 1)
{p 4 4 2}{cmd:display} displays the value of the symmetric gamma function for each level.
{title:Example}
{cmd:. gammasym -2.1 -1.5 -.7 -.6 1.1 1.00024 2.232447}
{cmd:. gammasym .3 .4 -.1 -.5 0 .6 , weight(1 2 1 1 2 1)}
{cmd:. gammasym 1 1 1 1 1 , w(1 2 3 4 5)}
{title:Outputs}
{p 4 4 2}{cmd:r(gamma#)} return the value of the symmetric gamma function at the level #.
{title:Author}
{p 4 8 2}Jean-Benoit Hardouin, Regional Health Observatory (ORS) - 1, rue Porte Madeleine - BP 2439 - 45032 Orleans Cedex 1 - France.
You can contact the author at {browse "mailto:jean-benoit.hardouin@neuf.fr":jean-benoit.hardouin@neuf.fr} and visit the websites {browse "http://anaqol.free.fr":AnaQol} and {browse "http://freeirt.free.fr":FreeIRT}

181
Modules/ado/personal/g/gausshermite.ado Normal file
View File

@ -0,0 +1,181 @@
*! version 2 15jan2013
************************************************************************************************************
* gausshermite : Estimate an integral of the form |f(x)g(x/mu,sigma)dx or f(x,y)g(x,y/mu,Sigma)dxdy where g(x/mu,sigma) is the distribution function
* of the gaussian distribution of mean mu and variance sigma^2 and g(x,y/mu,Sigma) is the distribution function
* of the bivariate normal distribution of mean mu and covariance matrix Sigma by Gauss Hermite quadratures
*
* Version 1 : May 5, 2005 (Jean-Benoit Hardouin)
* Version 1.1: June 14, 2012 /*name option*/ (Jean-Benoit Hardouin)
* Version 2: January 15, 2013 /*bivariate normal distribution*/ (Jean-Benoit Hardouin, Mohand-Larbi Feddag, Myriam Blanchin)
*
* Jean-Benoit Hardouin, jean-benoit.hardouin@univ-nantes.fr
* EA 4275 "Biostatistics, Pharmacoepidemiology and Subjectives Measures in Health"
* Faculty of Pharmaceutical Sciences - University of Nantes - France
* http://www.sphere-nantes.org
*
* News about this program : http://anaqol.free.fr
* FreeIRT Project : http://freeirt.free.fr
*
* Copyright 2005, 2013 Jean-Benoit Hardouin, Mohand-Larbi Feddag, Myriam Blanchin
*
* 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 gausshermite,rclass
version 7
syntax anything [, Sigma(real -1) Var(string) MU(string) Nodes(integer 12) Display Name(string)]
tempfile gauss
qui capture save `gauss',replace
local save=0
if _rc==0 {
qui save `gauss',replace
local save=1
}
tokenize `anything'
drop _all
tempname mean variance C
qui set obs `=`nodes'*`nodes''
if "`name'"=="" {
if `sigma'!=-1{
if "`var'"==""{
local name x
local nb=1
}
else{
di in red "{p}Please fill in the {hi:name} option{p_end}"
error 198
exit
}
}
else{
if "`var'"!=""{
local name1 x1
local name2 x2
local nb=2
}
else{
di in red "{p}Please fill in the {hi:name} option{p_end}"
error 198
exit
}
}
}
else {
local nb=wordcount("`name'")
if `nb'==2{
local name1=word("`name'",1)
local name2=word("`name'",2)
}
}
if `nb'==2{
capture confirm matrix `var'
if !_rc{
if colsof(`var')==2 & rowsof(`var')==2{
matrix `C'=cholesky(`var')
}
else{
di in red "{p}The covariance matrix in the {hi:var} option should be a 2x2 matrix for a bivariate distribution{p_end}"
error 198
exit
}
}
else{
matrix `variance'=(1,0\0,1)
matrix `C'=cholesky(`variance')
}
}
else{
if `sigma'==-1{
local sig=1
}
else{
local sig=`sigma'
}
}
capture confirm matrix `mu'
if !_rc{
if colsof(`mu')==1 & rowsof(`mu')==1{
local `mean'=`mu'[1,1]
}
else{
matrix `mean'=`mu'
}
}
else{
if "`mu'"==""{
if `nb'==1{
local `mean'=0
}
else{
matrix `mean'=(0,0)
}
}
else{
local `mean'=`mu'
}
}
tempname noeuds poids
qui ghquadm `nodes' `noeuds' `poids'
if `nb'==1{
qui gen `name'=.
qui gen poids=.
forvalues i=1/`nodes' {
qui replace `name'=`noeuds'[1,`i'] in `i'
qui replace poids=`poids'[1,`i'] in `i'
}
qui replace `name'=`name'*(sqrt(2)*`sig')+``mean''
qui gen f=poids/sqrt(_pi)*(`1')
*list `name' poids f in 1/5
}
else{
forvalues i=1/`nb'{
qui gen `name`i''=.
qui gen poids`i'=.
}
local line=1
forvalues i=1/`nodes' {
forvalues j=1/`nodes' {
qui replace `name1'=`noeuds'[1,`i'] *(sqrt(2)*`C'[1,1])+`mean'[1,1] in `line'
qui replace `name2'=`noeuds'[1,`i'] *(sqrt(2)*`C'[2,1])+`noeuds'[1,`j'] *(sqrt(2)*`C'[2,2])+`mean'[1,2] in `line'
qui replace poids1=`poids'[1,`i'] in `line'
qui replace poids2=`poids'[1,`j'] in `line'
local ++line
}
}
qui gen f=poids1*poids2*(`1')/(_pi)
*list `name1' `name2' poids1 poids2 f in 10/20
}
qui su f
return scalar int=r(sum)
if "`display'"!="" {
di in green "int_R (`1')g(`name'/sigma=`sig')d`name'=" in yellow %12.8f `r(sum)'
}
drop _all
if `save'==1 {
qui use `gauss',clear
}
end

43
Modules/ado/personal/g/gausshermite.hlp Normal file
View File

@ -0,0 +1,43 @@
{smcl}
{* 5may2005}{...}
{hline}
help for {hi:gausshermite}{right:Jean-Benoit Hardouin}
{hline}
{title:Estimation of integrals using Gauss Hermite quadratures}
{p 8 14 2}{cmd:gausshermite} {it:function} , {cmdab:s:igma}({it:#}) {cmd:mu}({it:#}) {cmdab:n:odes}({it:#}) {cmdab:d:isplay}
{title:Description}
{p 4 4 2}{cmd:gausshermite} approximates the integrals of the form f(x)g(x/mu,sigma) on all the reals where g(x/mu,sigma) is the gaussian distribution function with mean mu and variance sigma^2.
{title:Options}
{p 4 4 2}{it:function} defines f(x). For example, if f(x)=x^2, {it:function} is x^2. It is necessary to use x for the variable of integration.
{p 4 4 2}{cmd:mu} defines the mean of x (0 by default).
{p 4 4 2}{cmd:sigma} defines the standard deviation of x (1 by default).
{p 4 4 2}{cmd:nodes} defines the number of quadrature nodes (12 by default).
{p 4 4 2}{cmd:display} allows automatically displaying the estimation.
{p 4 4 2}Note that the quadrature nodes and the associated weights are computed using the {cmd:ghquadm} Stata command. Find this command with {stata findit ghquadm:findit ghquadm}.
{title:Example}
{cmd:. gausshermite x^2}
{cmd:. gausshermite x^4+exp(x)-2, sigma(1.5) mu(-.4) d n(10)}
{title:Outputs}
{p 4 4 2}The estimated value of the integral is saved in {cmd:r(int)}.
{title:Author}
{p 4 8 2}Jean-Benoit Hardouin, Regional Health Observatory (ORS) - 1, rue Porte Madeleine - BP 2439 - 45032 Orleans Cedex 1 - France.
You can contact the author at {browse "mailto:jean-benoit.hardouin@orscentre.org":jean-benoit.hardouin@orscentre.org} and visit the websites {browse "http://anaqol.free.fr":AnaQol} and {browse "http://freeirt.free.fr":FreeIRT}

72
Modules/ado/personal/g/gausshermite1.ado Normal file
View File

@ -0,0 +1,72 @@
*! version 1 5may2005
************************************************************************************************************
* gausshermite : Estimate an integral of the form |f(x)g(x/mu,sigma)dx where g(x/mu,sigma) is the distribution function
* of the gaussian distribution of mean mu and variance sigma^2 by Gauss Hermite quadratures
*
* Version 1: May 5, 2005
*
*
* Jean-benoit Hardouin, Regional Health Observatory of Orl<72>ans - France
* jean-benoit.hardouin@orscentre.org
*
* News about this program : http://anaqol.free.fr
* FreeIRT Project : http://freeirt.free.fr
*
* Copyright 2005 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 gausshermite,rclass
version 7
syntax anything [, Sigma(real 1) MU(real 0) Nodes(integer 12) Display]
tempfile gauss
qui capture save `gauss',replace
local save=0
if _rc==0 {
qui save `gauss',replace
local save=1
}
tokenize `anything'
drop _all
qui set obs 100
tempname noeuds poids
qui ghquadm `nodes' `noeuds' `poids'
qui gen x=.
qui gen poids=.
forvalues i=1/`nodes' {
qui replace x=`noeuds'[1,`i'] in `i'
qui replace poids=`poids'[1,`i'] in `i'
}
qui replace x=x*(sqrt(2)*`sigma')+`mu'
qui gen f=poids/sqrt(_pi)*(`1')
qui su f
return scalar int=r(sum)
if "`display'"!="" {
di in green "int_R (`1')g(x/sigma=`sigma')dx=" in yellow %12.8f `r(sum)'
}
drop _all
if `save'==1 {
qui use `gauss',clear
}
end

104
Modules/ado/personal/g/gausshermite2.ado Normal file
View File

@ -0,0 +1,104 @@
*! version 1 03 January 2012
************************************************************************************************************
* gausshermite2 : Estimate an integral of the form : f(x,y)g(x,y/mu,Sigma)dxdy where g(x,y/mu,Sigma) is the distribution function
* of the bivariate normal distribution of mean mu and covariance matrix Sigma by Gauss Hermite quadratures
*
* Version 1: 03 January 2012
*
*
* Mohand Feddag, University of Nantes - France
* Mohand-Larbi.Feddag@univ-nantes.fr *
* News about this program : http://anaqol.free.fr
* FreeIRT Project : http://freeirt.free.fr
*
* Copyright 2012 Mohand Feddag
*
* 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 gausshermite2,rclass
version 11
syntax anything [, Mu(string) Sigma(string) Nodes(integer 12) Display]
tempfile gauss2
qui capture save `gauss2',replace
local save=0
if _rc==0 {
qui save `gauss2',replace
local save=1
}
if "`mu'"=="" {
tempname mu
matrix `mu'=[0,0]
}
if "`sigma'"=="" {
tempname sigma
matrix `sigma'=[1,0\0,1]
}
tokenize `anything'
drop _all
qui set obs `=`nodes'*`nodes''
tempname noeuds poids
qui ghquadm `nodes' `noeuds' `poids'
* Cholesky transformation for the covariance matrix sigma
matrix C=cholesky(`sigma')
*matrix list C
local line=1
qui gen x1=.
qui gen x2=.
qui gen poids1=.
qui gen poids2=.
forvalues i=1/`nodes' {
forvalues j=1/`nodes' {
qui replace x1=`noeuds'[1,`i'] *(sqrt(2)*C[1,1])+`mu'[1,1] in `line'
qui replace x2=`noeuds'[1,`i'] *(sqrt(2)*C[2,1])+`noeuds'[1,`j'] *(sqrt(2)*C[2,2])+`mu'[1,2] in `line'
qui replace poids1=`poids'[1,`i'] in `line'
qui replace poids2=`poids'[1,`j'] in `line'
local ++line
}
}
* Double somme du produit poids[i]*poids[j]*f(x1,x2) qui est affect<63> a la variable sum
qui gen f=poids1*poids2*(`1')/(_pi)
*list x1 x2 poids1 f (sqrt(2)*_pi)
qui su f
local sum=r(sum)
return scalar int=`sum'
if "`display'"!="" {
di in green "int_R^2 (`1')g(x1,x2/mu=`mu',Sigma=`Sigma')dx1dx2=" in yellow %12.8f `sum'
}
drop _all
if `save'==1 {
qui use `gauss2',clear
}
end

134
Modules/ado/personal/g/gausshermite3.ado Normal file
View File

@ -0,0 +1,134 @@
*! version 2 11dec2012
************************************************************************************************************
* gausshermite : Estimate an integral of the form |f(x)g(x/mu,sigma)dx or f(x,y)g(x,y/mu,Sigma)dxdy where g(x/mu,sigma) is the distribution function
* of the gaussian distribution of mean mu and variance sigma^2 and g(x,y/mu,Sigma) is the distribution function
* of the bivariate normal distribution of mean mu and covariance matrix Sigma by Gauss Hermite quadratures
*
* Version 1: May 5, 2005
* Version 1.1: June 14, 2012 /*name option*/
* Version 2: December 11, 2012 /*bivariate normal distribution*/
*
* Jean-benoit Hardouin, Regional Health Observatory of Orl<72>ans - France
* jean-benoit.hardouin@orscentre.org
*
* News about this program : http://anaqol.free.fr
* FreeIRT Project : http://freeirt.free.fr
*
* Copyright 2005, 2012 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 gausshermite3,rclass
version 7
syntax anything [, Sigma(string) MU(string) Nodes(integer 12) Display Name(string)]
tempfile gauss
qui capture save `gauss',replace
local save=0
if _rc==0 {
qui save `gauss',replace
local save=1
}
tokenize `anything'
drop _all
qui set obs `=`nodes'*`nodes''
if "`name'"=="" {
local name x
local nb=1
}
else {
local nb=wordcount("`name'")
if `nb'==2{
di word("`name'",1)
local name1=word("`name'",1)
local name2=word("`name'",2)
}
}
capture confirm matrix `sigma'
if !_rc{
if colsof(`sigma')==1 & rowsof(`sigma')==1{
local sig=`sigma'[1,1]
}
else{
matrix C=cholesky(`sigma')
}
}
else{
local sig=`sigma'
}
capture confirm matrix `mu'
if !_rc{
if colsof(`mu')==1 & rowsof(`mu')==1{
local mean=`mu'[1,1]
}
else{
matrix mean=`mu'
}
}
else{
local mean=`mu'
}
tempname noeuds poids
qui ghquadm `nodes' `noeuds' `poids'
if `nb'==1{
qui gen `name'=.
qui gen poids=.
forvalues i=1/`nodes' {
qui replace `name'=`noeuds'[1,`i'] in `i'
qui replace poids=`poids'[1,`i'] in `i'
}
qui replace `name'=`name'*(sqrt(2)*`sig')+`mean'
qui gen f=poids/sqrt(_pi)*(`1')
*list `name' poids f in 1/5
}
else{
forvalues i=1/`nb'{
qui gen `name`i''=.
qui gen poids`i'=.
}
local line=1
forvalues i=1/`nodes' {
forvalues j=1/`nodes' {
qui replace `name1'=`noeuds'[1,`i'] *(sqrt(2)*C[1,1])+mean[1,1] in `line'
qui replace `name2'=`noeuds'[1,`i'] *(sqrt(2)*C[2,1])+`noeuds'[1,`j'] *(sqrt(2)*C[2,2])+mean[1,2] in `line'
qui replace poids1=`poids'[1,`i'] in `line'
qui replace poids2=`poids'[1,`j'] in `line'
local ++line
}
}
qui gen f=poids1*poids2*(`1')/(_pi)
*list `name1' `name2' poids1 poids2 f in 10/20
}
qui su f
return scalar int=r(sum)
if "`display'"!="" {
di in green "int_R (`1')g(`name'/sigma=`sigma')d`name'=" in yellow %12.8f `r(sum)'
}
drop _all
if `save'==1 {
qui use `gauss',clear
}
end

503
Modules/ado/personal/g/geekel2d.ado Normal file
View File

@ -0,0 +1,503 @@
*! version 4.3 18january2006
************************************************************************************************************
* GEEkel2d: GEE for estimation of unidimensional or 2-dimensional Latent Trait models (Kelderman and Rijkes 1994)
*
* Version 4.3: January 18, 2006 /*Faster version*/
*
* 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
* Version 4 (2004-06-08): Jean-Benoit Hardouin
* Version 4.1 (2005-04-02): Jean-Benoit Hardouin
* Version 4.2 (2005-07-02): Jean-Benoit Hardouin
*
* Use the ghquadm program (findit ghquadm)
*
* Jean-benoit Hardouin, Regional Health Observatory of Orl<72>ans - France
* jean-benoit.hardouin@neuf.fr
*
* News about this program : http://anaqol.free.fr
* FreeIRT Project : http://freeirt.free.fr
*
* Copyright 2003-2006 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=0.25
if `Q'==2 {
local sigmath22=.25
local sigmath12=0.125
}
tempname beta
matrix `beta'=J(`nbitems'+`Q'*(`Q'+1)/2,1,0)
forvalues i=1/`nbitems' {
qui count if ``i''==1
local si`i'=r(N)
local delta``i''=-log(`si`i''/(`N'-`si`i''))
matrix `beta'[`i',1]=`delta``i'''
forvalues j=`i'/`nbitems' {
quiet count if ``j''==1&``i''==1
local si`i'j`j'=r(N)
local si`j'j`i'=r(N)
if "`var'"=="" {
forvalues k=`j'/`nbitems' {
quiet count if ``i''==1&``j''==1&``k''==1
local si`i'j`j'k`k'=r(N)
local si`i'j`k'k`j'=r(N)
local si`j'j`i'k`k'=r(N)
local si`j'j`k'k`i'=r(N)
local si`k'j`i'k`j'=r(N)
local si`k'j`j'k`i'=r(N)
forvalues l=`k'/`nbitems' {
quiet count if ``i''==1&``j''==1&``k''==1&``l''==1
local si`i'j`j'k`k'l`l'=r(N)
local si`i'j`j'k`l'l`k'=r(N)
local si`i'j`k'k`j'l`l'=r(N)
local si`i'j`k'k`l'l`j'=r(N)
local si`i'j`l'k`j'l`k'=r(N)
local si`i'j`l'k`k'l`j'=r(N)
local si`j'j`i'k`k'l`l'=r(N)
local si`j'j`i'k`l'l`k'=r(N)
local si`j'j`k'k`i'l`l'=r(N)
local si`j'j`k'k`l'l`i'=r(N)
local si`j'j`l'k`i'l`k'=r(N)
local si`j'j`l'k`k'l`i'=r(N)
local si`k'j`i'k`j'l`l'=r(N)
local si`k'j`i'k`l'l`j'=r(N)
local si`k'j`j'k`i'l`l'=r(N)
local si`k'j`j'k`l'l`i'=r(N)
local si`k'j`l'k`i'l`j'=r(N)
local si`k'j`l'k`j'l`i'=r(N)
local si`l'j`i'k`j'l`k'=r(N)
local si`l'j`i'k`k'l`j'=r(N)
local si`l'j`j'k`i'l`k'=r(N)
local si`l'j`j'k`k'l`i'=r(N)
local si`l'j`k'k`i'l`j'=r(N)
local si`l'j`k'k`j'l`i'=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 H2i`j'=`B'[`j',1]^2*`sigmath11'+`B'[`j',1]*`B'[`j',2]*`sigmath12'+`B'[`j',2]^2*`sigmath22'
local H4i`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
}
else if `Q'==1 {
local H2i`j'=`B'[`j',1]^2*`sigmath11'
local H4i`j'=3*`B'[`j',1]^4*`sigmath11'^2
}
/* CALCUL DES MOMENTS D'ORDRE 1 ET 2 ET DE LA MATRICE V11*/
local mui`j'=`l1``j'''+`H2i`j''/2*`l3``j'''+`H4i`j''/24*`l5``j'''
local sigmai`j'j`j'=`l2``j'''+`H2i`j''/2*(`l3``j'''-2*`l1``j'''*`l3``j''')+`H4i`j''/24*(`l5``j'''-2*(`l3``j''')^2-2*`l1``j'''*`l5``j''')
}
matrix `V11'=J(`nbitems',`nbitems',0)
forvalues j=1/`nbitems' {
matrix `V11'[`j',`j']=`sigmai`j'j`j''
forvalues l=`=`j'+1'/`nbitems' {
if `Q'==2 {
local H2i`j'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 H4i`j'1j`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 H4i`j'2j`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 H4i`j'3j`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
}
else if `Q'==1 {
local H2i`j'j`l'=`B'[`j',1]*`B'[`l',1]*`sigmath11'
local H4i`j'1j`l'3=3*`B'[`j',1]*`B'[`l',1]^3*`sigmath11'^2
local H4i`j'2j`l'2=3*`B'[`j',1]^2*`B'[`l',1]^2*`sigmath11'^2
local H4i`j'3j`l'1=3*`B'[`l',1]*`B'[`j',1]^3*`sigmath11'^2
}
local H2i`l'j`j'=`H2i`j'j`l''
local H4i`l'1j`j'3=`H4i`j'1j`l'3'
local H4i`l'2j`j'2=`H4i`j'2j`l'2'
local H4i`l'3j`j'1=`H4i`j'3j`l'1'
local sigmai`j'j`l'=`H2i`j'j`l''*(`l2``j'''*`l2``l''')+`H4i`j'1j`l'3'/6*`l2``j'''*`l4``l'''+`H4i`j'3j`l'1'/6*`l4``j'''*`l2``l'''+(`H4i`j'2j`l'2'-`H2i`j''*`H2i`l'')/4*`l3``j'''*`l3``l'''
}
}
/* DEFINITION DE LA MATRICE COMPOCARRE*/
tempname compocarre m
local carre=`nbitems'*(`nbitems'-1)/2
matrix `compocarre'=J(2,`carre',0)
local m=0
forvalues j=1/`nbitems' {
forvalues l=`=`j'+1'/`nbitems' {
local m=`m'+1
matrix `compocarre'[1,`m']=`j'
matrix `compocarre'[2,`m']=`l'
}
}
/* CALCUL DE LA MATRICE V22*/
matrix `V22'=J(`carre',`carre',0)
forvalues k=1/`carre' {
local j=`compocarre'[1,`k']
local l=`compocarre'[2,`k']
matrix `V22'[`k',`k']=(1-2*`mui`j'')*(1-2*`mui`l'')*`sigmai`j'j`l''+`sigmai`j'j`j''*`sigmai`l'j`l''-`sigmai`j'j`l''^2
}
/* CALCUL DES MATRICES V12, V21 ET V*/
matrix `V12'=J(`nbitems',`carre',0)
forvalues k=1/`carre' {
local j=`compocarre'[1,`k']
local l=`compocarre'[2,`k']
matrix `V12'[`j',`k']=(1-2*`mui`j'')*`sigmai`j'j`l''
matrix `V12'[`l',`k']=(1-2*`mui`l'')*`sigmai`j'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)
forvalues j=1/`nbitems' {
matrix `D11'[`j',`j']=-`l2``j'''-`H2i`j''/2*`l4``j'''-`H4i`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'''
}
else 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)
forvalues k=1/`carre' {
local j=`compocarre'[1,`k']
local l=`compocarre'[2,`k']
matrix `D21'[`k',`j']=-`H2i`j'j`l''*`l3``j'''*`l2``l'''-`H4i`j'1j`l'3'/6*`l3``j'''*`l4``l'''-`H4i`j'3j`l'1'/6*`l5``j'''*`l2``l'''-(`H4i`j'2j`l'2'-`H2i`j''*`H2i`l'')/4*`l4``j'''*`l3``l'''
matrix `D21'[`k',`l']=-`H2i`j'j`l''*`l2``j'''*`l3``l'''-`H4i`j'3j`l'1'/6*`l4``j'''*`l3``l'''-`H4i`j'1j`l'3'/6*`l2``j'''*`l5``l'''-(`H4i`j'2j`l'2'-`H2i`j''*`H2i`l'')/4*`l3``j'''*`l4``l'''
tempname tmp1 tmp2 tmp3
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'
}
else 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)
forvalues j=1/`nbitems' {
matrix `chsi'[`j',1]=(`si`j''-`N'*`mui`j'')/`N'
}
forvalues k=1/`carre' {
local j=`compocarre'[1,`k']
local l=`compocarre'[2,`k']
local tmp=`nbitems'+`k'
matrix `chsi'[`tmp',1]=(`si`j'j`l''-`si`j''*`mui`l''-`si`l''*`mui`j''+`N'*`mui`j''*`mui`l''-`N'*`sigmai`j'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 sigm ath11=`beta'[`l',1]
if `Q'==2 {
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=`=`i'+1'/`nbitems' {
matrix `xicarreA'[`i',`j']=`si`i'j`j''-`si`i''*`mui`j''-`si`j''*`mui`i''+`N'*`mui`i''*`mui`j''
matrix `xicarreA'[`i',`j']=`xicarreA'[`j',`i']
}
forvalues col=1/`carre' {
local j=`compocarre'[1,`col']
local k=`compocarre'[2,`col']
matrix `xicarreB'[`i',`col']=`si`i'j`j'k`k''-`mui`i''*`si`j'j`k''-`mui`j''*`si`i'j`k''-`mui`k''*`si`i'j`j''+`mui`i''*`mui`j''*`si`k''+`mui`i''*`mui`k''*`si`j''+`mui`j''*`mui`k''*`si`i''-`N'*`mui`i''*`mui`j''*`mui`k''-`sigmai`j'j`k''*`si`i''+`N'*`mui`i''*`sigmai`j'j`k'''
}
}
forvalues row=1/`carre' {
forvalues col=`row'/`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']=`si`i'j`j'k`k'l`l''-`mui`i''*`si`j'j`k'l`l''-`mui`j''*`si`i'j`k'k`l''-`mui`k''*`si`i'j`j'k`l''-`mui`l''*`si`i'j`j'k`k''+`mui`i''*`mui`j''*`si`k'j`l''+`mui`i''*`mui`k''*`si`j'j`l''+`mui`i''*`mui`l''*`si`j'j`k''+`mui`j''*`mui`k''*`si`i'j`l''+`mui`j''*`mui`l''*`si`i'j`k''+`mui`k''*`mui`l''*`si`i'j`j''-`mui`i''*`mui`j''*`mui`k''*`si`l''-`mui`i''*`mui`j''*`mui`l''*`si`k''-`mui`i''*`mui`k''*`mui`l''*`si`j''-`mui`j''*`mui`k''*`mui`l''*`si`i''-`sigmai`i'j`j''*`si`k'j`l''-`sigmai`k'j`l''*`si`i'j`j''+`sigmai`i'j`j''*`mui`k''*`si`l''+`sigmai`i'j`j''*`mui`l''*`si`k''+`sigmai`k'j`l''*`mui`i''*`si`j''+`sigmai`k'j`l''*`mui`j''*`si`i''+`N'*`mui`i''*`mui`j''*`mui`k''*`mui`l''-`N'*`sigmai`i'j`j''*`mui`k''*`mui`l''-`N'*`sigmai`k'j`l''*`mui`i''*`mui`j''+`N'*`sigmai`i'j`j''*`sigmai`k'j`l''
matrix `xicarreC'[`col',`row']=`xicarreC'[`row',`col']
}
}
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 does not converge"
return scalar error=1
exit
}
if `variat'[1,1]<=`critconv'&`compteur'>0 {
noi di in green _col(8) "The algorithm converges 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 no more converges after the `compteur'th iteration"
}
di ""
if "`var'"=="" {
noi di in green _col(30) "Parameters" in green _col(43) "Standard errors"
forvalues 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])
tempname rho
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"
forvalues 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]
tempname rho
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'"!="" {
tempname noeuds poids
ghquadm `quad' `noeuds' `poids'
tempvar vrais logvrais P
qui gen `P'=0
qui gen `vrais'=0
if `Q'==1 {
forvalues u=1/`quad'{
tempvar vrais`u'
qui gen `vrais`u''=1/sqrt(_pi)
forvalues j=1/`nbitems' {
qui replace `P'=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 `P'=1-`P' if ``j''==0
qui replace `vrais`u''=`vrais`u''*`P'
}
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 sigma
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' {
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 replace `P'=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 `P'=1-`P' if ``j''==0
qui replace `vraisu`u'v`v''=`vraisu`u'v`v''*`P'
}
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 V `W'
}
matrix `beta'=`beta''
return matrix b= `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

78
Modules/ado/personal/g/geekel2d.hlp Normal file
View File

@ -0,0 +1,78 @@
{smcl}
{* 2july2005}{...}
{hline}
help for {hi:geekel2d}{right:Jean-Benoit Hardouin}
{hline}
{title:Estimation of the parameters of undimensional and bidimensional IRT models}
{p 8 14 12}{cmd:geekel2d} {it:varlist} [{cmd:,} {cmdab:coef}({it:matrixname}) {cmdab:nbit}({it:#}) {cmdab:critconv}({it:#}) {cmdab:ll} {cmdab:quad}({it:#}) {cmdab:novar}]
{p 8 14 12}{it:varlist} is a list of two existing dichotomous variables or more.
{title:Description}
{p 4 8 2}{cmd:geekel2d} estimates, by Generalized Estimating
Equations (GEE), the parameters of the model defined by Kelderman (1994) with
one or two dimensions and dichotomic items. This model includes the Rasch model
and the One Parameter Logistic Model (OPLM) for the unidimensional models, the
Multidimensional Generalized Rasch Model (MGRM) and the Multidimensional
Completely Sufficient Rasch Model (MMSRM) for the two-dimensional models.
{title:Options}
{p 4 8 2}{cmd:coef} is the name of a matrix which contains the coeficients B. This
matrix relies the items and the latent traits. Each row represents an item and
there is as many colmuns than the supposed number of latent traits (one or two).
The coefficients are choosen, in general, among the first intergers, but
{cmd:geekel2d} allows using real coefficients. By default, the Rasch model is
supposed (the matrix {cmd: coef} is a vector of 1).
{p 4 8 2}{cmd:nbit} defines the maximal number of iterations in the estimation
algorithm. By default, this number is fixed to 30.
{p 4 8 2}{cmd:critconv} is the value of the convergence criterion, calculated
as the square of the cross-product of the vector containing the difference
between two successive iterations of the parameters estimations. By default,
this criterion is fixed to 1e-15.
{p 4 8 2}{cmd:ll} estimates the marginal log-likelihood and the Akaike
Information Criterion (AIC) by Gauss-Hermite quadratures.
{p 4 8 2}{cmd:quad} defines the number of nodes of quadratures.
{p 4 8 2}{cmd:novar} avoids to compute the standards errors of the estimators (faster).
{title:Remarks}
{p 4 8 2}For detailed informations on the Kelderman model, see Kelderman and
Rijkes (1994) or Adams and al. (1997).
{p 4 8 2}{cmd:geekel2d} don't allows using of polytomous items.
{p 4 8 2}The {cmd:ghquadm} Stata module is needed (use {cmd:findit ghquadm} to obtain it).
{title:Example}
{p 4 8 2}{cmd:. geekel2d item1 item2 item3 item4} /*Rasch model*/
{p 4 8 2}{cmd:. matrix B=(1,0\1,0\0,1\0,1)}
{p 4 8 2}{cmd:. geekel2d item1 item2 item3 item4 , coef(B) nbit(50) critconv(1e-30)}
{title:References}
{p 4 8 2}Kelderman H. and Rijkes C. P. M., Loglinear multidimensional IRT models for polytomously scored items. {it:Psychometrika}, 1994, {it:59}, 149-176.
{p 4 8 2}Adams R. J., Wilson M. R. and Wang W., The multidimensional random coefficient multinomial logit model. {it:Applied Psychological Measurement}, 1997, {it:21}, 1-23.
{title:Author}
{p 4 8 2}Jean-Benoit Hardouin, Regional Health Observatory (ORS) - 1, rue Porte
Madeleine - BP 2439 - 45032 Orleans Cedex 1 - France. You can contact the
author at
{browse "mailto:jean-benoit.hardouin@orscentre.org":jean-benoit.hardouin@orscentre.org}
and visit the websites {browse "http://anaqol.free.fr":AnaQol}
and {browse "http://freeirt.free.fr":FreeIRT}

534
Modules/ado/personal/g/geekel2d4.ado Normal file
View File

@ -0,0 +1,534 @@
************************************************************************************************************
* 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<72>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

158
Modules/ado/personal/g/gengroup.ado Normal file
View File

@ -0,0 +1,158 @@
*! version 1.3 1 June 2013
*! Jean-Benoit Hardouin
*
************************************************************************************************************
* Stata program : genscore
* Generate groups of individals based on the values of an ordinal or continuous variable
*
* Historic
* Version 1 (2007-05-27): Jean-Benoit Hardouin
* Version 1.1 (2007-06-21): Jean-Benoit Hardouin /*Correction of a bug without -if- */
* Version 1.2 (2009-10-28): Jean-Benoit Hardouin /*-continuous- option*/
* Version 1.3 (2013-06-01): Jean-Benoit Hardouin /*Correction in presence of missing values*/
*
* Jean-benoit Hardouin, phD, Assistant Professor
* EA4275 - SPHERE
* Team of Biostatistics, Pharmacoepidemiology and Subjective Measures in Health Sciences
* University of Nantes - Faculty of Pharmaceutical Sciences
* France
* jean-benoit.hardouin@univ-nantes.fr
*
* News about this program :http://www.anaqol.org
*
* Copyright 2007, 2009, 2013 Jean-Benoit Hardouin
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
************************************************************************************************************
program define gengroup ,rclas
version 7.0
syntax varlist(numeric min=1 max=1) [if/] [in] [, NEWvariable(namelist min=1 max=1) REPlace MINsize(integer 30) DETails CONTinuous]
tempvar sort
qui gen `sort'=_n
if "`if'"!="" {
local if2="if `if'"
local if3="&(`if')"
}
marksample touse
if "`newvariable'"=="" {
local newvariable group
}
capture confirm new variable `newvariable'
if _rc!=0&"`replace'"=="" {
di in red "The variable {hi:`newvariable'} is already defined"
exit 198
}
else if _rc!=0&"`replace"!="" {
qui drop `newvariable'
}
if "`continuous'"=="" {
qui gen `newvariable'=`varlist' `if2' `in'
qui su `newvariable' `if2' `in'
local min=r(min)
local max=r(max)
local groupmin=`min'
local groupmax=`min'-1
local numgroup=1
local recode
local list
while (`groupmin'<`max'+1) {
local n=0
while (`n'<`minsize') {
local groupmax=`groupmax'+1
qui count if `newvariable'>=`groupmin'&`newvariable'<=`groupmax'`if3' `in'
local n=r(N)
if `groupmax'>`max' {
local n=`minsize'+1
local numgr
}
}
if `groupmax'<`max' {
local list `list' `groupmax'
}
local recode `recode' `groupmin'/`groupmax'=`numgroup'
if "`details'"!="" {
di in gr "Group " in ye `numgroup' in gr ": Values " in ye `groupmin' in gr " to " in ye `groupmax'
}
local groupmin=`groupmax'+1
local groupmax=`groupmin'-1
local numgroup=`numgroup'+1
}
qui recode `newvariable' `recode' `if2' `in'
qui count if `newvariable'==`numgroup'-1`if3' `in'
local dernier=r(N)
if `dernier'<`minsize' {
qui recode `newvariable' `=`numgroup'-1'=`=`numgroup'-2' `if2' `in'
if "`details'"!="" {
di in gr "The group " in ye `=`numgroup'-1' in gr " is recoded in " in ye `=`numgroup'-2'
}
local list2
forvalues i=1/`=`numgroup'-3' {
local w:word `i' of `list'
local list2 `list2' `w'
}
local list `list2'
}
}
else {
local list
qui sort `varlist'
qui tempvar sort2
qui gen `sort2'=_n
qui gen `newvariable'=0 `if2' `in'
qui replace `newvariable'=. if `varlist'==.
qui count if `varlist'!=.`if3' `in'
local nbind=r(N)
local tmpgp=floor(`nbind'/`minsize')
local minsize=floor(`nbind'/`tmpgp')
local nbused=0
tempvar used
qui gen `used'=0 `if2' `in'
local num=1
while (`=`nbused'+`minsize''<`nbind'+1) {
qui su `varlist' if `sort2'==`=`nbused'+`minsize''`if3' `in'
local mean=r(mean)
local list `list' `mean'
local mean=round(`mean',0.0000001)+0.0000001
if "`details'"!="" {
di in gr "The values inferior to " in ye `mean' in gr " are recoded in " in ye `num'
}
qui replace `newvariable'=`num' if `varlist'<=`mean'&`used'==0
qui replace `used'=1 if `newvariable'!=0&`newvariable'!=.
qui count if `used'==1
local nbused=`r(N)'
local num=`num'+1
}
qui replace `newvariable'=`num'-1 if `newvariable'==0
if "`details'"!="" {
qui su `varlist'
local max=r(max)
if `max'>`mean' {
di in gr "Individuals between " in ye `mean' in gr "and " in ye `max' " are recoded in " in ye `=`num'-1'
}
}
qui sort `sort'
}
return local list `list'
end

51
Modules/ado/personal/g/gengroup.hlp Normal file
View File

@ -0,0 +1,51 @@
{smcl}
{* 6May2013}{...}
{hline}
help for {hi:gengroup}{right:Jean-benoit Hardouin}
{hline}
{title:Module to generate group of individuals based on a ordinal or continuous variable}
{p 8 14 2}{cmd:gengroup} {it:varname} [{cmd:if} {it:exp}] [{cmd:in} {it:range}] [{cmd:,} {cmdab:new:variable}({it:newvarname}) {cmdab:rep:lace} {cmdab:min:size}(#) {cmdab:det:ails} {cmdab:cont:inuous}]
{title:Description}
{p 4 8 2}{cmd:gengroup} creates groups of individuals by using the values of an ordinal or continuous variable.
The module creates groups by recoding several adjacent values of the variable, until obtaining groups with more than individuals than the number defined in the {cmd:minsize} option.
{title:Options}
{p 4 8 2}{cmd:newvariable} defines the name of the new variable ({it:group} by default).
{p 4 8 2}{cmd:replace} replaces the variable defined in the {cmd:newvariable} option if it already exists.
{p 4 8 2}{cmd:minsize} defines the minimal number of individuals in each group (30 by default).
{p 4 8 2}{cmd:details} diplays the composition of each group.
{p 4 8 2}{cmd:continuous} allows handling a continuous variable.
{title:Examples}
{p 4 8 2}{inp:. gengroup score}
{p 4 8 2}{inp:. gengroup score, newvariable(grouptocreate) replace minsize(80)}
{p 4 8 2}{inp:. gengroup score, details}
{title:Author}
{p 4 8 2}Jean-Benoit Hardouin, PhD, assistant professor{p_end}
{p 4 8 2}Team of Biostatistics, Clinical Research and Subjective Measures in Health Sciences{p_end}
{p 4 8 2}University of Nantes - Faculty of Pharmaceutical Sciences{p_end}
{p 4 8 2}1, rue Gaston Veil - BP 53508{p_end}
{p 4 8 2}44035 Nantes Cedex 1 - FRANCE{p_end}
{p 4 8 2}Email:
{browse "mailto:jean-benoit.hardouin@univ-nantes.fr":jean-benoit.hardouin@univ-nantes.fr}{p_end}
{p 4 8 2}Website {browse "http://www.anaqol.org":AnaQol}
{title:Also see}
{p 4 13 2}Online: help for {help egen}, {help generate} and {help genscore} if installed.{p_end}

133
Modules/ado/personal/g/gengutt.ado Normal file
View File

@ -0,0 +1,133 @@
*! Version 1.1 24 March 2015
*! Jean-Benoit Hardouin
************************************************************************************************************
* Stata program : gengutt
* Generation of Guttman's error per individuals
* Release 1.1 : March 24, 2015 /*no required -order- option*/
*
*
* Jean-benoit Hardouin, phD, Assistant Professor
* Team of Biostatistics, Pharmacoepidemiology and Subjective Measures in Health Sciences (UPRES EA 4275 SPHERE)
* University of Nantes - Faculty of Pharmaceutical Sciences
* France
* jean-benoit.hardouin@anaqol.org
*
* Release 1 : March 20, 2015 /*initial version with only -order- option*/
*
* News about this program :http://www.anaqol.org
*
* Copyright 2015 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 gengutt , rclass
syntax varlist(min=2 numeric) [if] , Gen(string) [ Replace Order(string) PAIRWise] /*[,PAIRWise PAIR ppp pmm noADJust GENERror(string) REPlace GRaph MONOtonicity(string) NIRESTscore(string) NIPmatrix(string) HTML(string)]*/
preserve
tokenize `varlist'
local nbitems : word count `varlist'
local nbs=0
forvalues i=1/`nbitems' {
qui su ``i''
local modamax`i'=r(max)
local nbs=`nbs'+`modamax`i''
local testmodamax`i'=0
local ordermax`i'=0
}
if "`order'"=="" {
tempname mattmp
matrix `mattmp'=J(3,`nbs',0)
local col=1
forvalues i=1/`nbitems' {
forvalues s=1/`modamax`i'' {
qui count if ``i''>=`s'
matrix `mattmp'[3,`col']=r(N)
matrix `mattmp'[1,`col']=`i'
matrix `mattmp'[2,`col']=`s'
local ++col
}
}
local list
forvalues j=1/`nbs' {
local max=0
local col=0
local item=0
forvalues i=1/`nbs' {
if `mattmp'[3,`i']>`max' {
local max=`mattmp'[3,`i']
local col=`i'
local item=`mattmp'[1,`i']
}
}
local list `list' `item'
matrix `mattmp'[3,`col']=0
}
local order `list'
}
local nbsteps: word count `order'
tempname matorder
matrix `matorder'=J(3,`nbsteps',0)
local defordre=1
forvalues i=1/`nbsteps' {
local tmp: word `i' of `order'
local ++testmodamax`tmp'
matrix `matorder'[1,`i']=`tmp'
local ordermax`tmp'=`ordermax`tmp''+1
matrix `matorder'[2,`i']=`ordermax`tmp''
local i`tmp'_m`ordermax`tmp''=`defordre'
local s`defordre' i`tmp'_m`ordermax`tmp''
tempname `s`defordre''
qui gen ``s`defordre'''=``tmp''>=`ordermax`tmp''
local ++defordre
}
local error=0
forvalues i=1/`nbitems' {
if `modamax`i''>`testmodamax`i'' {
di in red "The item `i' presents more answer categories than defined in the -order- option "
local error=1
}
}
if `error'!=0 {
error 198
}
if "`replace'"!="" {
qui capture drop `gen'
}
qui gen `gen'=0
forvalues i=1/`nbitems' {
if "`replace'"!="" {
qui capture drop `gen'_``i''
}
qui gen `gen'_``i''=0
}
forvalues i=1/`nbsteps' {
forvalues j=`=`i'+1'/`nbsteps' {
qui replace `gen'=`gen'+1 if ``s`i'''<``s`j'''&``s`i'''!=.&``s`j'''!=.
local iti=`matorder'[1,`i']
local itj=`matorder'[1,`j']
qui replace `gen'_``iti''=`gen'_``iti''+1 if ``s`i'''<``s`j'''&``s`i'''!=.&``s`j'''!=.
qui replace `gen'_``itj''=`gen'_``itj''+1 if ``s`i'''<``s`j'''&``s`i'''!=.&``s`j'''!=.
}
}
if "`replace'"!="" {
capture drop `gen'_gp
}
egen `gen'_gp=concat(``s1''-``s`nbsteps''')
capture restore, not
end

86
Modules/ado/personal/g/genscore.ado Normal file
View File

@ -0,0 +1,86 @@
*! version 1.4 27december2005
*! Jean-Benoit Hardouin
*
************************************************************************************************************
* Stata program : genscore
* Generate scores from a list of variables
* Version 1.4 : December 27, 2005 /*corrects a bug with the mean option*/
*
* Historic
* Version 1.2 (2005-10-01): Jean-Benoit Hardouin
* Version 1.3 (2005-12-09): Jean-Benoit Hardouin /*centered and standardized options*/
*
* Jean-benoit Hardouin, Regional Health Observatory of Orl<72>ans - France
* jean-benoit.hardouin@orscentre.org
*
* News about this program : http://anaqol.free.fr
* FreeIRT Project : http://freeirt.free.fr
*
* Copyright 2005 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 genscore
version 7.0
syntax varlist(min=1) [if] [in] [fweight] [, SCore(namelist min=1 max=1) CENTered STAndardized MEan MIssing(string) REPlace]
marksample touse
if "`score'"=="" {
local score score
}
local nbitems:word count `varlist'
tokenize `varlist'
if "`missing'"=="" {
local missing .
}
capture confirm new variable `score'
quietly {
if _rc!=0&"`replace"=="" {
di in red "The variable {hi:`score'} already defined"
exit 198
}
else if _rc!=0&"`replace"!="" {
drop `score'
}
forvalues i=1/`nbitems' {
tempname var`i'
local sd=1
local moy=0
if "`standardized'"!=""|"`centered'"!="" {
su ``i'' [`weight'`exp']
local moy=r(mean)
local sd=r(sd)
if "`standardized'"=="" {
local sd=1
}
}
gen `var`i''=(``i''-`moy')/`sd'
}
gen `score'=0 if `touse'
forvalues i=1/`nbitems' {
replace `score'=`score'+`var`i'' if `touse'
replace `score'=. if `touse'&``i''==`missing'&``i''>=.
}
if "`mean'"!="" {
replace `score'=`score'/`nbitems' if `touse'
}
}
end

54
Modules/ado/personal/g/genscore.hlp Normal file
View File

@ -0,0 +1,54 @@
{smcl}
{* 9december2005}{...}
{hline}
help for {hi:genscore}{right:Jean-benoit Hardouin}
{hline}
{title:Module to generate a score}
{p 8 14 2}{cmd:genscore} {it:varlist} [{cmd:if} {it:exp}] [{cmd:in} {it:range}] [{it:weights}] [{cmd:,} {cmdab:sc:ore}({it:newvarname}) {cmdab:me:an} {cmdab:cent:ered} {cmdab:sta:ndardized} {cmdab:mi:ssing}(string) {cmdab:rep:lace}]
{title:Description}
{p 4 8 2}{cmd:genscore} allows creating a new variable containing the score
computed as the sum (or the mean) of the responses to the variables defined in
{it:varlist}. Variables of {it:varlist} can be centered or standardized. The new
variable will contain a missing value for each individual with at least one
missing values among the variables of {it:varlist}.
{p 4 8 2}Only {cmd:fweights} are allowed, and the weights have an influence only with
the {cmd:centered} or {cmd:standardized} options.
{title:Options}
{p 4 8 2}{cmd:score} defines the name of the new variable ({it:score} by default).
{p 4 8 2}{cmd:mean} allows obtaining a mean score (the sum of the responses is divided by the number of variables).
{p 4 8 2}{cmd:centered} computes the scores with centered variables.
{p 4 8 2}{cmd:standardized} computes the scores with standardized variables.
{p 4 8 2}{cmd:missing} defines a specific modality of the variables which can be considered as a missing value.
{p 4 8 2}{cmd:replace} allows replacing the variable defined in the {cmd:score} option if it already exists.
{title:Examples}
{p 4 8 2}{cmd:. genscore turn trunk mpg [fweight=w],standardized mean}
{p 4 8 2}{cmd:. genscore item1-item20, score(totscore) missing(9)}
{p 4 8 2}{cmd:. genscore var*, score(scorevar) replace}
{title:Author}
{p 4 8 2}Jean-Benoit Hardouin, Regional Health Observatory (ORS) - 1, rue Porte
Madeleine - BP 2439 - 45032 Orleans Cedex 1 - France. You can contact the author at
{browse "mailto:jean-benoit.hardouin@orscentre.org":jean-benoit.hardouin@orscentre.org}
and visit the websites {browse "http://anaqol.free.fr":AnaQol} and
{browse "http://freeirt.free.fr":FreeIRT}
{title:Also see}
{p 4 8 2}Online: help for {help egen}, {help generate} and {help mvdecode}{p_end}

181
Modules/ado/personal/g/guttmax.ado Normal file
View File

@ -0,0 +1,181 @@
*! Version 1 19 November 2008
*! Jean-Benoit Hardouin
************************************************************************************************************
* Stata program : guttmax
* Research of the maximal number of Guttman Errors for a specific score
*
* Historic :
* Version 1 (November 19, 2008) [Jean-Benoit Hardouin]
*
* Jean-benoit Hardouin, phD, Assistant Professor
* Team of Biostatistics, Clinical Research and Subjective Measures in Health Sciences (UPRES EA 4275)
* University of Nantes - Faculty of Pharmaceutical Sciences
* France
* jean-benoit.hardouin@anaqol.org
*
* Requiered Stata modules:
* -anaoption- (version 1)
* -traces- (version 3.2)
* -gengroup- (version 1)
*
* News about this program :http://www.anaqol.org
* FreeIRT Project website : http://www.freeirt.org
*
* Copyright 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
*
***********************************************************************************************************
program define guttmax , rclass
version 10
syntax anything [,Score(real -1) Response(string)]
if `score'!=-1&"`response'"!="" {
di in red "You cannot defined in the same time the score and the response options";
error 198
}
local step:word count `anything'
local nbitems=0
forvalues i=1/`step' {
local step`i':word `i' of `anything'
if `step`i''>`nbitems' {
local nbitems=`step`i''
}
}
di in green "Number of steps: " in ye `step'
di in green "Number of items: " in ye `nbitems'
tempname mate
matrix `mate'=J(1,`step',0)
if `score'!=-1 {
if `score'>`step' {
di in red "You cannot defined a score greater to the maximal possible score";
error 198
}
di in green "Score: " in ye `score'
local maxstep=0
forvalues i=1/`nbitems' {
local nstep`i'=0
forvalues j=1/`step' {
if `step`j''==`i' {
local ++nstep`i'
if `nstep`i''>`maxstep' {
local maxstep=`nstep`i''
}
}
}
}
tempname calcul
matrix `calcul'=J(`nbitems',`maxstep',0)
*matrix list `calcul'
forvalues i=1/`nbitems' {
*di "item `i'"
local n=1
forvalues j=1/`step' {
if `step`j''==`i' {
forvalues s=`n'/`nstep`i'' {
matrix `calcul'[`i',`s']=`calcul'[`i',`s']+`j'
matrix list `calcul'
}
local ++n
}
}
forvalues j=2/`nstep`i'' {
matrix `calcul'[`i',`j']=`calcul'[`i',`j']/`j'
di "j=`j'"
matrix list `calcul'
}
}
*matrix list `calcul'
while (`score'>0) {
local max=0
forvalues i=1/`nbitems' {
local s=min(`score',`nstep`i'')
*di "forvalues j=1/`s' {"
forvalues j=1/`s' {
if `calcul'[`i',`j']>`max' {
local maxi=`i'
local maxj=`j'
local max=`calcul'[`i',`j']
}
}
}
*di "maxi=`maxi' maxj=`maxj'"
local d=0
forvalues l=1/`step' {
if `step`l''==`maxi'&`mate'[1,`l']==0&`d'<=`maxj' {
matrix `mate'[1,`l']=1
local ++d
}
}
forvalues l=1/`maxj' {
matrix `calcul'[`maxi',`l']=0
}
local score=`score'-`maxj'
*matrix list `calcul'
}
}
else { /*if reponse*/
local nbrep: word count `response'
if `nbrep'!=`nbitems' {
di in red "You cannot defined a number of responses different of the number of items";
error 198
}
tempname cpt
qui matrix `cpt'=J(`nbitems',1,0)
forvalues i=1/`nbitems' {
local rep`i':word `i' of `response'
qui matrix `cpt'[`i',1]=`rep`i''
}
forvalues i=1/`step' {
if `cpt'[`step`i'',1]>0 {
qui matrix `mate'[1,`i']=1
qui matrix `cpt'[`step`i'',1]=`cpt'[`step`i'',1]-1
}
else {
qui matrix `mate'[1,`i']=0
}
}
}
local emax=0
forvalue i=1/`step' {
forvalues j=`=`i'+1'/`step' {
if `mate'[1,`i']==0&`mate'[1,`j']==1 {
local ++emax
}
}
}
if `score'!=-1 {
di in green "Responses profile generating the most important number of Guttman error"
matrix list `mate' ,noheader nonames
di in green "Max number of Guttman errors : " in ye `emax'
}
else {
di in green "Responses profile"
matrix list `mate' ,noheader nonames
di in green "Number of Guttman errors : " in ye `emax'
}
return scalar maxegutt =`emax'
end

136
Modules/ado/personal/g/guttmax1.ado Normal file
View File

@ -0,0 +1,136 @@
*! Version 1 19 November 2008
*! Jean-Benoit Hardouin
************************************************************************************************************
* Stata program : guttmax
* Research of the maximal number of Guttman Errors for a specific score
*
* Historic :
* Version 1 (November 19, 2008) [Jean-Benoit Hardouin]
*
* Jean-benoit Hardouin, phD, Assistant Professor
* Team of Biostatistics, Clinical Research and Subjective Measures in Health Sciences (UPRES EA 4275)
* University of Nantes - Faculty of Pharmaceutical Sciences
* France
* jean-benoit.hardouin@anaqol.org
*
* Requiered Stata modules:
* -anaoption- (version 1)
* -traces- (version 3.2)
* -gengroup- (version 1)
*
* News about this program :http://www.anaqol.org
* FreeIRT Project website : http://www.freeirt.org
*
* Copyright 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
*
***********************************************************************************************************
program define guttmax1 , rclass
version 10
syntax anything [,Score(int 0)]
local step:word count `anything'
local nbitems=0
forvalues i=1/`step' {
local step`i':word `i' of `anything'
if `step`i''>`nbitems' {
local nbitems=`step`i''
}
}
di in green "Number of steps: " in ye `step'
di in green "Number of items: " in ye `nbitems'
di in green "Score: " in ye `score'
tempname mate
matrix `mate'=J(1,`step',0)
local maxstep=0
forvalues i=1/`nbitems' {
local nstep`i'=0
forvalues j=1/`step' {
if `step`j''==`i' {
local ++nstep`i'
if `nstep`i''>`maxstep' {
local maxstep=`nstep`i''
}
}
}
}
tempname calcul
matrix `calcul'=J(`nbitems',`maxstep',0)
*matrix list `calcul'
forvalues i=1/`nbitems' {
*di "item `i'"
local n=1
forvalues j=1/`step' {
if `step`j''==`i' {
forvalues s=`n'/`nstep`i'' {
matrix `calcul'[`i',`s']=`calcul'[`i',`s']+`j'
}
local ++n
}
}
forvalues j=2/`nstep`i'' {
matrix `calcul'[`i',`j']=`calcul'[`i',`j']/`j'
}
}
*matrix list `calcul'
while (`score'>0) {
local max=0
forvalues i=1/`nbitems' {
local s=min(`score',`nstep`i'')
*di "forvalues j=1/`s' {"
forvalues j=1/`s' {
if `calcul'[`i',`j']>`max' {
local maxi=`i'
local maxj=`j'
local max=`calcul'[`i',`j']
}
}
}
*di "maxi=`maxi' maxj=`maxj'"
local d=0
forvalues l=1/`step' {
if `step`l''==`maxi'&`mate'[1,`l']==0&`d'<=`maxj' {
matrix `mate'[1,`l']=1
local ++d
}
}
forvalues l=1/`maxj' {
matrix `calcul'[`maxi',`l']=0
}
local score=`score'-`maxj'
*matrix list `calcul'
}
local emax=0
forvalue i=1/`step' {
forvalues j=`=`i'+1'/`step' {
if `mate'[1,`i']==0&`mate'[1,`j']==1 {
local ++emax
}
}
}
di in green "Responses profile generating the most important number of Guttman error"
matrix list `mate' ,noheader nonames
di in green "Max number of Guttman errors : " in ye `emax'
return scalar maxegutt =`emax'
end