|
|
*! version 1.5 : July 11th, 2011
|
|
|
*! Jean-Benoit Hardouin
|
|
|
************************************************************************************************************
|
|
|
* 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
|
|
|
*
|
|
|
* Jean-benoit Hardouin, Faculty of Pharmaceutical Sciences - 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 2010-2011 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 raschpower,rclass
|
|
|
syntax [varlist] [, n0(int 100) n1(int 100) gamma(real .5) d(string) var(real 1) fast nodata gammafix EXPectedpower(real -1)]
|
|
|
version 11
|
|
|
|
|
|
tempfile raschpowerfile
|
|
|
capture qui save "`raschpowerfile'",replace
|
|
|
if "`d'"=="" {
|
|
|
tempname d
|
|
|
matrix `d'=[-1\-.5\0\.5\1]
|
|
|
}
|
|
|
|
|
|
local nbitems=rowsof(`d')
|
|
|
|
|
|
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: " _c
|
|
|
tempname dd
|
|
|
matrix `dd'=`d''
|
|
|
matrix list `dd',noblank nohalf nonames noheader
|
|
|
|
|
|
matrix `dd'=`d'/*sqrt(`var')*/
|
|
|
local gamma=`gamma'/*sqrt(`var')*/
|
|
|
|
|
|
if "`data'"=="" {
|
|
|
clear
|
|
|
local temp=2^(`nbitems')
|
|
|
qui range x 0 `=`temp'-1' `temp'
|
|
|
qui g t=x
|
|
|
loc i=1
|
|
|
qui count if t>0
|
|
|
loc z=r(N)
|
|
|
qui while `z'>0 {
|
|
|
qui g item`i'=mod(t,2^`i')==2^`=`i'-1'
|
|
|
qui replace t=t-item`i'*2^`=`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''
|
|
|
qui gen mean=-`n1'*`gamma'/(`n0'+`n1') if group==0
|
|
|
qui replace mean=`n0'*`gamma'/(`n0'+`n1') if group==1
|
|
|
|
|
|
if "`fast'"=="" {
|
|
|
qui gen proba=.
|
|
|
forvalues i=1/`=2*`temp'' {
|
|
|
local int=1
|
|
|
forvalues j=1/`nbitems' {
|
|
|
qui su item`j' in `i'
|
|
|
local rep=r(mean)
|
|
|
local diff=`d'[`j',1]
|
|
|
local int "`int'*exp(`rep'*(x-`diff'))/(1+exp(x-`diff'))"
|
|
|
}
|
|
|
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'
|
|
|
}
|
|
|
}
|
|
|
else {
|
|
|
qui gen proba=1
|
|
|
forvalues i=1/`nbitems' {
|
|
|
qui gen eps`i'=exp(mean-`d'[`i',1])
|
|
|
qui replace proba=proba*eps`i'^item`i'/(1+eps`i')
|
|
|
}
|
|
|
}
|
|
|
qui gen eff=.
|
|
|
forvalues i=0/1 {
|
|
|
qui replace eff=proba*`n`i'' if group==`i'
|
|
|
}
|
|
|
qui replace eff=proba
|
|
|
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'
|
|
|
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*
|
|
|
gen res=proba*50
|
|
|
*list item* group efftmp eff2 proba
|
|
|
|
|
|
qui expand eff2
|
|
|
qui drop proba eff eff2
|
|
|
}
|
|
|
qui gen i=_n
|
|
|
|
|
|
*su
|
|
|
qui alpha item*
|
|
|
local alpha=r(alpha)
|
|
|
|
|
|
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 groupc=group-.5
|
|
|
matrix est=(`gamma',`=sqrt(`var')')
|
|
|
|
|
|
|
|
|
|
|
|
if "`gammafix'"=="" {
|
|
|
constraint 1 _cons=`=ln(`var')'
|
|
|
qui xtlogit rep groupc ,nocons i(i) offset(offset) constraint(1)
|
|
|
}
|
|
|
else {
|
|
|
qui gllamm rep groupc, nocons i(i) offset(offset) iterate(0) fam(bin) link(logit) from(est) copy
|
|
|
}
|
|
|
|
|
|
|
|
|
tempname b V
|
|
|
matrix `b'=e(b)
|
|
|
matrix `V'=e(V)
|
|
|
local gammaest=`b'[1,1]/*sqrt(`var')*/
|
|
|
local se=`V'[1,1]^.5/*sqrt(`var')*/
|
|
|
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'/*sqrt(`var')*//`se')+normal(-1.96-`gamma'*sqrt(`var')/`se')
|
|
|
local poweruni=1-normal(1.96-`gamma'/*sqrt(`var')*//`se')
|
|
|
local clpower=normal(sqrt(`n1'*`gamma'^2/(`n1'/`n0'+1))-1.96)
|
|
|
/*si on ne n<>glige pas le deuxi<78>me terme, la bonne puissance est*/
|
|
|
*di in green "Estimation of the power" _col(60) in ye %6.4f `power' _col(86) in ye %6.4f `clpower'
|
|
|
di in green "Estimation of the power" _col(60) in ye %6.4f `poweruni' _col(86) in ye %6.4f `clpower'
|
|
|
/*si on ne n<>glige pas le deuci<63>me terme, la bonne puissance est*/
|
|
|
*local clnsn=(`n1'/`n0'+1)/((`n1'/`n0')*(`gamma'/sqrt(`var'))^2)*(1.96-invnorm(1-`power'))^2
|
|
|
local clnsn=(`n1'/`n0'+1)/((`n1'/`n0')*(`gamma'/sqrt(`var'))^2)*(1.96-invnorm(1-`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'
|
|
|
di in green "Ratio of the number of patients" in ye %6.2f _col(68)`=(`n0'+`n1')/(2*`clnsn')'
|
|
|
if `expectedpower'!=-1 {
|
|
|
qui sampsi `=-`gamma'/2' `=`gamma'/2', sd1(`=sqrt(`var')') sd2(`=sqrt(`var')') alpha(0.05) power(`expectedpower')
|
|
|
local expn=r(N_1)
|
|
|
local expn2=`expn'*`=(`n0'+`n1')/(2*`clnsn')'
|
|
|
di in green "Number of patients for a power of" %6.2f `=`expectedpower'*100' "%" _col(51) in ye %7.2f `expn2' "/" %7.2f `expn2' _col(77) in ye %7.2f `expn' "/" %7.2f `expn'
|
|
|
}
|
|
|
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=`=`n0'/`clnsn''
|
|
|
return scalar CronbachAlpha=`alpha'
|
|
|
|
|
|
|
|
|
|
|
|
capture qui use `raschpowerfile',clear
|
|
|
|
|
|
end
|