*! version 3.2 : February 6th, 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 * * 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