corentinchoisy
/
simul_confusion
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
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.
master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'master'
${ noResults }
simul_confusion/Modules/ado/personal/r/raschpower54.ado

821 lines
27 KiB
Plaintext
Raw Permalink Normal View History Unescape Escape

Initial commit
7 months ago
*! version 5.4 : March 28th, 2014
*! 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 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) : inequal variance between groups
* version 5.2 : January 22th, 2014 (Jean-Benoit Hardouin, Myriam Blanchin) : inequal 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
*
* 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
* http://www.sphere-nantes.org
*
* News about this program : http://www.anaqol.org
* FreeIRT Project : http://www.freeirt.org
*
* Copyright 2010-2014 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 raschpower54,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)]
version 11
if "`html'" != "" {
di "<!-- SphereCalc start of response -->"
di "<pre>"
}
tempfile raschpowerfile
capture qui save "`raschpowerfile'",replace
tempname db d dlong matvar
/*******************************************************************************
DEFINITION DES PARAMETRES
*******************************************************************************/
if "`graph'" != "" {
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 in red "{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'
}
/*******************************************************************************
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 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<6C>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'/(`n1'+`n0')') deltagroup(`gamma') nbobs(1000000)
qui simirt, clear pcm(`d') cov(`var1') mu(`=`n0'/(`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'/(`n1'+`n0')') nbobs(`=round(`n0'/(`n1'+`n0')*1000000)')
qui gen group=0
qui append using `save1'
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'"=="" {
if "`dixj'"!="10" {
di ".." _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)
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 {
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'"=="" {
if "`dixj'"!="10" {
di ".." _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)
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
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
gsort item*
qui expand eff2
}
qui drop proba eff eff2
}
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 `var0'==`var1' {
if "`freevar'"=="" {
local fv "fixedvar(`var0')"
}
else {
local fv
}
}
else {
if "`freevar'"=="" {
local fv "fixedvargroupc(`var')"
}
else {
local fv "vargroupc"
}
}
if "`freeitems'"=="" {
local fi "fixed(`db') fixedmu"
local row=1
}
else {
local fi
local row=(`nbmodat'-1)*`nbitems'+1
}
*di "pcm item*, `fi' covariates(groupc) `fv'"
noi pcm 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{
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 clpower=normal(sqrt(`n1'*`gamma'^2/((`n1'/`n0'+1)*`var'))-1.96)
local clnsn=(`n1'/`n0'+1)/((`n1'/`n0')*(`gamma'/sqrt(`var'))^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'' _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'" != "" {
matrix res=J(15,5,.)
matrix colnames res= n0 n1 gamma variance power
local l=1
if "`longitudinal'"==""{
foreach var in 0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 5 6 7 8 9 {
qui raschpower54, 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(0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 5 6 7 8 9) xtitle(Variance of the latent variable) 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 {
twoway (connected power variance,sort), title("n0=`n0', n1=`n1', {&gamma}=`gamma'") ylabel(0(0.2)1) xlabel(0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 5 6 7 8 9) xtitle(Variance of the latent variable) ytitle(Power) name(planif_graph,replace)
}
}
else {
if `gcorr' != -1 {
local l=1
foreach var in 0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 5 6 7 8 9 {
local cov = `gcorr'*`var'
local matv "`var',`cov'" " \ `cov',`var'"
qui raschpower54, 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(0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 5 6 7 8 9) xtitle(Variance of the latent variable) 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 {
twoway (connected power variance,sort), title("n=`n0', {&gamma}=`gamma', corr=`gcorr'") ylabel(0(0.2)1) xlabel(0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 5 6 7 8 9) xtitle(Variance of the latent variable) ytitle(Power) name(planif_graph,replace)
}
}
if `gvar' != -1 {
matrix res=J(15,5,.)
matrix colnames res= n0 n1 gamma corr power
local l=1
forvalues corr = 0.1(0.1)0.9 {
local cov = `corr'*`gvar'
local matv "`gvar',`cov'" " \ `cov',`gvar'"
qui raschpower54, 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 variance,sort), title("n=`n0', {&gamma}=`gamma', var=`gvar'") ylabel(0(0.2)1) xlabel(0.4 0.7 0.9) xtitle(Correlation) 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 {
twoway (connected power corr,sort), title("n=`n0', {&gamma}=`gamma', var=`gvar'") ylabel(0(0.2)1) xlabel(0.1[0.1]0.9) xtitle(Correlation) ytitle(Power) name(planif_graph,replace)
}
}
}
//di in gr "Method: " in ye "`method'"
if "`html'" != "" {
di "</pre>"
}
}
capture qui use `raschpowerfile',clear
end