Added a few .ado (raschpower)
parent
b16a752e24
commit
ef1fbfbd4f
@ -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
|
||||
|
||||
|
||||
@ -0,0 +1,85 @@
|
||||
program define ghquadm
|
||||
* stolen from gllamm6 who stole it from rfprobit (Bill Sribney)
|
||||
version 4.0
|
||||
|
||||
parse "`*'", parse(" ")
|
||||
local n = `1'
|
||||
if `n' + 2 > _N {
|
||||
di in red /*
|
||||
*/ "`n' + 2 observations needed to compute quadrature points"
|
||||
exit 2001
|
||||
}
|
||||
|
||||
tempname x w xx ww a b
|
||||
local i 1
|
||||
local m = int((`n' + 1)/2)
|
||||
matrix x = J(1,`m',0)
|
||||
matrix w = x
|
||||
while `i' <= `m' {
|
||||
if `i' == 1 {
|
||||
scalar `xx' = sqrt(2*`n'+1)-1.85575*(2*`n'+1)^(-1/6)
|
||||
}
|
||||
else if `i' == 2 { scalar `xx' = `xx'-1.14*`n'^0.426/`xx' }
|
||||
else if `i' == 3 { scalar `xx' = 1.86*`xx'-0.86*x[1,1] }
|
||||
else if `i' == 4 { scalar `xx' = 1.91*`xx'-0.91*x[1,2] }
|
||||
else {
|
||||
local im2 = `i' -2
|
||||
scalar `xx' = 2*`xx'-x[1,`im2']
|
||||
}
|
||||
hermite `n' `xx' `ww'
|
||||
matrix x[1,`i'] = `xx'
|
||||
matrix w[1,`i'] = `ww'
|
||||
local i = `i' + 1
|
||||
}
|
||||
if mod(`n', 2) == 1 { matrix x[1,`m'] = 0}
|
||||
/* start in tails */
|
||||
matrix `b' = (1,1)
|
||||
matrix w = w#`b'
|
||||
matrix w = w[1,1..`n']
|
||||
matrix `b' = (1,-1)
|
||||
matrix x = x#`b'
|
||||
matrix x = x[1,1..`n']
|
||||
/* other alternative (left to right) */
|
||||
/*
|
||||
above: matrix x = J(1,`n',0)
|
||||
while ( `i'<=`n'){
|
||||
matrix x[1, `i'] = -x[1, `n'+1-`i']
|
||||
matrix w[1, `i'] = w[1, `n'+1-`i']
|
||||
local i = `i' + 1
|
||||
}
|
||||
*/
|
||||
matrix `2' = x
|
||||
matrix `3' = w
|
||||
end
|
||||
|
||||
|
||||
program define hermite /* integer n, scalar x, scalar w */
|
||||
* stolen from gllamm6 who stole it from rfprobit (Bill Sribney)
|
||||
version 4.0
|
||||
local n "`1'"
|
||||
local x "`2'"
|
||||
local w "`3'"
|
||||
local last = `n' + 2
|
||||
tempname i p
|
||||
matrix `p' = J(1,`last',0)
|
||||
scalar `i' = 1
|
||||
while `i' <= 10 {
|
||||
matrix `p'[1,1]=0
|
||||
matrix `p'[1,2] = _pi^(-0.25)
|
||||
local k = 3
|
||||
while `k'<=`last'{
|
||||
matrix `p'[1,`k'] = `x'*sqrt(2/(`k'-2))*`p'[1,`k'-1] /*
|
||||
*/ - sqrt((`k'-3)/(`k'-2))*`p'[1,`k'-2]
|
||||
local k = `k' + 1
|
||||
}
|
||||
scalar `w' = sqrt(2*`n')*`p'[1,`last'-1]
|
||||
scalar `x' = `x' - `p'[1,`last']/`w'
|
||||
if abs(`p'[1,`last']/`w') < 3e-14 {
|
||||
scalar `w' = 2/(`w'*`w')
|
||||
exit
|
||||
}
|
||||
scalar `i' = `i' + 1
|
||||
}
|
||||
di in red "hermite did not converge"
|
||||
exit 499
|
||||
end
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@ -0,0 +1,908 @@
|
||||
*! version 5.12 : Novembre 19, 2021
|
||||
*! Jean-Benoit Hardouin, Myriam Blanchin, Bastien Perrot
|
||||
************************************************************************************************************
|
||||
* raschpower: Estimation of the power of the Wald test in order to compare the means of the latent trait in two groups of individuals
|
||||
*
|
||||
* Version 1 : January 25th, 2010 (Jean-Benoit Hardouin)
|
||||
* Version 1.1 : January 26th, 2010 (Jean-Benoit Hardouin)
|
||||
* Version 1.2 : November 1st, 2010 (Jean-Benoit Hardouin)
|
||||
* version 1.3 : May 2nd, 2011 (Jean-Benoit Hardouin)
|
||||
* version 1.4 : July 7th, 2011 (Jean-Benoit Hardouin) : minor corrections
|
||||
* version 1.5 : July 11th, 2011 (Jean-Benoit Hardouin) : minor corrections
|
||||
* version 2 : August 30th, 2011 (Jean-Benoit Hardouin, Myriam Blanchin) : corrections
|
||||
* version 3 : October 18th, 2011 (Jean-Benoit Hardouin, Myriam Blanchin) : Extension to the PCM, -method- option, -nbpatterns- options, changes in the presentation of the results
|
||||
* version 3.1 : October 25th, 2011 (Jean-Benoit Hardouin, Myriam Blanchin) : POPULATION+GH method
|
||||
* version 3.2 : February 6th, 2012 (Jean-Benoit Hardouin, Myriam Blanchin) : minor corrections
|
||||
* version 4 : January 30th, 2013 (Jean-Benoit Hardouin, Myriam Blanchin) : Extension to longitudinal design
|
||||
* version 4.1 : June 3rd, 2013 (Jean-Benoit Hardouin, Myriam Blanchin) : detail option
|
||||
* version 5 : October 22th, 2013 (Jean-Benoit Hardouin, Myriam Blanchin) : Extension to longitudinal design with polytomous items
|
||||
* version 5.1 : January 21th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin) : unequal variance between groups
|
||||
* version 5.2 : January 22th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin) : unequal variance between groups
|
||||
* version 5.3 : January 29th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin) : correction of a bug
|
||||
* version 5.4 : April 4th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin, Bastien Perrot) : HTML option, graph option, minor corrections for results display
|
||||
* version 5.5 : April 9th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin, Bastien Perrot) : graph option for longitudinal design
|
||||
* version 5.6 : April 10th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin, Bastien Perrot) : corrections of bugs
|
||||
* version 5.7 : June 21th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin, Bastien Perrot) : minor corrections
|
||||
* version 5.8 : July 9th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin, Bastien Perrot) : correction of a bug when the variances are different in the 2 groups (transversal)
|
||||
* version 5.9 : September 26th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin, Bastien Perrot) : fixed of a bug with the graph option (correlation could be greater than one)
|
||||
* version 5.10 : July 18th, 2019 (Jean-Benoit Hardouin) : filesave and dirsave options
|
||||
* version 5.11 : July 24th, 2019 (Jean-Benoit Hardouin) : pcm->pcmold
|
||||
* version 5.12 : July 24th, 2019 (Myriam Blanchin) : correction of bug (essai.dta)
|
||||
|
||||
* Jean-benoit Hardouin, jean-benoit.hardouin@univ-nantes.fr
|
||||
* Myriam Blanchin, myriam.blanchin@univ-nantes.fr
|
||||
* Bastien Perrot, Bastien.perrot@univ-nantes.fr
|
||||
* INSERM UMR 1246-SPHERE "Methods in Patient Centered Outcomes and Health Research", Nantes University, University of Tours
|
||||
* http://www.sphere-nantes.org
|
||||
*
|
||||
* News about this program : http://www.anaqol.org - anaqol@sphere-nantes.fr
|
||||
*
|
||||
* Copyright 2010-2014, 2019 Jean-Benoit Hardouin, Myriam Blanchin, Bastien Perrot
|
||||
*
|
||||
* 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 raschpower,rclass
|
||||
syntax [varlist] [, n0(int 100) n1(int 100) Gamma(real .5) Difficulties(string) Var(string) Method(string) NBPatterns(int 2) nodata EXPectedpower(real -1) LONGitudinal freevar freeitems DETail HTML(string) graph gvar(real -1) gcorr(real -1) min(real 0.25) step(real 1) max(real 9) FILESave DIRSave(string) pro]
|
||||
version 11
|
||||
|
||||
if "`pro'"!="" {
|
||||
di "START"
|
||||
}
|
||||
|
||||
if "`html'" != "" {
|
||||
set scheme sj
|
||||
local htmlregion "graphregion(fcolor(white) ifcolor(white))"
|
||||
di "<!-- SphereCalc start of response -->"
|
||||
di "<pre>"
|
||||
}
|
||||
|
||||
tempfile raschpowerfile
|
||||
capture qui save "`raschpowerfile'",replace
|
||||
tempname db d dlong matvar
|
||||
|
||||
/*******************************************************************************
|
||||
DEFINITION DES PARAMETRES
|
||||
*******************************************************************************/
|
||||
|
||||
if "`graph'" != "" & "`longitudinal'" != "" {
|
||||
if `gvar' != -1 & `gcorr' != -1 {
|
||||
di in red "You cannot use both {hi:gvar} and {hi:gcorr}"
|
||||
error 198
|
||||
}
|
||||
if `gvar' == -1 & `gcorr' == -1 {
|
||||
di in red "You must precise values for {hi:gvar} or {hi:gcorr} if you use the {hi:graph} option"
|
||||
error 198
|
||||
}
|
||||
}
|
||||
|
||||
if "`difficulties'"=="" {
|
||||
matrix `d'=[-1\-0.5\0\0.5\1]
|
||||
}
|
||||
else {
|
||||
matrix `d'=`difficulties'
|
||||
}
|
||||
local nbitems=rowsof(`d')
|
||||
local nbmodat=colsof(`d')+1
|
||||
|
||||
if "`longitudinal'"==""{ /*CAS TRANSVERSAL*/
|
||||
local nbt=1
|
||||
local nbtotitems=`nbitems'
|
||||
local nbpatmax=2*`nbmodat'^`nbitems'
|
||||
if "`var'"==""{
|
||||
local var=1
|
||||
local var0=1
|
||||
local var1=1
|
||||
}
|
||||
else{
|
||||
capture confirm number `var'
|
||||
if !_rc {
|
||||
local var0=`var'
|
||||
local var1=`var'
|
||||
}
|
||||
else{
|
||||
local nbw:word count `var'
|
||||
local grp=1
|
||||
if `nbw'==2 {
|
||||
local t1: word 1 of `var'
|
||||
capture confirm number `t1'
|
||||
if !_rc {
|
||||
local t2: word 2 of `var'
|
||||
capture confirm number `t2'
|
||||
if !_rc {
|
||||
local var0=`t1'
|
||||
local var1=`t2'
|
||||
local grp=2
|
||||
}
|
||||
}
|
||||
}
|
||||
if `nbw'!=2|`grp'==1 {
|
||||
capture confirm matrix `var'
|
||||
if !_rc & colsof(`matvar')==1 & rowsof(`matvar')==1{
|
||||
matrix `matvar'=`var'
|
||||
local var0=`matvar'[1,1]
|
||||
local var1=`matvar'[1,1]
|
||||
}
|
||||
else{
|
||||
di in red "{p}The {hi:var} option should contain a number or a 1x1 matrix for transversal studies{p_end}"
|
||||
error 198
|
||||
exit
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
matrix `dlong'=`d' /*DEMANDER A MYRIAM A QUOI CA SERT CAR ON EST EN TRANSVERSAL*/
|
||||
local varc=((`n0'-1)*`var0'+(`n1'-1)*`var1')/(`n0'+`n1'-2)
|
||||
|
||||
local sumd=0
|
||||
forvalues numit=1/`nbitems'{
|
||||
local sumd=`sumd'+`d'[`numit',1]
|
||||
}
|
||||
local sumd=`sumd'/`nbitems'
|
||||
if `=abs(`sumd')'>=`=2*sqrt(`varc')'{
|
||||
if "`html'" == "" {
|
||||
di "{p}WARNING: It is not recommended to use Raschpower when the gap between the global mean of the latent trait (fixed to 0) and the mean of the item parameters is greater than or equal to 2 standard deviation of the latent trait. {p_end}"
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
else{ /*CAS LONGITUDINAL*/
|
||||
local nbt=2
|
||||
local nbtotitems=2*`nbitems'
|
||||
local nbpatmax=`nbmodat'^(`nbitems'*2)
|
||||
local n1=0
|
||||
local mean1=0
|
||||
local mean2=`mean1'+`gamma'
|
||||
if "`var'"==""{
|
||||
matrix `matvar'=(1,0\0,1)
|
||||
}
|
||||
else{
|
||||
matrix `matvar'=`var'
|
||||
capture confirm matrix `matvar'
|
||||
if _rc | colsof(`matvar')!=2 | rowsof(`matvar')!=2{
|
||||
di in red "{p}The {hi:var} option should contain a 2x2 matrix for longitudinal studies{p_end}"
|
||||
error 198
|
||||
exit
|
||||
}
|
||||
}
|
||||
matrix `dlong'=J(`nbtotitems',`=`nbmodat'-1',.)
|
||||
matrix `dlong'[1,1]=`d'
|
||||
matrix `dlong'[`=`nbitems'+1',1]=`d'
|
||||
}
|
||||
|
||||
/*
|
||||
if "`max'" == "" {
|
||||
local max = 9 + `step'/2
|
||||
}
|
||||
else {
|
||||
local max = `max' + `step'/2
|
||||
}*/
|
||||
|
||||
forval v = `min'(`step')`max' {
|
||||
local w = `v'
|
||||
}
|
||||
if `w' < `max' {
|
||||
local max = `max' + `step'
|
||||
}
|
||||
if "`longitudinal'"!="" & `gvar'!=-1 {
|
||||
local max = min(`max',0.99999) /* pour avoir des corrélations < 1 */
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
DEFINITION DE LA METHODE
|
||||
*******************************************************************************/
|
||||
|
||||
if substr("`method'",1,3)=="pop" | substr("`method'",1,3)=="POP"{
|
||||
local method POPULATION+GH
|
||||
}
|
||||
|
||||
if "`method'"=="MEAN+GH"&`nbpatterns'*(`n1'+`n0')>=`nbpatmax' {
|
||||
di in gr "The MEAN+GH will be inefficient compared to GH since the maximal number of pattern's responses"
|
||||
di in gr "is lesser than the number of pattern retained by the MEAN+GH method."
|
||||
di in gr "The -method- option is replaced by GH."
|
||||
local method GH
|
||||
}
|
||||
else if ("`method'"=="MEAN+GH" | "`method'"=="MEAN") & `nbpatmax'>1000000{
|
||||
di in red "The number of patterns is too high for the chosen method"
|
||||
exit
|
||||
}
|
||||
else if "`method'"=="" {
|
||||
if `nbpatmax'<2000 {
|
||||
local method GH
|
||||
}
|
||||
else {
|
||||
local method POPULATION+GH
|
||||
}
|
||||
}
|
||||
|
||||
if "`method'"!="MEAN+GH" & "`method'"!="MEAN" & "`method'"!="GH" & "`method'"!="POPULATION+GH"{
|
||||
di in red "Invalid method name"
|
||||
exit
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
AFFICHAGE DES PARAMETRES
|
||||
*******************************************************************************/
|
||||
|
||||
if "`longitudinal'"==""{
|
||||
di in gr "Number of individuals in the first group: " in ye `n0'
|
||||
di in gr "Number of individuals in the second group: " in ye `n1'
|
||||
di in green "Group effect: " in ye `gamma'
|
||||
di in gr "Variance of the latent trait in the first group: " in ye `var0'
|
||||
di in gr "Variance of the latent trait in the second group: " in ye `var1'
|
||||
}
|
||||
else{
|
||||
di in gr "Number of individuals at each time: " in ye `n0'
|
||||
di in green "Time effect: " in ye `gamma'
|
||||
di in gr "Variance matrix of the latent trait: "
|
||||
matrix list `matvar',noheader
|
||||
di
|
||||
}
|
||||
di in gr "Number of items: " in ye `nbitems'
|
||||
di in green "Difficulties parameters of the items: "
|
||||
tempname dd
|
||||
matrix `dd'=`d''
|
||||
local items
|
||||
forvalues i=1/`nbitems' {
|
||||
local items "`items' item`i'"
|
||||
}
|
||||
local modalities
|
||||
forvalues i=1/`=`nbmodat'-1' {
|
||||
local modalities "`modalities' delta_`i'"
|
||||
}
|
||||
|
||||
matrix colnames `dd'=`items'
|
||||
matrix rownames `dd'=`modalities'
|
||||
|
||||
|
||||
di
|
||||
tempname ddt
|
||||
matrix `ddt' = `dd''
|
||||
matrix colnames `ddt'=`modalities'
|
||||
matrix rownames `ddt'=`items'
|
||||
if "`html'" != "" {
|
||||
matrix list `ddt',noblank nohalf noheader /* affiche la matrice des difficultés dans l'autre sens */
|
||||
}
|
||||
else {
|
||||
matrix list `dd',noblank nohalf noheader
|
||||
}
|
||||
di
|
||||
|
||||
if "`detail'"!="" & "`html'" == ""{
|
||||
di in gr "Method: " in ye "`method'"
|
||||
di in gr "Number of studied response's patterns: " in ye `nbpatmax'
|
||||
}
|
||||
matrix `dd'=`d'
|
||||
local gamma=`gamma'
|
||||
|
||||
local tmp=1
|
||||
qui matrix `db'=J(`=`nbitems'*(`nbmodat'-1)',1,.)
|
||||
forvalues j=1/`nbitems' {
|
||||
forvalues m=1/`=`nbmodat'-1' {
|
||||
qui matrix `db'[`tmp',1]=`d'[`j',`m']
|
||||
local ++tmp
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/*******************************************************************************
|
||||
CREATION DU DATASET ATTENDU
|
||||
*******************************************************************************/
|
||||
|
||||
if "`data'"=="" {
|
||||
clear
|
||||
if "`method'"!="POPULATION+GH"{
|
||||
local temp=`nbmodat'^(`nbtotitems')
|
||||
qui range x 0 `=`temp'-1' `temp'
|
||||
qui g t=x
|
||||
loc i=`nbtotitems'
|
||||
qui count if t>0
|
||||
loc z=r(N)
|
||||
qui while `z'>0 {
|
||||
qui gen item`=`nbtotitems'-`i'+1'=floor(t/`nbmodat'^`=`i'-1')
|
||||
qui replace t=mod(t,`nbmodat'^`=`i'-1')
|
||||
qui count if t>0
|
||||
loc z=r(N)
|
||||
loc i=`i'-1
|
||||
}
|
||||
drop t
|
||||
if "`longitudinal'"==""{
|
||||
qui expand 2
|
||||
qui gen group=0 in 1/`temp'
|
||||
qui replace group=1 in `=`temp'+1'/`=2*`temp''
|
||||
}
|
||||
}
|
||||
else { /*METHODE POPULATION (SIMULATION)*/
|
||||
if "`longitudinal'"==""{
|
||||
*qui simirt, clear pcm(`d') cov(`var') group(`=`n1'*`gamma'/(`n1'+`n0')') deltagroup(`gamma') nbobs(1000000)
|
||||
qui simirt, clear pcm(`d') cov(`var1') mu(`=`n0'*`gamma'/(`n1'+`n0')') nbobs(`=round(`n1'/(`n1'+`n0')*1000000)')
|
||||
qui gen group=1
|
||||
tempfile save1
|
||||
qui save `save1',replace
|
||||
qui simirt, clear pcm(`d') cov(`var0') mu(`=-`n1'*`gamma'/(`n1'+`n0')') nbobs(`=round(`n0'/(`n1'+`n0')*1000000)')
|
||||
qui gen group=0
|
||||
qui append using `save1'
|
||||
*qui save d:\tmp,replace
|
||||
bysort group:su lt
|
||||
qui drop lt1
|
||||
qui contract item* group, freq(freq)
|
||||
qui count
|
||||
local patobs=r(N)
|
||||
if `patobs'>=`=(`n0'+`n1')*`nbpatterns''{
|
||||
qui gen keep=0
|
||||
qui gsort +group -freq
|
||||
qui replace keep=1 in 1/`=`nbpatterns'*`n0''
|
||||
qui gsort -group -freq
|
||||
qui replace keep=1 in 1/`=`nbpatterns'*`n1''
|
||||
}
|
||||
else{
|
||||
qui gen keep=1
|
||||
}
|
||||
}
|
||||
else{
|
||||
qui simirt, clear pcm(`dlong') covmat(`matvar') mu(`mean1' `mean2') dim(`nbitems' `nbitems') nbobs(1000000)
|
||||
forvalues j=1/`nbitems'{
|
||||
rename itemA`j' item`j'
|
||||
rename itemB`j' item`=`nbitems'+`j''
|
||||
}
|
||||
qui drop lt1 lt2
|
||||
qui contract item*, freq(freq)
|
||||
local patobs=r(N)
|
||||
if `patobs'>=`=`n0'*`nbpatterns''{
|
||||
qui gen keep=0
|
||||
qui gsort -freq
|
||||
qui replace keep=1 in 1/`=`nbpatterns'*`n0''
|
||||
}
|
||||
else{
|
||||
qui gen keep=1
|
||||
}
|
||||
}
|
||||
qui keep if keep==1
|
||||
qui count
|
||||
local retain=r(N)
|
||||
di "Number of kept patterns:`retain'"
|
||||
local method GH
|
||||
}
|
||||
if "`longitudinal'"==""{
|
||||
qui gen mean1=-`n1'*`gamma'/(`n0'+`n1') if group==0
|
||||
qui replace mean1=`n0'*`gamma'/(`n0'+`n1') if group==1
|
||||
qui gen var1=`var0' if group==0
|
||||
qui replace var1=`var1' if group==1
|
||||
}
|
||||
else{
|
||||
qui gen mean1=`mean1'
|
||||
qui gen mean2=`mean2'
|
||||
}
|
||||
/*CALCUL DES PROBAS PAR PATTERNS DE REPONSE*/
|
||||
if "`method'"=="GH" {
|
||||
local temp=`nbmodat'^(`nbtotitems')
|
||||
local diff0=0
|
||||
qui gen proba=.
|
||||
local dixj=10
|
||||
qui count
|
||||
local tmp=r(N)
|
||||
|
||||
forvalues i=1/`tmp' {
|
||||
local dix=floor(`tmp'/10)
|
||||
if mod(`i',`dix')==0 & "`html'"=="" {
|
||||
qui su proba
|
||||
local tmpprob=r(sum)
|
||||
if "`dixj'"!="10" {
|
||||
di ".." _c
|
||||
}
|
||||
*di "`dixj'% (`tmpprob')" _c
|
||||
di "`dixj'%" _c
|
||||
local dixj=`dixj'+10
|
||||
}
|
||||
forvalues t=1/`nbt'{
|
||||
local int`t'=1
|
||||
forvalues j=`=(`t'-1)*`nbitems'+1'/`=`t'*`nbitems'' {
|
||||
qui su item`j' in `i'
|
||||
local rep=r(mean)
|
||||
local diff0=0
|
||||
local diff1=`dlong'[`j',1]
|
||||
local sum "1+exp(x`t'-`diff1')"
|
||||
forvalues m=2/`=`nbmodat'-1' {
|
||||
local diff`m'=`diff`=`m'-1''+`dlong'[`j',`m']
|
||||
local sum "`sum'+exp(`m'*x`t'-`diff`m'')"
|
||||
}
|
||||
local int`t' "(`int`t''*exp(`rep'*x`t'-`diff`rep''))/(`sum')"
|
||||
}
|
||||
}
|
||||
if "`longitudinal'"==""{
|
||||
qui su mean1 in `i'
|
||||
local mean=r(mean)
|
||||
qui su var1 in `i'
|
||||
local vart=r(mean)
|
||||
*di "gausshermite `int1',mu(`mean') sigma(`=sqrt(`vart')') display name(x1)"
|
||||
*prout
|
||||
qui gausshermite `int1',mu(`mean') sigma(`=sqrt(`vart')') display name(x1)
|
||||
qui replace proba=r(int) in `i'
|
||||
}
|
||||
else{
|
||||
local int "`int1'*`int2'"
|
||||
matrix mean=(`mean1',`mean2')
|
||||
qui gausshermite `int',mu(mean) var(`matvar') display name(x1 x2)
|
||||
qui replace proba=r(int) in `i'
|
||||
}
|
||||
}
|
||||
di
|
||||
}
|
||||
else { /*if "`method'"!="GH"*/
|
||||
qui gen proba=1
|
||||
forvalues t=1/`nbt'{
|
||||
forvalues i=`=(`t'-1)*`nbitems'+1'/`=`t'*`nbitems'' {
|
||||
local diff0=0
|
||||
local diff1=`dlong'[`i',1]
|
||||
qui gen eps0=1
|
||||
qui gen eps1=exp(mean`t'-`diff1')
|
||||
qui gen d=eps0+eps1
|
||||
forvalues m=2/`=`nbmodat'-1' {
|
||||
local diff`m'=`diff`=`m'-1''+`dlong'[`i',`m']
|
||||
qui gen eps`m'=exp(`m'*mean`t'-`diff`m'')
|
||||
qui replace d=d+eps`m'
|
||||
}
|
||||
local listeps
|
||||
forvalues m=0/`=`nbmodat'-1' {
|
||||
qui replace proba=proba*eps`m'/d if item`i'==`m'
|
||||
local listeps `listeps' eps`m'
|
||||
}
|
||||
qui drop `listeps' d
|
||||
}
|
||||
}
|
||||
if "`method'"=="MEAN+GH" {
|
||||
set tracedepth 1
|
||||
if "`longitudinal'"==""{
|
||||
qui gen keep=0
|
||||
qui gsort -group -proba
|
||||
local min=min(`=`nbmodat'^`nbitems'',`=`n1'*`nbpatterns'')
|
||||
qui replace keep=1 in 1/`min'
|
||||
qui gsort +group -proba
|
||||
local min=min(`=`nbmodat'^`nbitems'',`=`n0'*`nbpatterns'')
|
||||
qui replace keep=1 in 1/`min'
|
||||
qui keep if keep==1
|
||||
qui su proba if group==0
|
||||
local sumproba0=r(sum)*100
|
||||
qui su proba if group==1
|
||||
local sumproba1=r(sum)*100
|
||||
}
|
||||
else{
|
||||
qui gen keep=0
|
||||
qui gsort -proba
|
||||
local min=min(`nbpatmax',`=`n0'*`nbpatterns'')
|
||||
qui replace keep=1 in 1/`min'
|
||||
qui keep if keep==1
|
||||
}
|
||||
|
||||
qui drop keep proba
|
||||
local diff0=0
|
||||
qui gen proba=.
|
||||
qui count
|
||||
local nnew=r(N)
|
||||
di in gr "Number of studied response's patterns for the GH step: " in ye `nnew'
|
||||
if "`longitudinal'"==""{
|
||||
di in gr "(" in ye %6.2f `sumproba0' in gr "% of the group 0 and " in ye %6.2f `sumproba1' in gr "% of the group 1)"
|
||||
}
|
||||
local dixj=10
|
||||
forvalues i=1/`nnew' {
|
||||
local dix=floor(`nnew'/10)
|
||||
if mod(`i',`dix')==0 & "`html'"=="" {
|
||||
qui su proba
|
||||
local sum=r(sum)
|
||||
if "`dixj'"!="10" {
|
||||
di ".." _c
|
||||
}
|
||||
di "`dixj'% (`sum')" _c
|
||||
local dixj=`dixj'+10
|
||||
}
|
||||
forvalues t=1/`nbt'{
|
||||
local int`t'=1
|
||||
forvalues j=`=(`t'-1)*`nbitems'+1'/`=`t'*`nbitems'' {
|
||||
qui su item`j' in `i'
|
||||
local rep=r(mean)
|
||||
local diff0=0
|
||||
local diff1=`dlong'[`j',1]
|
||||
local sum "1+exp(x`t'-`diff1')"
|
||||
forvalues m=2/`=`nbmodat'-1' {
|
||||
local diff`m'=`diff`=`m'-1''+`dlong'[`j',`m']
|
||||
local sum "`sum'+exp(`m'*x`t'-`diff`m'')"
|
||||
}
|
||||
local int`t' "(`int`t''*exp(`rep'*x`t'-`diff`rep''))/(`sum')"
|
||||
}
|
||||
}
|
||||
if "`longitudinal'"==""{
|
||||
qui su mean1 in `i'
|
||||
local mean=r(mean)
|
||||
qui su var1 in `i'
|
||||
local vart=r(mean)
|
||||
*di "qui gausshermite `int1',mu(`mean') sigma(`=sqrt(`vart')') display name(x1)"
|
||||
qui gausshermite `int1',mu(`mean') sigma(`=sqrt(`vart')') display name(x1)
|
||||
qui replace proba=r(int) in `i'
|
||||
}
|
||||
else{
|
||||
local int "`int1'*`int2'"
|
||||
matrix mean=(`mean1',`mean2')
|
||||
qui gausshermite `int',mu(mean) var(`matvar') display
|
||||
qui replace proba=r(int) in `i'
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
qui gen eff=proba
|
||||
qui gen eff2=.
|
||||
|
||||
if "`longitudinal'"==""{
|
||||
qui keep item* eff* group proba
|
||||
local p1=1/`n1'
|
||||
local p0=1/`n0'
|
||||
qui replace eff2=floor(eff/`p1') if group==1
|
||||
qui replace eff2=floor(eff/`p0') if group==0
|
||||
qui replace eff=eff-eff2*(`p1'*group+`p0'*(1-group))
|
||||
qui su eff2 if group==1
|
||||
local aff1=r(sum)
|
||||
qui su eff2 if group==0
|
||||
local aff0=r(sum)
|
||||
|
||||
local unaff1=`n1'-`aff1'
|
||||
local unaff0=`n0'-`aff0'
|
||||
qui gen efftmp=eff2
|
||||
qui gsort + group - eff
|
||||
qui replace eff2=eff2+1 in 1/`unaff0'
|
||||
qui gsort - group - eff
|
||||
qui replace eff2=eff2+1 in 1/`unaff1'
|
||||
|
||||
qui drop if eff2==0
|
||||
qui gsort group item*
|
||||
qui expand eff2
|
||||
}
|
||||
else{
|
||||
qui keep item* eff* proba
|
||||
local p0=1/`n0'
|
||||
qui replace eff2=floor(eff/`p0')
|
||||
qui replace eff=eff-eff2*`p0'
|
||||
qui su eff2
|
||||
local aff0=r(sum)
|
||||
|
||||
local unaff0=`n0'-`aff0'
|
||||
qui gen efftmp=eff2
|
||||
qui gsort - eff
|
||||
qui replace eff2=eff2+1 in 1/`unaff0'
|
||||
|
||||
qui drop if eff2==0
|
||||
qui gsort item*
|
||||
qui expand eff2
|
||||
}
|
||||
qui drop proba eff eff2
|
||||
}
|
||||
/*
|
||||
if "`html'" == "" {
|
||||
qui save "d:\essai.dta", replace
|
||||
}*/
|
||||
qui alpha item*
|
||||
local alpha=r(alpha)
|
||||
|
||||
if "`longitudinal'"==""{ /*TRANSVERSAL*/
|
||||
qui gen groupc=-`n1'/(`n0'+`n1') if group==0
|
||||
qui replace groupc=`n0'/(`n0'+`n1') if group==1
|
||||
if `nbmodat'==2 { /*DICHOTOMOUS CASE*/
|
||||
if "`freeitems'"=="" {
|
||||
local offset "offset(offset)"
|
||||
local items
|
||||
}
|
||||
else {
|
||||
local offset
|
||||
local items "item1-item`nbitems'"
|
||||
}
|
||||
qui gen i=_n
|
||||
|
||||
tempname diff
|
||||
matrix `diff'=`dd''
|
||||
|
||||
qui reshape long item, i(i)
|
||||
qui rename item rep
|
||||
qui rename _j item
|
||||
|
||||
qui gen offset=0
|
||||
forvalues i=1/`nbitems' {
|
||||
qui replace offset=-`diff'[1,`i'] if item==`i'
|
||||
qui gen item`i'=item==`i'
|
||||
qui replace item`i'=-item`i'
|
||||
}
|
||||
if `var0'==`var1' { /*EQUAL VARIANCES*/
|
||||
*constraint 1 _cons=`=ln(`var')'
|
||||
*qui xtlogit rep groupc ,nocons i(i) offset(offset) constraint(1)
|
||||
if "`freevar'"=="" {
|
||||
constraint 1 _cons=`=sqrt(`var0')'
|
||||
local constvar "constraint(1)"
|
||||
}
|
||||
qui gllamm rep groupc `items',nocons i(i) `offset' `constvar' fam(bin) link(logit) nrf(1) trace iterate(100)
|
||||
}
|
||||
else {/*UNEQUAL VARIANCES*/
|
||||
tempvar G0 G1
|
||||
qui gen `G0'=group==0
|
||||
qui gen `G1'=group==1
|
||||
qui eq B0: `G0'
|
||||
qui eq B1: `G1'
|
||||
constraint 1 _cons=0
|
||||
if "`freevar'"=="" { /*FIXED UNEQUAL VARIANCE*/
|
||||
constraint 2 `G0'=`=sqrt(`var0')'
|
||||
constraint 3 `G1'=`=sqrt(`var1')'
|
||||
constraint 1 _cons=`=sqrt(`var0')'
|
||||
local constvar "constraint(1 2 3)"
|
||||
}
|
||||
else { /*FREE UNEQUAL VARIANCES*/
|
||||
local constvar "constraint(1)"
|
||||
}
|
||||
qui gllamm rep groupc `items',nocons i(i) `offset' `constvar' fam(bin) link(logit) eqs(B0 B1) nrf(2) trace iterate(100)
|
||||
}
|
||||
tempname b V/*LONGITUDINAL*/
|
||||
local row = 1
|
||||
}
|
||||
else { /*POLYTOMOUS CASE*/
|
||||
matrix `db'=`db''
|
||||
if "`freeitems'"=="" {
|
||||
local fi "fixed(`db') fixedmu"
|
||||
local row=1
|
||||
}
|
||||
else {
|
||||
local fi
|
||||
local row=(`nbmodat'-1)*`nbitems'+1
|
||||
}
|
||||
if `var0'==`var1' {
|
||||
if "`freevar'"=="" {
|
||||
local fv "fixedvar(`var0')"
|
||||
}
|
||||
else {
|
||||
local fv
|
||||
*noi pcm item*, `fi' covariates(groupc)
|
||||
|
||||
}
|
||||
}
|
||||
else {
|
||||
if "`freevar'"=="" {
|
||||
local fv "fixedvargroupc(`var')"
|
||||
}
|
||||
else {
|
||||
local fv "vargroupc"
|
||||
}
|
||||
}
|
||||
*di "jk pcm item*, `fi' covariates(groupc) `fv'"
|
||||
/********************************ESSAI 23 mai 2014*/
|
||||
qui simirt, clear pcm(`d') cov(`var1') mu(`=`n0'*`gamma'/(`n1'+`n0')') nbobs(`n1')
|
||||
qui gen group=1
|
||||
tempfile save1
|
||||
qui save `save1',replace
|
||||
qui simirt, clear pcm(`d') cov(`var0') mu(`=-`n1'*`gamma'/(`n1'+`n0')') nbobs(`n0')
|
||||
qui gen group=0
|
||||
qui append using `save1'
|
||||
qui gen groupc=group-0.5
|
||||
*di "ON A SIMULE"
|
||||
**********************************FIN ESSAI 23 mai 2014*/
|
||||
|
||||
qui pcmold item*, `fi' covariates(groupc) `fv'
|
||||
}
|
||||
tempname b V
|
||||
matrix `b'=e(b)
|
||||
matrix `V'=e(V)
|
||||
local gammaest=`b'[1,`row']
|
||||
local se=`V'[`row',`row']^.5
|
||||
}
|
||||
else{ /*LONGITUDINAL*/
|
||||
qui raschlong item*, nbt(`nbt') diff(`d') var(`matvar')
|
||||
tempname b V
|
||||
matrix `b'=e(b)
|
||||
matrix `V'=e(V)
|
||||
local gammaest=`b'[1,1]
|
||||
local se=`V'[1,1]^.5
|
||||
local n1=`n0'
|
||||
}
|
||||
local poweruni=1-normal(1.96-`gamma'/`se')
|
||||
|
||||
if "`html'" == "" {
|
||||
di
|
||||
di
|
||||
}
|
||||
if "`detail'"!=""{
|
||||
di in gr "{hline 91}"
|
||||
di _col(60) "Estimation with the "
|
||||
di _col(50) "Cramer-Rao bound" _col(75) "classical formula"
|
||||
di in gr "{hline 91}"
|
||||
}
|
||||
|
||||
if "`longitudinal'"==""{
|
||||
local power=1-normal(1.96-`gamma'/`se')+normal(-1.96-`gamma'/`se')
|
||||
local nbw:word count `var'
|
||||
if `nbw' == 2 { /* variance commune */
|
||||
local varc = ((`n0'-1)*`var0'+(`n1'-1)*`var1')/(`n0'+`n1'-2)
|
||||
}
|
||||
else {
|
||||
local varc = `var'
|
||||
}
|
||||
local clpower=normal(sqrt(`n1'*`gamma'^2/((`n1'/`n0'+1)*`varc'))-1.96)
|
||||
local clnsn=(`n1'/`n0'+1)/((`n1'/`n0')*(`gamma'/sqrt(`varc'))^2)*(1.96+invnorm(`poweruni'))^2
|
||||
local ratio=(`n0'+`n1')/(`clnsn'*(1+`n1'/`n0'))
|
||||
if "`detail'"==""{
|
||||
di in gr "{hline 65}"
|
||||
di in green "Estimation of the variance of the group effect" _col(56) in ye %10.4f `=`se'^2'
|
||||
di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni'
|
||||
di in gr "{hline 65}"
|
||||
}
|
||||
else{
|
||||
if "`gammafixed'"=="" {
|
||||
di in green "Estimated value of the group effect" _col(59) in ye %7.2f `gammaest'
|
||||
}
|
||||
di in green "Estimation of the s.e. of the group effect" _col(59) in ye %7.2f `se'
|
||||
di in green "Estimation of the variance of the group effect" _col(56) in ye %10.4f `=`se'^2'
|
||||
di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni' _col(86) in ye %6.4f `clpower'
|
||||
di in green "Number of patients for a power of" %6.2f `=`poweruni'*100' "%" _col(59) in ye `n0' "/" `n1' _col(77) in ye %7.2f `clnsn' "/" %7.2f `=`clnsn'*`n1'/`n0''
|
||||
di in green "Ratio of the number of patients" in ye %6.2f _col(68)`ratio'
|
||||
di in gr "{hline 91}"
|
||||
}
|
||||
}
|
||||
else{
|
||||
local clpower=normal(sqrt(`n0'*`gamma'^2)/(2*(`=`matvar'[1,1]'-`=`matvar'[2,1]'))-1.96)
|
||||
local clnsn=2*(`=`matvar'[1,1]'-`=`matvar'[2,1]')*(1.96+invnorm(`poweruni'))^2/(`gamma'^2)
|
||||
local ratio=(`n0'+`n1')/(`clnsn'*(1+`n1'/`n0'))
|
||||
if "`detail'"==""{
|
||||
di in gr "{hline 65}"
|
||||
di in green "Estimation of the variance of the group effect" _col(56) in ye %10.4f `=`se'^2'
|
||||
di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni'
|
||||
di in gr "{hline 65}"
|
||||
}
|
||||
else{
|
||||
if "`gammafixed'"=="" {
|
||||
di in green "Estimated value of the time effect" _col(59) in ye %7.2f `gammaest'
|
||||
}
|
||||
di in green "Estimation of the s.e. of the time effect" _col(59) in ye %7.2f `se'
|
||||
di in green "Estimation of the variance of the time effect" _col(56) in ye %10.4f `=`se'^2'
|
||||
di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni' _col(86) in ye %6.4f `clpower'
|
||||
di in green "Number of patients for a power of" %6.2f `=`poweruni'*100' "%" _col(63) in ye `n0' _col(85) in ye %7.2f `clnsn'
|
||||
di in green "Ratio of the number of patients" in ye %6.2f _col(68)`ratio'
|
||||
di in gr "{hline 91}"
|
||||
}
|
||||
}
|
||||
|
||||
if `expectedpower'!=-1 {
|
||||
qui sampsi `=-`gamma'/2' `=`gamma'/2', sd1(`=sqrt(`var')') sd2(`=sqrt(`var')') alpha(0.05) power(`expectedpower') ratio(`=`n1'/`n0'')
|
||||
local expn_1=r(N_1)
|
||||
local expn_2=r(N_2)
|
||||
local expn2=`expn_1'*`ratio'
|
||||
di in green "Number of patients for a power of" %6.2f `=`expectedpower'*100' "%" _col(51) in ye %7.2f `expn2' "/" %7.2f `=`expn2'*`n1'/`n0'' _c
|
||||
if "`detail'"!="" {
|
||||
di _col(77) in ye %7.2f `expn_1' "/" %7.2f `expn_2'
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
return scalar EstGamma=`gammaest'
|
||||
return scalar CRbound=`=`se'^2'
|
||||
return scalar CRPower=`poweruni'
|
||||
return scalar ClPower=`clpower'
|
||||
return scalar ClSS=`clnsn'
|
||||
return scalar Ratio=`ratio'
|
||||
return scalar CronbachAlpha=`alpha'
|
||||
|
||||
|
||||
if "`graph'" != "" {
|
||||
local x = "`min'(`step')`max'"
|
||||
matrix res=J(15,5,.)
|
||||
matrix colnames res= n0 n1 gamma variance power
|
||||
local l=1
|
||||
|
||||
if "`longitudinal'"==""{
|
||||
forval var = `x' {
|
||||
qui raschpower, n0(`n0') n1(`n1') gamma(`gamma') difficulties(`d') var(`var') html(`html')
|
||||
qui matrix res[`l',1]=(`n0', `n1', `gamma',`var',r(CRPower))
|
||||
local ++l
|
||||
}
|
||||
clear
|
||||
qui svmat res,names(col)
|
||||
|
||||
if "`html'" != "" {
|
||||
qui local saving "saving(`c(tmpdir)'/`html'_planif_graph,replace) nodraw"
|
||||
qui graph twoway (connected power variance,sort), title("n0=`n0', n1=`n1', {&gamma}=`gamma'") ylabel(0(0.2)1) xlabel(`x') xtitle(Variance of the latent variable, margin(top)) ytitle(Power) name(planif_graph,replace) xscale(range(`min' `=`max'+`step'')) `saving'
|
||||
qui graph use `c(tmpdir)'/`html'_planif_graph.gph
|
||||
qui graph export `c(tmpdir)'/`html'_planif_graph.eps, replace
|
||||
di "<br />"
|
||||
di "<img src=" _char(34) "/data/`html'_planif_graph.png" _char(34)
|
||||
di " class=" _char(34) "resgraph" _char(34) " alt=" _char(34) "planif_graph" _char(34) " title= " _char(34) "Power as a function of the variance of the latent variable - click to enlarge" _char(34) " width=" _char(34) "350" _char(34) " height=" _char(34) "240" _char(34) " >"
|
||||
}
|
||||
else {
|
||||
if "`filesave'"!="" {
|
||||
qui local saving "saving(`dirsave'//planif,replace) nodraw"
|
||||
}
|
||||
else {
|
||||
qui local saving
|
||||
}
|
||||
twoway (connected power variance,sort), title("n0=`n0', n1=`n1', {&gamma}=`gamma'") ylabel(0(0.2)1) xlabel(`x') xtitle(Variance of the latent variable, margin(top)) ytitle(Power) name(planif_graph,replace) xscale(range(`min' `=`max'+`step'')) `saving'
|
||||
}
|
||||
}
|
||||
|
||||
else {
|
||||
if `gcorr' != -1 {
|
||||
local l=1
|
||||
forval var = `x' {
|
||||
local cov = `gcorr'*`var'
|
||||
local matv "`var',`cov'" " \ `cov',`var'"
|
||||
qui raschpower, n0(`n0') n1(`n1') gamma(`gamma') difficulties(`d') var("`matv'") html(`html') `longitudinal'
|
||||
qui matrix res[`l',1]=(`n0', `n1', `gamma',`var',r(CRPower))
|
||||
local ++l
|
||||
}
|
||||
clear
|
||||
qui svmat res,names(col)
|
||||
if "`html'" != "" {
|
||||
qui local saving "saving(`c(tmpdir)'/`html'_planif_graph,replace) nodraw"
|
||||
qui graph twoway (connected power variance,sort), title("n=`n0', {&gamma}=`gamma', corr=`gcorr'") ylabel(0(0.2)1) xlabel(`x') xtitle(Variance of the latent variable, margin(top)) ytitle(Power) name(planif_graph,replace) `saving'
|
||||
qui graph use `c(tmpdir)'/`html'_planif_graph.gph
|
||||
qui graph export `c(tmpdir)'/`html'_planif_graph.eps, replace
|
||||
di "<br />"
|
||||
di "<img src=" _char(34) "/data/`html'_planif_graph.png" _char(34)
|
||||
di " class=" _char(34) "resgraph" _char(34) " alt=" _char(34) "planif_graph" _char(34) " title= " _char(34) "Power as a function of the variance of the latent variable - click to enlarge" _char(34) " width=" _char(34) "350" _char(34) " height=" _char(34) "240" _char(34) " >"
|
||||
}
|
||||
else {
|
||||
if "`filesave'"!="" {
|
||||
qui local saving "saving(`dirsave'//planif,replace) nodraw"
|
||||
}
|
||||
else {
|
||||
qui local saving
|
||||
}
|
||||
twoway (connected power variance,sort), title("n=`n0', {&gamma}=`gamma', corr=`gcorr'") ylabel(0(0.2)1) xlabel(`x') xtitle(Variance of the latent variable, margin(top)) ytitle(Power) name(planif_graph,replace) xscale(range(`min' `=`max'+`step'')) `saving'
|
||||
}
|
||||
}
|
||||
|
||||
if `gvar' != -1 {
|
||||
matrix res=J(9,5,.)
|
||||
matrix colnames res= n0 n1 gamma corr power
|
||||
local l=1
|
||||
forvalues corr = `x' {
|
||||
local cov = `corr'*`gvar'
|
||||
local matv "`gvar',`cov'" " \ `cov',`gvar'"
|
||||
qui raschpower, n0(`n0') n1(`n1') gamma(`gamma') difficulties(`d') var("`matv'") html(`html') `longitudinal'
|
||||
qui matrix res[`l',1]=(`n0', `n1', `gamma',`corr',r(CRPower))
|
||||
local ++l
|
||||
}
|
||||
clear
|
||||
qui svmat res,names(col)
|
||||
if "`html'" != "" {
|
||||
qui local saving "saving(`c(tmpdir)'/`html'_planif_graph,replace) nodraw"
|
||||
qui graph twoway (connected power corr,sort), title("n=`n0', {&gamma}=`gamma', var=`gvar'") ylabel(0(0.2)1) xlabel(`x') xtitle(Correlation, margin(top)) ytitle(Power) name(planif_graph,replace) `saving'
|
||||
qui graph use `c(tmpdir)'/`html'_planif_graph.gph
|
||||
qui graph export `c(tmpdir)'/`html'_planif_graph.eps, replace
|
||||
di "<br />"
|
||||
di "<img src=" _char(34) "/data/`html'_planif_graph.png" _char(34)
|
||||
di " class=" _char(34) "resgraph" _char(34) " alt=" _char(34) "planif_graph" _char(34) " title= " _char(34) "Power as a function of the variance of the latent variable - click to enlarge" _char(34) " width=" _char(34) "350" _char(34) " height=" _char(34) "240" _char(34) " >"
|
||||
}
|
||||
else {
|
||||
if "`filesave'"!="" {
|
||||
qui local saving "saving(`dirsave'//planif,replace) nodraw"
|
||||
}
|
||||
else {
|
||||
qui local saving
|
||||
}
|
||||
twoway (connected power corr,sort), title("n=`n0', {&gamma}=`gamma', var=`gvar'") ylabel(0(0.2)1) xlabel(`x') xtitle(Correlation) ytitle(Power) name(planif_graph,replace) xscale(range(`min' `=`max'+`step''))
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//di in gr "Method: " in ye "`method'"
|
||||
|
||||
if "`html'" != "" {
|
||||
di "</pre>"
|
||||
}
|
||||
}
|
||||
|
||||
capture qui use `raschpowerfile',clear
|
||||
|
||||
end
|
||||
@ -0,0 +1,78 @@
|
||||
{smcl}
|
||||
{* 6february2012}{...}
|
||||
{hline}
|
||||
help for {hi:raschpower}{right:Jean-Benoit Hardouin, Myriam Blanchin}
|
||||
{hline}
|
||||
|
||||
{title:Estimation of the power of the Wald test in order to compare the means of the latent trait in two groups of individuals}
|
||||
|
||||
{p 8 14 2}{cmd:raschpower} [{cmd:,} {cmdab:n0}({it:#}) {cmdab:n1}({it:#})
|
||||
{cmdab:gamma}({it:#}) {cmdab:var}({it:# [#]}) {cmdab:d:ifficulties}({it:matrix}) {cmdab:m:ethod}({it:method}) {cmdab:det:ail} {cmdab:exp:ectedpower}({it:#})]
|
||||
|
||||
|
||||
{title:Description}
|
||||
|
||||
{p 4 8 2}{cmd:raschpower} allows estimating the power of the Wald test comparing the means of two groups of patients in the context of the Rasch model or the partial-credit model. The estimation is based
|
||||
on the estimation of the variance of the difference of the means based on the Cramer-Rao bound.
|
||||
|
||||
{title:Options}
|
||||
|
||||
{p 4 8 2}{cmd:n0} and {cmd:n1} indicates the numbers of individuals in the two groups [100 by default].
|
||||
|
||||
{p 4 8 2}{cmd:gamma} indicates the group effect (difference between the two means) [0.5 by default].
|
||||
|
||||
{p 4 8 2}{cmd:var} is the expected values of the variances of the latent trait (1 by default): if this option contains only one value, variances are considered as equal between the two groups; if this option contains two values, variances are considered as unequal between the two groups.
|
||||
|
||||
{p 4 8 2}{cmdab:d:ifficulties} is a matrix containing the item parameters [one row per item, one column per positive modality - (-1.151, -0.987\-0.615, -0.325\-0.184, -0.043\0.246, 0.554\0.782, 1.724) by default].
|
||||
|
||||
{p 4 8 2}{cmd:method}({it:method}) indicates the method for constructing data. ({it:method}) may be GH, MEAN, MEAN+GH or POPULATION+GH [default is method(GH) if number of patterns<500, method(MEAN+GH) if 500<=number of patterns<10000,
|
||||
method(MEAN) if 10000<=number of patterns<1000000, method(POPULATION+GH) otherwise].
|
||||
|
||||
{p 8 14 2} {bf:GH}: The probability of all possible response patterns is estimated by Gauss-Hermite quadratures.
|
||||
|
||||
{p 8 14 2} {bf:MEAN}: The mean of the latent trait for each group is used instead of Gauss-Hermite quadratures.
|
||||
|
||||
{p 8 14 2} {bf:MEAN+GH}: In a first step, the MEAN method is used to determine the most probable patterns. In a second step, the probability of response patterns is estimated by Gauss-Hermite quadratures on the most probable patterns.
|
||||
|
||||
{p 8 14 2} {bf:POPULATION+GH}: The most frequent response patterns are selected from a simulated population of 1,000,000 of individuals. The probability of the selected response patterns is estimated by Gauss-Hermite quadratures.
|
||||
|
||||
{p 4 8 2}{cmd:detail} allows a comparison with the classical formula (for manifest variable).
|
||||
|
||||
{p 4 8 2}{cmd:expectedpower} allows searching for a sample size in order to reach a fixed level of power (the obtained sample sizes take into account the ratio between $n0$ and $n1$).
|
||||
|
||||
|
||||
{title:Example}
|
||||
|
||||
{p 4 8 2}{cmd:. raschpower}
|
||||
|
||||
{p 4 8 2}{cmd:. raschpower, n0(200) n1(200) gamma(0.4) var(1.3)}
|
||||
|
||||
{p 4 8 2}{cmd:. matrix diff=(-1.47\-0.97\-.23\-0.12\0.02\0.1)}{p_end}
|
||||
{p 4 8 2}{cmd:. raschpower, n0(127) n1(134) gamma(0.23) d(diff) var(2.58)}{p_end}
|
||||
|
||||
{p 4 8 2}{cmd:. matrix diff=(-0.328,-0.811,0.329\ 0.556,0.818,1.409\ 1.394,1.049,1.288\ 0.560,0.363,0.950)}{p_end}
|
||||
{p 4 8 2}{cmd:. raschpower, d(diff) n0(167) n1(205) gamma(0.178) var(0.77) exp(0.8)}{p_end}
|
||||
|
||||
|
||||
{title:References}
|
||||
|
||||
{p 4 8 2}Hardouin J.B., Amri S., Feddag M., Sébille V. (2012) Towards Power And Sample Size Calculations For The Comparison Of Two Groups Of Patients With Item Response Theory Models. Statistics in Medicine, 31(11): 1277-1290.{p_end}
|
||||
{p 4 8 2}Blanchin M., Hardouin J.B., Guillemin F., Falissard B., Sébille V. (2013) Power and sample size determination for the group comparison of patient-reported outcomes with Rasch family models. PLoS ONE, 8(2): e57279.{p_end}
|
||||
{p 4 8 2}Feddag M.L., Blanchin M., Hardouin J.B., Sébille V. (2014) Power analysis on the time effect for the longitudinal Rasch model. Journal of Applied Measurement: (under press).{p_end}
|
||||
{p 4 8 2}Guilleux A., Blanchin M., Hardouin J.B., Sébille V. (2014) Power and sample size determination in the Rasch model: Evaluation of the robustness of a numerical method to non-normality of the latent trait. Plos One, 9(1): e0083652.{p_end}
|
||||
|
||||
|
||||
{title:Author}
|
||||
|
||||
{p 4 8 2}Jean-Benoit Hardouin, PhD, assistant professor{p_end}
|
||||
{p 4 8 2}Myriam Blanchin, PhD, research assistant{p_end}
|
||||
{p 4 8 2}Bastien Perrot, biostatistician{p_end}
|
||||
{p 4 8 2}EA4275 "Biostatistics, Pharmacoepidemiology 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}Emails:
|
||||
{browse "mailto:jean-benoit.hardouin@univ-nantes.fr":jean-benoit.hardouin@univ-nantes.fr}{p_end}
|
||||
{p 13 8 2}{browse "mailto:myriam.blanchin@univ-nantes.fr":myriam.blanchin@univ-nantes.fr}{p_end}
|
||||
{p 4 8 2}Websites {browse "http://www.anaqol.org":AnaQol}
|
||||
and {browse "http://www.freeirt.org":FreeIRT}
|
||||
@ -0,0 +1,405 @@
|
||||
*! version 2.0.3 SRH 1 June 2003
|
||||
program define remcor
|
||||
version 6.0
|
||||
* takes b and transform set of r.effs at same level using choleski
|
||||
* returns bi where the r.effs are independent and correlation equations
|
||||
* have been removed
|
||||
* also takes exponential of sd parameters
|
||||
* and evaluates contributions to linear predictor that don't change
|
||||
args b stop
|
||||
|
||||
tempname b2 s1 cov t r d dd u denom mean mzps
|
||||
tempvar junk
|
||||
gen double `junk'=0
|
||||
global HG_error=0
|
||||
|
||||
disp "*********in remcor:"
|
||||
|
||||
/* fixed effects $HG_xb1, $HG_xb2 etc. (tempnames stored in global macros-list) */
|
||||
qui disp "fixed parameters: "
|
||||
qui matrix list `b'
|
||||
if $HG_const==1{
|
||||
matrix `b2' = `b'*M_T' + M_a
|
||||
matrix coleq `b2' = $HG_cole
|
||||
matrix colnames `b2' = $HG_coln
|
||||
noi matrix list `b2'
|
||||
}
|
||||
else{
|
||||
matrix `b2' = `b'
|
||||
}
|
||||
local nffold=0
|
||||
local ff = 1
|
||||
local nxt = 1
|
||||
while(`ff' <= $HG_tpff){
|
||||
matrix `b2' = `b2'[1, `nxt'...]
|
||||
local np = M_nffc[1, `ff'] - `nffold'
|
||||
qui disp "np = " `np'
|
||||
if `np'>0{
|
||||
local nxt = `np' + 1
|
||||
local nffold = M_nffc[1,`ff']
|
||||
matrix `s1' = `b2'[1,1..`np']
|
||||
qui matrix list `s1'
|
||||
qui disp "tempname: ${HG_xb`ff'}"
|
||||
matrix score double ${HG_xb`ff'} = `s1' /* nontemp */
|
||||
}
|
||||
else{
|
||||
qui gen double ${HG_xb`ff'} = 0
|
||||
}
|
||||
qui disp ${HG_xb`ff'}[$which]
|
||||
local ff=`ff'+1
|
||||
}
|
||||
if "$HG_off"~=""{qui replace $HG_xb1=$HG_xb1+$HG_off}
|
||||
qui disp "HG_xb1 = " $HG_xb1[$which]
|
||||
|
||||
if $HG_ethr{
|
||||
local j = 1
|
||||
local ii = 1
|
||||
while `ii'<=$HG_nolog{
|
||||
local j = `j' + 1
|
||||
local jm = `j' + M_nresp[1,`ii']-3
|
||||
while `j' <= `jm'{
|
||||
local jp = `j' + 1
|
||||
* disp in re "replace HG_xb`jp' = HG_xb`j' + exp(HG_xb`jp')"
|
||||
qui replace ${HG_xb`jp'} = ${HG_xb`j'} + exp(${HG_xb`jp'})
|
||||
local j = `j' + 1
|
||||
}
|
||||
local ii = `ii' + 1
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* random effects */
|
||||
/* level 1 */
|
||||
local np = M_nbrf[1,1]
|
||||
if `np'>0{
|
||||
matrix `b2' = `b2'[1, `nxt'...]
|
||||
local nxt = 1
|
||||
matrix `s1' = `b2'[1,1..`np']
|
||||
qui matrix list `s1'
|
||||
matrix score double $HG_s1 = `s1'
|
||||
if $HG_nats{
|
||||
qui replace $HG_s1 = abs($HG_s1)
|
||||
}
|
||||
else{
|
||||
qui replace $HG_s1=exp($HG_s1)
|
||||
}
|
||||
qui disp "s1 = $HG_s1 = " $HG_s1[$which]
|
||||
local nxt = `nxt' + `np'
|
||||
}
|
||||
local lev = 2
|
||||
local rf = 2
|
||||
local nrfold = M_nrfc[2,1]
|
||||
|
||||
/* MASS POINTS */
|
||||
if($HG_free){
|
||||
tempname pdenom
|
||||
gen double `pdenom' = 1.0
|
||||
while(`lev'<=$HG_tplv&`rf'<=$HG_tprf){
|
||||
local j1 = M_nrfc[2, `lev']
|
||||
local nrf = `j1' - `nrfold'
|
||||
local nip = M_nip[1, `lev']
|
||||
scalar `denom' = 1 /* =exp(0) */
|
||||
qui replace `pdenom' = 1.0
|
||||
matrix M_zps`lev' = J(1,`nip',0)
|
||||
local k = 1
|
||||
qui disp "`nip' integration points at level `lev'"
|
||||
local nloc = `nip'
|
||||
local npar = M_np[1,`lev']
|
||||
if $HG_cip|`nip'>1{ local nloc = `nloc'-1}
|
||||
while `k' <= `nloc' {
|
||||
local j = `nrfold'+1
|
||||
while `j'<=`j1'{
|
||||
qui disp "level `lev', class `k' and random effect `j'"
|
||||
qui disp " nxt = " `nxt'
|
||||
matrix `b2' = `b2'[1, `nxt'...]
|
||||
qui matrix list `b2'
|
||||
local nxt = 1
|
||||
if `k'==1{
|
||||
/* linear predictors come before first masspoint */
|
||||
local np = M_nbrf[1,`j']-1
|
||||
if `np'>0 {
|
||||
qui disp "extracting coefficients for r.eff"
|
||||
matrix `s1' = `b2'[1,1..`np']
|
||||
matrix score double ${HG_s`j'} = `s1'
|
||||
qui disp "HG_s`j' = ${HG_s`j'} = " ${HG_s`j'}[$which]
|
||||
local nxt = `nxt' + `np'
|
||||
}
|
||||
/* first coeff fixed at one */
|
||||
if M_frld[1,`rf']~=1{
|
||||
matrix `s1' = (1)
|
||||
local lab: colnames `b2'
|
||||
local lab: word `nxt' of `lab'
|
||||
matrix colnames `s1'=`lab'
|
||||
qui matrix list `s1'
|
||||
capture drop `junk'
|
||||
matrix score double `junk' = `s1'
|
||||
if `np'>0{
|
||||
qui replace ${HG_s`j'}=${HG_s`j'}+`junk'
|
||||
}
|
||||
else{
|
||||
qui gen double ${HG_s`j'}=`junk'
|
||||
}
|
||||
qui matrix list `s1'
|
||||
}
|
||||
qui disp "HG_s`j' = ${HG_s`j'} = " ${HG_s`j'}[$which]
|
||||
qui disp "making M_zlc`j'"
|
||||
matrix M_zlc`j' = J(1,`nip',0)
|
||||
}
|
||||
matrix M_zlc`j'[1,`k'] = `b2'[1,`nxt']
|
||||
local nxt = `nxt' + 1
|
||||
local j = `j' + 1
|
||||
}
|
||||
if `k'<`nip'{
|
||||
if `npar'>0{
|
||||
qui disp "extract probability parameters for HG_p`cl' (npar=`npar')"
|
||||
matrix `s1' = `b2'[1,`nxt'..(`nxt'+`npar'-1)]
|
||||
qui matrix list `s1'
|
||||
capture drop `junk'
|
||||
matrix score double `junk' = `s1'
|
||||
qui gen double ${HG_p`lev'`k'} = `junk'
|
||||
local nxt = `nxt' + `npar' - 1
|
||||
qui replace `pdenom' = `pdenom' + exp(${HG_p`lev'`k'})
|
||||
}
|
||||
|
||||
scalar `mzps' = exp(`b2'[1,`nxt'])
|
||||
if `mzps' == . {
|
||||
global HG_error=1
|
||||
exit
|
||||
}
|
||||
matrix M_zps`lev'[1,`k'] = `b2'[1,`nxt']
|
||||
local nxt = `nxt' + 1
|
||||
scalar `denom' = `denom' + `mzps'
|
||||
}
|
||||
local k = `k' + 1
|
||||
}
|
||||
|
||||
if `npar'>0{
|
||||
qui gen double ${HG_p`lev'`nip'} = 0.0
|
||||
local k = 1
|
||||
while `k' <= `nip'{
|
||||
qui disp "divide HG_p`lev'`k'"
|
||||
qui replace ${HG_p`lev'`k'} = ${HG_p`lev'`k'} - ln(`pdenom')
|
||||
local k = `k' + 1
|
||||
}
|
||||
}
|
||||
|
||||
local k = 1
|
||||
while `k' <= `nip'{
|
||||
matrix M_zps`lev'[1,`k'] = M_zps`lev'[1,`k'] - ln(`denom')
|
||||
local k = `k' + 1
|
||||
}
|
||||
local j = `nrfold' + 1
|
||||
while `j' <= `j1'{ /* define last location */
|
||||
if $HG_cip == 1{
|
||||
local k = 1
|
||||
scalar `mean' = 0
|
||||
while `k'<`nip'{
|
||||
scalar `mean' = `mean' + M_zlc`j'[1,`k']*exp(M_zps`lev'[1,`k'])
|
||||
local k = `k' + 1
|
||||
}
|
||||
scalar `mzps' = exp(M_zps`lev'[1,`nip'])
|
||||
matrix M_zlc`j'[1,`nip'] = -`mean'/`mzps'
|
||||
}
|
||||
else if `nip'>1{
|
||||
matrix M_zlc`j'[1,`nip'] = `b2'[1,`nxt']
|
||||
local nxt = `nxt' + 1
|
||||
}
|
||||
qui disp "M_zlc`j'"
|
||||
qui matrix list M_zlc`j'
|
||||
local j = `j' + 1
|
||||
}
|
||||
qui disp "M_zps`lev'"
|
||||
qui matrix list M_zps`lev'
|
||||
local nrfold = `j1'
|
||||
local lev = `lev' + 1
|
||||
}
|
||||
}/*endif HG_free */
|
||||
else{
|
||||
/* ST. DEVS */
|
||||
qui disp "random parameters: "
|
||||
if $HG_tprf>1{matrix CHmat = J($HG_tprf-1,$HG_tprf-1,0)}
|
||||
while(`lev'<=$HG_tplv&`rf'<=$HG_tprf){
|
||||
local np = M_nbrf[1,`rf']
|
||||
qui disp "np = " `np'
|
||||
local nrf = M_nrfc[2, `lev'] - `nrfold'
|
||||
matrix `t' = J(`nrf',`nrf',0)
|
||||
local i = 1
|
||||
while (`i' <= `nrf'){
|
||||
qui disp " nxt = " `nxt'
|
||||
matrix `b2' = `b2'[1, `nxt'...]
|
||||
local nxt = 1
|
||||
qui matrix list `b2'
|
||||
local np = M_nbrf[1, `rf'] - 1
|
||||
qui disp `np' " loadings at random effect " `rf' ", level " `lev'
|
||||
if `np'>0{
|
||||
matrix `s1' = `b2'[1,1..`np']
|
||||
/*
|
||||
* fudge: exponentiate s1
|
||||
local ij = 1
|
||||
while `ij'<=`np'{
|
||||
matrix `s1'[1,`ij'] = exp(`s1'[1,`ij'])
|
||||
local ij = `ij' + 1
|
||||
}
|
||||
* end fudge
|
||||
*/
|
||||
|
||||
qui matrix list `s1'
|
||||
matrix score double ${HG_s`rf'} = `s1'
|
||||
local nxt = `nxt' + `np'
|
||||
}
|
||||
/* first (single non-) loading fixed at one, label in st. dev */
|
||||
if M_frld[1,`rf']~=1{
|
||||
matrix `s1' = (1)
|
||||
local lab: colnames `b2'
|
||||
local lab: word `nxt' of `lab'
|
||||
matrix colnames `s1' = `lab'
|
||||
capture drop `junk'
|
||||
tempname junk
|
||||
matrix score double `junk' = `s1'
|
||||
if `np'>0{
|
||||
qui replace ${HG_s`rf'} = ${HG_s`rf'} + `junk'
|
||||
}
|
||||
else{
|
||||
matrix score double ${HG_s`rf'} = `s1'
|
||||
*qui replace ${HG_s`rf'} = `junk'
|
||||
}
|
||||
}
|
||||
|
||||
qui disp "HG_s`rf' = ${HG_s`rf'} = " ${HG_s`rf'}[$which]
|
||||
* extract standard deviation
|
||||
* fudge: take exponential
|
||||
* matrix `t'[`i',`i'] = exp(`b2'[1, `nxt'])
|
||||
matrix `t'[`i',`i'] = `b2'[1, `nxt']
|
||||
matrix CHmat[`rf'-1,`rf'-1]=`t'[`i',`i']
|
||||
local nxt = `nxt' + 1
|
||||
local i = `i' + 1
|
||||
local rf = `rf' + 1
|
||||
}
|
||||
if (`nrf'>1&$HG_cor==1){ /* deal with correlations */
|
||||
/* extract correlation parameters */
|
||||
local i = 2
|
||||
while (`i' <= `nrf'){
|
||||
local k = `i' + `nrfold' - 1
|
||||
local j = 1
|
||||
while (`j' < `i'){
|
||||
local l = `j' + `nrfold' - 1
|
||||
qui disp "i = " `i' " j = " `j' " nxt = " `nxt'
|
||||
matrix `t'[`i',`j'] = `b2'[1,`nxt']
|
||||
matrix CHmat[`k',`l'] = `t'[`i',`j']
|
||||
local j = `j' + 1
|
||||
local nxt = `nxt' + 1
|
||||
}
|
||||
local i = `i' + 1
|
||||
}
|
||||
}
|
||||
qui matrix list `t'
|
||||
matrix M_chol = `t'
|
||||
/* unpacked parameters */
|
||||
local nrfold = M_nrfc[2,`lev']
|
||||
local lev = `lev' + 1
|
||||
} /* loop through levels */
|
||||
}/*endelse HG_free */
|
||||
if "`stop'"~=""{
|
||||
exit
|
||||
}
|
||||
local nrfold = M_nrfc[2,1]
|
||||
/* use B-matrix */
|
||||
if $HG_tprf>1&$HG_bmat==1{
|
||||
qui disp "dealing with B-matrix"
|
||||
local i = 1
|
||||
matrix Bmat = J($HG_tprf-1,$HG_tprf-1,0)
|
||||
while `i'<$HG_tprf{
|
||||
local j = 1
|
||||
while `j' < $HG_tprf{
|
||||
if M_b[`i',`j']>0{
|
||||
matrix Bmat[`i',`j']=`b2'[1,`nxt']
|
||||
local nxt = `nxt' + 1
|
||||
}
|
||||
local j = `j' + 1
|
||||
}
|
||||
local i = `i' + 1
|
||||
}
|
||||
qui matrix list Bmat
|
||||
|
||||
/* only works if B-matrix is upper diagonal */
|
||||
local i=2
|
||||
while `i'<$HG_tprf{
|
||||
local k = `i' + `nrfold'
|
||||
local j = 1
|
||||
qui disp "making s`k'"
|
||||
while `j'<`i'{
|
||||
local l = `j' + `nrfold'
|
||||
qui replace ${HG_s`k'} = ${HG_s`k'} + Bmat[`j',`i']*${HG_s`l'}
|
||||
qui disp " adding Bmat[`j',`i']s`l'"
|
||||
local j = `j' + 1
|
||||
}
|
||||
local i = `i' + 1
|
||||
}
|
||||
}
|
||||
|
||||
/* deal with geqs */
|
||||
|
||||
|
||||
if $HG_ngeqs>0{
|
||||
qui disp "dealing with geqs"
|
||||
local i = 1
|
||||
while `i'<=$HG_ngeqs{
|
||||
local k = M_ngeqs[1,`i']
|
||||
local n = M_ngeqs[2,`i']
|
||||
qui disp "random effect `k' has `n' covariates"
|
||||
local nxt2 = `nxt'+`n'-1
|
||||
matrix `s1' = `b2'[1,`nxt'..`nxt2']
|
||||
qui matrix list `s1'
|
||||
local nxt = `nxt2' + 1
|
||||
capture drop `junk'
|
||||
matrix score double `junk' = `s1'
|
||||
qui disp "multiply " `junk'[$which] " by HG_s`k' and add to HG_xb1"
|
||||
qui replace $HG_xb1 = $HG_xb1 + `junk'*${HG_s`k'}
|
||||
qui disp "HG_xb1:" $HG_xb1[$which]
|
||||
local i = `i' + 1
|
||||
}
|
||||
}
|
||||
|
||||
/* use inter */
|
||||
|
||||
if $HG_inter~=0{
|
||||
local k = $HG_l + 1
|
||||
local j = $HG_r + 1
|
||||
qui disp "HG_s`k' = HG_s`k'*HG_s`j'
|
||||
qui replace ${HG_s`k'} = ${HG_s`k'}*${HG_s`j'}
|
||||
}
|
||||
|
||||
/* use CHmat */
|
||||
if $HG_free==0&$HG_tprf>1{
|
||||
qui disp "dealing with Cholesky matrix"
|
||||
qui matrix list CHmat
|
||||
local i = 1
|
||||
while (`i'<$HG_tprf){
|
||||
local k = `i' + `nrfold'
|
||||
qui replace `junk'=0
|
||||
local j = `i'
|
||||
qui disp "making s`k'"
|
||||
while `j'<$HG_tprf{
|
||||
local l = `j' + `nrfold'
|
||||
qui replace `junk' = `junk' + CHmat[`j',`i']*${HG_s`l'}
|
||||
qui disp " adding CHmat[`j',`i']s`l'"
|
||||
local j = `j' + 1
|
||||
}
|
||||
qui replace ${HG_s`k'}=`junk'
|
||||
qui disp "s`k' = ${HG_s`k'} = " ${HG_s`k'}[$which]
|
||||
local i = `i' + 1
|
||||
}
|
||||
}
|
||||
|
||||
* label M_znow
|
||||
local i=2
|
||||
local lab
|
||||
while `i'<=$HG_tprf{
|
||||
local lab "`lab' ${HG_s`i'}"
|
||||
local i = `i' + 1
|
||||
}
|
||||
matrix colnames M_znow=`lab'
|
||||
qui disp "M_znow:"
|
||||
qui matrix list M_znow
|
||||
end
|
||||
|
||||
Binary file not shown.
Loading…
Reference in New Issue