You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
372 lines
13 KiB
Plaintext
372 lines
13 KiB
Plaintext
8 months ago
|
*! version 3.3 : October 2nd, 2012
|
||
|
|
*! Jean-Benoit Hardouin, Myriam Blanchin
|
||
|
|
************************************************************************************************************
|
||
|
|
* 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 25, 2010 (Jean-Benoit Hardouin)
|
||
|
|
* Version 1.1 : January 26, 2010 (Jean-Benoit Hardouin)
|
||
|
|
* Version 1.2 : November 1st, 2010 (Jean-Benoit Hardouin)
|
||
|
|
* version 1.3 : May 2th, 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 3.3 : October 2nd, 2012 (Jean-Benoit Hardouin, Myriam Blanchin) : minor corrections
|
||
|
|
*
|
||
|
|
* Jean-benoit Hardouin, jean-benoit.hardouin@univ-nantes.fr
|
||
|
|
* Myriam Blanchin, myriam.blanchin@univ-nantes.fr
|
||
|
|
* EA 4275 "Biostatistics, Pharmacoepidemiology and Subjectives Measures in Health"
|
||
|
|
* Faculty of Pharmaceutical Sciences - University of Nantes - France
|
||
|
|
*
|
||
|
|
* News about this program : http://www.anaqol.org
|
||
|
|
* FreeIRT Project : http://www.freeirt.org
|
||
|
|
*
|
||
|
|
* Copyright 2010-2012 Jean-Benoit Hardouin, 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 raschpower,rclass
|
||
|
|
syntax [varlist] [, n0(int 100) n1(int 100) Gamma(real .5) Difficulties(string) Var(real 1) Method(string) NBPatterns(int 2) nodata EXPectedpower(real -1)]
|
||
|
|
version 11
|
||
|
|
|
||
|
|
tempfile raschpowerfile
|
||
|
|
capture qui save "`raschpowerfile'",replace
|
||
|
|
tempname db d
|
||
|
|
if "`difficulties'"=="" {
|
||
|
|
matrix `d'=[-1.151, -0.987\-0.615, -0.325\-0.184, -0.043\0.246, 0.554\0.782, 1.724]
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
matrix `d'=`difficulties'
|
||
|
|
}
|
||
|
|
local nbitems=rowsof(`d')
|
||
|
|
local nbmodat=colsof(`d')+1
|
||
|
|
if "`method'"=="MEAN+GH"&`nbpatterns'*(`n1'+`n0')>=`=`nbmodat'^`nbitems'*2' {
|
||
|
|
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'"=="" {
|
||
|
|
if `nbmodat'^`nbitems'*2<1000 {
|
||
|
|
local method GH
|
||
|
|
}
|
||
|
|
else if `nbmodat'^`nbitems'<10000 {
|
||
|
|
local method MEAN+GH
|
||
|
|
}
|
||
|
|
else if `nbmodat'^`nbitems'<1000000 {
|
||
|
|
local method MEAN
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
local method POPULATION+GH
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
di in gr "Method: " in ye "`method'"
|
||
|
|
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 ye `var'
|
||
|
|
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'
|
||
|
|
matrix list `dd',noblank nohalf noheader
|
||
|
|
di in gr "Number of studied response's patterns: " in ye `=`nbmodat'^`nbitems'*2'
|
||
|
|
|
||
|
|
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
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if "`data'"=="" {
|
||
|
|
clear
|
||
|
|
if "`method'"!="POPULATION+GH"{
|
||
|
|
local temp=`nbmodat'^(`nbitems')
|
||
|
|
qui range x 0 `=`temp'-1' `temp'
|
||
|
|
qui g t=x
|
||
|
|
loc i=`nbitems'
|
||
|
|
qui count if t>0
|
||
|
|
loc z=r(N)
|
||
|
|
qui while `z'>0 {
|
||
|
|
qui gen item`=`nbitems'-`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
|
||
|
|
qui expand 2
|
||
|
|
qui gen group=0 in 1/`temp'
|
||
|
|
qui replace group=1 in `=`temp'+1'/`=2*`temp''
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
qui simirt, clear pcm(`difficulties') cov(`var') group(`=`n1'/(`n1'+`n0')') deltagroup(`gamma') nbobs(1000000)
|
||
|
|
qui drop lt1
|
||
|
|
qui contract item* group, freq(freq)
|
||
|
|
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''
|
||
|
|
qui keep if keep==1
|
||
|
|
qui count
|
||
|
|
local tmp=r(N)
|
||
|
|
di "Number of kept patterns:`tmp'"
|
||
|
|
local method GH
|
||
|
|
}
|
||
|
|
qui gen mean=-`n1'*`gamma'/(`n0'+`n1') if group==0
|
||
|
|
qui replace mean=`n0'*`gamma'/(`n0'+`n1') if group==1
|
||
|
|
|
||
|
|
if "`method'"=="GH" {
|
||
|
|
local temp=`nbmodat'^(`nbitems')
|
||
|
|
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 {
|
||
|
|
if "`dixj'"!="10" {
|
||
|
|
di ".." _c
|
||
|
|
}
|
||
|
|
di "`dixj'%" _c
|
||
|
|
local dixj=`dixj'+10
|
||
|
|
}
|
||
|
|
local int=1
|
||
|
|
forvalues j=1/`nbitems' {
|
||
|
|
qui su item`j' in `i'
|
||
|
|
local rep=r(mean)
|
||
|
|
local diff0=0
|
||
|
|
local diff1=`d'[`j',1]
|
||
|
|
local sum "1+exp(x-`diff1')"
|
||
|
|
forvalues m=2/`=`nbmodat'-1' {
|
||
|
|
local diff`m'=`diff`=`m'-1''+`d'[`j',`m']
|
||
|
|
local sum "`sum'+exp(`m'*x-`diff`m'')"
|
||
|
|
}
|
||
|
|
local int "(`int'*exp(`rep'*x-`diff`rep''))/(`sum')"
|
||
|
|
}
|
||
|
|
qui su mean in `i'
|
||
|
|
local mean=r(mean)
|
||
|
|
qui gausshermite `int',mu(`mean') sigma(`=sqrt(`var')') display
|
||
|
|
qui replace proba=r(int) in `i'
|
||
|
|
}
|
||
|
|
di
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
qui gen proba=1
|
||
|
|
forvalues i=1/`nbitems' {
|
||
|
|
local diff0=0
|
||
|
|
local diff1=`d'[`i',1]
|
||
|
|
qui gen eps0=1
|
||
|
|
qui gen eps1=exp(mean-`diff1')
|
||
|
|
qui gen d=eps0+eps1
|
||
|
|
forvalues m=2/`=`nbmodat'-1' {
|
||
|
|
local diff`m'=`diff`=`m'-1''+`d'[`i',`m']
|
||
|
|
qui gen eps`m'=exp(`m'*mean-`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
|
||
|
|
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
|
||
|
|
|
||
|
|
|
||
|
|
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'
|
||
|
|
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 {
|
||
|
|
if "`dixj'"!="10" {
|
||
|
|
di ".." _c
|
||
|
|
}
|
||
|
|
di "`dixj'%" _c
|
||
|
|
local dixj=`dixj'+10
|
||
|
|
}
|
||
|
|
local int=1
|
||
|
|
forvalues j=1/`nbitems' {
|
||
|
|
qui su item`j' in `i'
|
||
|
|
local rep=r(mean)
|
||
|
|
local diff0=0
|
||
|
|
local diff1=`d'[`j',1]
|
||
|
|
local sum "1+exp(x-`diff1')"
|
||
|
|
forvalues m=2/`=`nbmodat'-1' {
|
||
|
|
local diff`m'=`diff`=`m'-1''+`d'[`j',`m']
|
||
|
|
local sum "`sum'+exp(`m'*x-`diff`m'')"
|
||
|
|
}
|
||
|
|
local int "(`int'*exp(`rep'*x-`diff`rep''))/(`sum')"
|
||
|
|
}
|
||
|
|
qui su mean in `i'
|
||
|
|
local mean=r(mean)
|
||
|
|
qui gausshermite `int',mu(`mean') sigma(`=sqrt(`var')') display
|
||
|
|
qui replace proba=r(int) in `i'
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
qui gen eff=.
|
||
|
|
forvalues i=0/1 {
|
||
|
|
qui replace eff=proba*`n`i'' if group==`i'
|
||
|
|
}
|
||
|
|
qui replace eff=proba
|
||
|
|
qui keep item* eff group proba
|
||
|
|
|
||
|
|
local p1=1/`n1'
|
||
|
|
local p0=1/`n0'
|
||
|
|
qui gen eff2=.
|
||
|
|
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
|
||
|
|
gsort group item*
|
||
|
|
qui expand eff2
|
||
|
|
qui drop proba eff eff2
|
||
|
|
}
|
||
|
|
|
||
|
|
qui alpha item*
|
||
|
|
local alpha=r(alpha)
|
||
|
|
qui gen groupc=group-.5
|
||
|
|
if `nbmodat'==2 {
|
||
|
|
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'
|
||
|
|
}
|
||
|
|
|
||
|
|
constraint 1 _cons=`=ln(`var')'
|
||
|
|
qui xtlogit rep groupc ,nocons i(i) offset(offset) constraint(1)
|
||
|
|
tempname b V
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
matrix `db'=`db''
|
||
|
|
*di "qui pcm item*, fixed(`db') covariates(groupc) fixedmu fixedvar(`var')"
|
||
|
|
qui pcm item*, fixed(`db') covariates(groupc) fixedmu fixedvar(`var')
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
tempname b V
|
||
|
|
matrix `b'=e(b)
|
||
|
|
matrix `V'=e(V)
|
||
|
|
local gammaest=`b'[1,1]
|
||
|
|
local se=`V'[1,1]^.5
|
||
|
|
|
||
|
|
|
||
|
|
di
|
||
|
|
di
|
||
|
|
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 "`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'
|
||
|
|
local power=1-normal(1.96-`gamma'/`se')+normal(-1.96-`gamma'/`se')
|
||
|
|
local poweruni=1-normal(1.96-`gamma'/`se')
|
||
|
|
local clpower=normal(sqrt(`n1'*`gamma'^2/((`n1'/`n0'+1)*`var'))-1.96)
|
||
|
|
di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni' _col(86) in ye %6.4f `clpower'
|
||
|
|
local clnsn=(`n1'/`n0'+1)/((`n1'/`n0')*(`gamma'/sqrt(`var'))^2)*(1.96+invnorm(`poweruni'))^2
|
||
|
|
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''
|
||
|
|
local ratio=(`n0'+`n1')/(`clnsn'*(1+`n1'/`n0'))
|
||
|
|
di in green "Ratio of the number of patients" in ye %6.2f _col(68)`ratio'
|
||
|
|
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'' _col(77) in ye %7.2f `expn_1' "/" %7.2f `expn_2'
|
||
|
|
}
|
||
|
|
di in gr "{hline 91}"
|
||
|
|
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'
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
capture qui use `raschpowerfile',clear
|
||
|
|
|
||
|
|
end
|